Mimic therapeutic actions against keloid by thermostatted kinetic theory methods

This paper deals with the modeling of a wound-healing disease under a therapeutic action by employing the methods of the thermostatted kinetic theory for active particles. In particular, in order to test therapeutic actions against keloid formation and the possible development of a cancer, an external force field coupled to a Gaussian thermostat is introduced into a mathematical model recently proposed. Specifically the model depicts the competition of the immune system cells with a virus, the mutated fibroblast cells, and the cancer cells. Employing a computational analysis, the effects of three different external forces mimic therapeutic actions are analyzed: A vaccine for the virus, the PUVA therapy for the keloid and a vaccine for the cancer. The results are in agreement with the evidence that the sole action of the immune system is not sufficient to obtain a total depletion of keloid thus requiring the definition of a therapy. Further refinements and developments of the model are also discussed in the paper.


Introduction
Wound healing is a dynamic process consisting of four continuous, overlapping, and precisely programmed phases: Hemostasis, the inflammation phase, the proliferation phase, and the remodeling phase, see the review paper [1]. The events of each phase must happen in a precise and regulated manner. Interruptions, aberrancies, or prolongation in the process can lead to delayed wound healing or a non-healing chronic wound, such as keloid. Keloid is a hyperproliferative response of connective tissue in response to skin trauma. The causes which trigger this phenomenon are poorly understood and currently no successful treatment has been developed [2]. In particular the defective control mechanisms in keloid result in the excessive cell proliferation and extracellular matrix synthesis: keloid-derived fibroblasts have a greater proliferative capacity than normal derived fibroblasts [3]. Even if fibroblasts have a major role, other cells like keratinocytes and melanocytes can be involved [2], and the causes can be also the presence of a virus [4] and a generic susceptibility [5]. The possible therapy for keloid strictly depends on the location, size, depth of the lesion and the age of the patient. Therapeutic treatment includes occlusive dressings, compression therapy, intralesional corticosteroid injections, cryosurgery, excision, radiation therapy, laser therapy, interferon therapy, 5-fluorouracil, doxorubicin, bleomycin, verapamil, retinoic acid. Other promising therapies include antiangiogenic factors, including vascular endothelial growth factor inhibitors, phototherapy, tumor necrosis factor (TNF)-alpha inhibitors, and recombinant human interleukin, which are directed at decreasing collagen synthesis.
Recently the defective induction of stress-induced premature senescence during wound healing has been proposed as a possible mechanism of keloid formation [6]. Specifically keloid fibroblasts undergo senescence at a rate lower than that of normal scars, thus depositing collagen and other extracellular matrix proteins beyond that expected in normal wound healing. The proposed mechanism could lead to new treatment possibilities for keloid, e.g. a therapy that induces senescence could be used to prevent the formation of keloid, and consecutively enable the formation of a normal scar.

The thermostatted kinetic framework
This section is concerned with the underlying thermostatted kinetic framework that acts as a general paradigm for the derivation of specific models. Specifically the framework aims at modeling a complex biological system composed of a large number of cells (active particles) that interact in non-linear manner. A constant force field F is assumed to exert an action on the cells.
The microscopic state of a cell is the variable u ∈ D u ⊂ R, which means that the cell at time t ∈ [0, ∞) is able to express a strategy modeled by the variable u (activity variable). Cells expressing the same function are grouped into a subsystem, called functional subsystem. The evolution of each functional subsystem is depicted by the distribution function f i (t, u) : [0, ∞) × D u → R + , for i ∈ {1, 2, . . . , n}, and such that, for any fixed time t, the quantity f i (t, u)du represents the density of cells in the elementary volume du centered at u. In general, the macroscopic variables are defined as moments weighted of the distribution function f i . Specifically the pth-order moment of f i reads where, in particular, the density (mass), the activation momentum (linear momentum), and the activation energy (kinetic energy) of the functional subsystem f i is obtained for p = 0, p = 1 and p = 2, respectively.
. . , f n (t, u)) be the vector, whose components are the distribution functions of the functional subsystems and The evolution equation for the i-th functional subsystem is obtained by equating the time derivative of f i to the balance of the inlet and outlet flows in the elementary volume du. Accordingly we have where F i denotes the external force acting on the cells of the i-th functional subsystem, J i [f ](t, u) denotes the following operator that models the gain-loss of cells due to transitions in the activity variable: where -η ij models the probability that a cell of the i-th functional subsystem, with activity u * , interacts instantaneously with a cell of the jth functional subsystem, with activity u * .
is the density function modeling the probability that cells of the i-th functional subsystem, with activity u * , interacting with cells of the jth functional subsystem, with activity u * , reach the activity u. In particular A ij (u * , u * , u) satisfies the following identity: In what follows, the transport term in (3) that models the introduction of the external force coupled to the Gaussian thermostat [17,18], will be denoted as In particular (5) is a damping operator adjusted to control the global activation energy: Du Remark. If the external force depends on the activity variable, namely F i = F i (u), then the thermostatted term (5) reads The framework (3) can be further generalized by introducing the role of non-conservative processes. Specifically, interactions among the cells may generate proliferation/destruction of other cells (birth-death process). This type of interaction is modeled by the following operator: where μ ij denotes the net proliferative/destructive rate. Cell mutations can occur because of DNA corruptions, thus generating a new functional subsystem. Accordingly, mutations are modeled by the following operator: where ϕ i hk denotes the net mutative rate into the i-th functional subsystem, due to interactions that occur with rate η hk between the cells, with activity u * , of the h-th functional subsystem and the cells, with activity u * , of the k-th functional subsystem.
Bearing all the above in mind, the thermostatted kinetic framework with proliferative/destructive and mutative interactions reads: It is worth stressing that the parameters of the model can be a function of the activity, namely η ij (u * , u * ), μ ij (u, u * ) and ϕ i hk (u * , u * , u). However in order to simplify our model these functions will be assumed as constants. From the mathematical point of view, the Cauchy problem related to the general framework (9) has been analyzed and the existence and uniqueness of the solution has been proved, including the existence of stationary solutions [19]. Therefore we are allowed to develop the appropriate computational methods to obtain simulations of a specific model.

A thermostatted kinetic model for the treatment of keloid
This section is concerned with the derivation of a specific thermostatted kinetic model for the treatment of keloid by means of a therapy. According to the general framework (9), the first step is the definition of the functional subsystems. Following [8], we assume that keloid formation involves the following five interacting functional subsystems (see table 1): 1) Fibroblast cells (Fc): The activity variable represents the proliferation ability and the distribution function is denoted by f 1 (t, u). 2) Virus (V): The activity variable is a magnitude of their aggressiveness related to their proliferation ability. The distribution function is denoted by f 2 (t, u).

3) Keloid-fibroblast cells (KFc):
The activity variable represents the proliferation ability. The keloid-fibroblast cell is a fibroblast cell that, because of a mutation, has acquired a significant advantage with respect to the proliferation rate. The distribution function is denoted by f 3 (t, u).

4) Cancer cells (Cc):
The activity variable is a magnitude of their progression ability and the distribution function is denoted by f 4 (t, u). 5) Immune system cells (ISc): The activity variable u is a magnitude of the activation and thus of the response to foreign agents. The distribution function is denoted by In what follows the cells of the functional subsystems 2, 3, and 4 will be called non-self-cells.
Bearing formula (1) in mind, the local density and the local activation of the i-th functional subsystem, for i ∈ {1, 2, 3, 4, 5}, read, respectively, The term A ij defined in (4) is assumed to be defined by a delta Dirac function (deterministic output d ij (u * , u * ) of a pair interaction) depending on the microscopic state of the interacting pairs. Specifically, let α be a positive parameter which models the heterogeneity rate of the KFc, and ε < 1 a scale parameter that is introduced to evaluate the difference among the heterogeneity rates of the V and the Cc with respect to the KFc, then A ij reads where The assumption (12) means that the functional subsystems V, KFc, and Cc are the only subsystems subject to transitions into the activity variable, namely are derived under the assumption that the non-self-cells have a tendency to increase their microscopic state. Bearing all the above in mind we have Let β be the proliferation rate of the KFc, and β i the proliferation rate of the ISc. The modeling of proliferative events is based on the following assumption for the proliferation rate of cells: According to (14), the V proliferates because of encounters with the Fc and the ISc; the KFc proliferate because of encounters with the Fc and the V; the Cc proliferate because of encounters with the V; the ISc proliferate because of encounters with V, KFc, and Cc. It is worth stressing that, according to (14), the proliferation rate of the KFc is greater than the proliferation rates of Fc, V, and Cc; the proliferation rate of the ISc, when encounter the V and the Cc, is greater than the proliferation rate of the ISc when encounter the KFc. The proliferation term thus writes: Let δ i be the inhibition (destruction) rate of the ISc because of encounter with the V and the Cc, and δ the destruction rate of the V and the Cc because of the action of the ISc. In order to model the interactions that have as result the destruction of cells, we define the following destruction rate: According to (16), the Fc are destroyed because of interactions with the V; the V and the Cc are destroyed because of interactions with the ISc; the KFc are destroyed because of interactions with the V and the ISc; the ISc are destroyed because of interactions with V, KFc, and Cc; the ISc are able to eliminate the V and the Cc more efficiently than the KFc; the latter are destroyed by the V less efficiently than the Fc.
The assumption (16) takes into account that the non-self-cells with a high level of activity have the ability to inhibit or destroy the ISc (immune suppression or immune-subversion). In particular we assume that the Cc and the V have a more efficiently action of inhibition on the ISc with respect to the KFc (inhibited immune cells do not play a relevant role in the competition and may be equivalently assumed as eliminated). Accordingly the destructive terms read: Bearing all the above in mind, we have Let γ be the mutation rate of a fibroblast cell, which measures the possibility that a Fc undergoes a mutation and becomes a keloid-fibroblast cell, and λ the mutation rate that measures the possibility that a keloid-fibroblast cell undergoes a mutation and becomes a cancer cell. The genetic mutations are modeled by defining the following mutative rate: According to (18), we assume that the Fc may mutate into KFc because of encounters with the Fc; the Fc may mutate into KFc because of encounter with the V; the KFc may mutate into Cc because of encounters with the V.
In particular, we have assumed that is more likely that the Fc become KFc when they encounter the V. It is worth noting that (18) implies that the microscopic state of the cells does not change during the mutation. The mutation term thus writes: Finally we assume that a constant external force field F i , for i = {2, 3, 4}, acts on the i-th functional subsystem and mimics a specific therapy. The thermostatted kinetic model, which consists in a system of evolution equations for each distribution function f i , for i ∈ {1, 2, 3, 4, 5}, thus reads: The model (20) is characterized by 11 parameters which are positive constants (eventually equal to zero), small with respect to unity and having the following biological meaning: -α is the heterogeneity rate of the KFc; -β is the proliferation rate of the KFc; -β i is the proliferation rate of the ISc; -δ is the destruction rate of the V and the Cc by the ISc; -δ i is the destruction rate of the ISc by the V and the Cc; -γ is the mutation rate of the Fc into KFc; -λ is mutation rate of the KFc into Cc; -is the scale factor; -F i , for i = {2, 3, 4}, is the external force that acts on the i-th functional subsystem and mimics a therapy.
It is worth pointing out that the α-parameter refers to transition into the activity variable, the β-parameters refer to proliferative events, the-δ parameters refer to destructive interactions, the parameters γ and λ refer to mutations, and ε is a scale parameter, see table 2. All parameters have to be tuned by suitable experiments.

Computational analysis: Mimic therapeutic actions
This section is concerned with the computational analysis of the model (20) and specifically with the evolution of the different functional subsystems when an external action acts on the system as a therapeutic action. The main aim is to simulate the prompt response against the formation and evolution of keloid and the possible onset of cancer. The computational analysis is thus addressed to analyze the effects of three different therapeutic actions: An action which mimics a vaccine against the virus, an action which mimics a vaccine against cancer cells, and an action which can also mimic surgery on keloid. It is important to note that, according to our model, the onset of cancer cells is a consequence of mutations in the keloid cells because of the virus, the latter in part also responsible for mutations in the fibroblasts cells which generate the keloid (remember that according to our assumptions, the genetic susceptibility is also responsible for keloid formation). Therefore the main role of the external actions should be to act on the functional subsystems of virus and keloid cells. However in order to have a more global view of what the model developed in the present paper is able to reproduce, we will consider also an external action on the cancer cells. It is worth stressing that our model is an exploratory model. Thus, at this stage, we are only interested in the emerging phenomena that the model is able to reproduce and not in the tuning with an in vivo or in vitro experiment or empirical data, which can be considered as subject of further investigations. Accordingly, the computational analysis focuses on the model (20) when only one external action is applied. The dynamics of the model when the three actions can act at the same moment will be straightforward. However to think that the three different actions may be applied at the same time cannot be suitable then we believe that our model will fit well if a combination of the three different actions is performed at different times. Specifically the first step is to act against the virus (F 2 = 0, F 3 = F 4 = 0); when the virus is killed, the reached non-equilibrium stationary state will be the initial state for the model (20) with the action of a therapy for keloid (F 3 = 0, F 4 = 0); when keloid cells are eliminated, the new reached non-equilibrium stationary state will be the initial state for the model (20) with the therapeutic action against cancer cells (F 4 = 0). In particular the intermediate step can be performed thanks to the introduction of the thermostat that allows the existence of a non-equilibrium stationary state.
The computational analysis will be addressed by fixing nine of the eleven parameters and performing a sensitivity analysis on the parameter α and on the therapeutical action F i , for i ∈ {2, 3, 4}. In the analysis we will show the evolution of the density, the behavior of the distribution function and the activation energy. The computational scheme is that of the well-known generalized collocation method where the variable u is discretized into a suitable set of collocation points. The integral terms are approximated by means of algebraic weighted sums in the nodal points of the discretization. The particularization of the evolution equations in each node and the enforcing of the initial conditions transform the model that is a system of integro-differential equations into a systems of ordinary differential equations, describing the evolution of the values of the distribution functions in the node of the collocation, see sect. 2 of paper [9] for all the details.
The choice of the distribution functions at time t = 0 is based on the assumption that, before the formation of keloid and the related onset of cancer, the virus infects the Fc. Our analysis starts when the number of Fc in the wound is equal to the number of V and a number of immune system cells have reached the wound. Accordingly we assume non-zero initial conditions only for the functional subsystem of Fc, V, and ISc.
The test case is based on the following choice of the parameters: γ = 0.4 (the rate of mutation of a fibroblast cell into a keloid-fibroblast cell is not negligible), δ = 0.3 (the ability of the immune system cells to inhibit the non-self-cells is quite low), δ i = 0.5 (the non-self-cells have an intermediate ability to inhibit the response of the immune system), β = 0.4 (the non-self-cells have an intermediate ability to proliferate), β i = 0.35 (the rate of proliferation of the ISc is quite low), λ = 0.5 (the rate of mutation of a keloid-fibroblast cell into a cancer cell is not negligible) and = 0.5.
It is worth stressing that the computational analysis performed in the present paper does not cover the whole variety of conceivable dynamics but represents a useful test case.

Simulating the effects of a vaccine against the virus
This subsection deals with the computational analysis for the model (20) when a constant external force, mimic a vaccine against the virus, is introduced. Following the interaction rules proposed into the model, the onset of keloid is a consequence of the virus, then it is expected that if the vaccine is able to reduce the action of the virus then keloid formation and cancer can be prevented. Accordingly we have F 2 = 0, F 3 = F 4 = 0 and we let the magnitude of the parameter α vary from low to higher values, namely from low to higher heterogeneity. The results of the computational analysis are summarized as follows. Analysis for low values of α. The effect of the external action is evident by comparing the dynamics depicted by the model (20) when F 2 = 0 with the dynamics when F 2 = 0.
The dynamics of the model (20) for F 2 = F 3 = F 4 = 0 has been widely debated in [8], therefore the new computational analysis will be focused on the case F 2 = 0 and F 3 = F 4 = 0, and more precisely F 2 ∈ {0, 0.0005, 0.01}. Looking at the left panel of fig. 1, for α = 0.2 (low rate of heterogeneity) the heterogeneity of the virus is bounded, then the main effect of the external force (vaccine for the virus) is to decrease the maximum of the density of V when the magnitude of the vaccine increase. Therefore, thanks also to the action of the immune system, we can see a faster depletion of the V and consequently of the KFc and the Cc. However, as the right panel of fig. 1 shows, the keloid cells start to increase again; this is a consequence of the genetic susceptibility that thus requires the definition of a therapeutical action for keloid. The global energy of the system even in the presence of the external action is preserved thanks to the thermostat. Analysis for intermediate values of α. Setting α = 0.5, the heterogeneity of the non-self-cells becomes less negligible with respect the previous case. Therefore the role of the external force is now fundamental for avoiding keloid formation and onset of cancer. The computational analysis shows that high values of the vaccine are required in order to reduce the number density of keloid cells with high values of activity. Consequently the number density of the Cc is reduced. The vaccine for the virus has an important role in the proliferation of the immune system cells. Indeed as fig. 2 shows, in absence of a vaccine for the V the immune system is inhibited; the vaccine helps the immune system cells to proliferate again. Analysis for high values of α. The heterogeneity of the non-self-cells in now very high thus keloid formation and cancer need to be inspected. As fig. 3 shows the vaccine is able to control the evolution of the virus and the depletion of the virus with high levels of aggressiveness depends on the magnitude of the vaccine. However high values of the vaccine are necessary to eliminate keloid and cancer, see fig. 4, then the definition of a therapeutic action against keloid formation and cancer is now fundamental.

Simulating the effects of a therapy against the keloid
This subsection deals with the computational analysis for the model (20) when a constant external force, mimic a therapeutical action against the keloid, is considered. Following [6], a therapy that induces senescence could be used to prevent the formation of a normal scar. Specifically the treatment proposed in [6] is based on photodynamic and PUVA therapy, which are capable to induce cell senescence. The therapy, based on a combination of a photosensitizer and a light source, produces oxidative stress and thus produce higher quantity of senescent cells with the minor apoptotic and necrotic effect.    Bearing all above in mind, we assume that F 3 = 0, F 2 = F 4 = 0 and again we let the magnitude of the parameter α vary from lower to higher values. Specifically for low values of α, the therapy and the immune system are capable to perform a prompt action, see the left panel of fig. 5, where the density of the KFc is depicted for different values of the external force. Moreover, as the right panel of fig. 5 shows, the global activation energy is controlled by thermostated. The effects of the therapy are visualized in fig. 6, where the distribution function of the KFc is depicted in absence of therapy and in the presence of the therapy. As fig. 6 shows, the therapy allows the depletion of cells at different stages of mutation and a bounded magnitude of therapy is required in order to have the total depletion of keloid. Moreover, differently from the case of a vaccine for the virus, in this case the increase of new keloid cells, after the elimination of keloid, is not observed. For intermediate values of α, the general effect of the therapy is that of acting homogeneously on the KFc and the Cc and thus to obtain a distribution more uniform with respect to the degree of mutation. In particular the Fc are able to proliferate again. For high values of α, the rate of mutation allows the onset of keloid cells with high proliferation rate and of cancer cells with high values of progression. In particular even if the therapy allows the immune system to have a prompt response (see fig. 7) against the cancer cells (and virus), the immune system cells are not able to prevent tumor formation at tissue scale. Therefore a high magnitude of the therapy is required. This is the case where surgery is the only therapy that can work.

Simulating the effects of a vaccine against the cancer
This subsection deals with the computational analysis for the model (20) when a constant external force, mimic a therapeutical action against the cancer, acts on the system. In the case of cancer, the vaccine strictly depends on the type of tumor developed. Therefore the main aim of this subsection is to show that our model is able to consider the introduction of a vaccine and to simulate the evolution of the cancer. In particular we will focus on the case of high values of α. According to our model, the development of cancer is the final result of the keloid formation, then the vaccine against the cancer has a fundamental action only on the Cc and the ISc. As fig. 8 shows, the vaccine is able to inhibit the formation of a tumor at the macroscopic scale. Moreover, the immune system is able to proliferate again, see the right panel of fig. 9. Finally, as the right panel of fig. 9 shows, the action of the therapy allows to maintain bounded the global activation energy of system (thanks to the introduction of the thermostat).

Conclusions and research perspectives
The present paper has been devoted to test the capability of a new thermostatted kinetic framework for the active particles to model a therapy against keloid formation and the possible development of cancer. By employing different therapies on the different functional subsystems that compose the system, we have shown that some therapies proposed in the pertinent literature can be well described by our model. The analysis performed in this paper has been of computational kind and has been focused on the reaching of the emerging phenomena that are typical of keloid formation. The role of the immune system has been taken into account as a whole system without specifying who are the cells involved in the competition; this assumption is based on the fact that its action on keloid is limited by the genetic susceptibility that does not allow to the immune system to recognize these cells as foreign cells and thus acting for a possible total depletion.
The model developed and analyzed in the present paper is based on the introduction of an external force field, which mimics a therapy, at the macroscopic scale but in the absence of external interactions at the microscopic scale. The therapy acts directly on the number of foreign cells without performing an interaction at the microscopic level (in the activity variable). The modeling of external therapy at the microscopic (cellular) scale can be performed by representing the external therapy as a functional subsystem that have the ability to modify the state u of the cell by a particular action related to the variable ω ∈ D u . Assuming that the ith inner functional subsystem interacts with the rth external agent, for r ∈ {1, 2, . . . , m}, and denoting by g ir = g ir (t, ω): [0, ∞) × D u → R + the related distribution function (known function of its arguments), the microscopic external actions are modeled by the following operator where g i = (g i1 , . . . , g im ) and -η e ir is the inner-outer encounter rate between the rth external agent, with state ω * , and the cell of the ith population, with state u * .
-B ir (u * , ω * , u) is the inner-outer transition probability density which describes the probability density that a cell of the i-th population, with state u * , falls into the state u after an interaction with the r-th external agent, whose state is ω * . The density B ir satisfies, for all r ∈ {1, 2, . . . , m} and i ∈ {1, 2, . . . , n}, the following condition: Bearing all the above in mind, the thermostatted kinetic framework for open systems with proliferative/destructive and mutative interactions now reads: It is worth stressing that usually the interaction domain of the cell with state u * is not the whole domain D u but a subset Ω u * ⊆ D u , which contains the cells with activity u * ∈ Ω u * that are able to interact with the particles with activity u * . This is a phenomenon that is typical in tumor dynamics and specifically when the immune system is not capable to interact with all the tumor cells (tumor escape), see the review paper [20]. Specifically the immune response fails to completely eliminate the tumor, and the interact process results in the selection of tumor cell variants that are able to resist, avoid, or suppress the antitumor immune response, leading to the escape phase. Accordingly, a positive function ω(u * , u * ) can be introduced to weight the interactions among the cells; this function is assumed normalized with respect to integration over u * and it has a compact support in the domain of influence Ω u * ⊆ D u of the interactions. Moreover, Du ω(u * , u * ) du * = Ωu * ω(u * , u * ) du * = 1.
The development of a model which includes a therapy against the keloid and that takes into account the above-raised issues is object of future research directions. Moreover the possibility to perform an asymptotic analysis that allows the derivation of the dynamics at the tissue scale is a further research perspective. The asymptotic analysis can be performed by employing the methods developed in papers [21,22]. The final goal is the identification of the parameters with the aim to compare the emerging phenomena with the experimental data.