Mathematical Model of Stem Cell Differentiation and Tissue Regeneration with Stochastic Noise

Differentiation and self-renewal of stem cells is an essential process for the maintenance of tissue composition. The promise of novel medical therapies combined with the complexity of this process encourage us to employ numerical and mathematical methods. This will allow us to understand better the mechanisms which regulate stem cell behaviour. Perturbations to the cellular environment may have an influence on the death rate, proliferation rate and on the fraction of self-renewal at every stage of differentiation. In this paper, we present mathematical study of the effect of stochastic noise on the process of tissue regeneration. Here, a system of Itô stochastic differential equations with linear diffusion coefficients that is based on a deterministic model of multistage cell lineages is investigated. Numerical simulations of the stochastic model are shown for a different number of stages of differentiation. Interactions between the noise, added to the different stages, are characterised using numerical simulations. The long-time behaviour of the two-dimensional version of the model is fully characterised; asymptotic stability of the related Markov semigroup is proved using the theory of the Markov semigroups and the method of the Khasminskií function.


Introduction
Tissues are composed of a large number of cells with different functions and morphology. However, most of mature cells are not able to proliferate. Instead, mature cells are replenished by the presence of population of adult stem cells which possess P. R. Paździorek (B) Institute of Mathematics of the Polish Academy of Sciences, Warsaw, Poland e-mail: p.pazdziorek@gmail.com the ability to self renew and to differentiate into various cell lineages. The key role of stem cells is to supply tissue with new cells in homeostatic manner.
The process of maintaining the number of mature cells is a well regulated, multistep process, consisting of differentiation of stem cells and their progeny and of the selfrenewal of the stem cell population. The process of differentiation may be divided into several stages. Both stem cells and their progeny can self-renew (Marciniak-Czochra et al. 2009). However, when the progeny of the stem cells enter to the higher stages of differentiation, their ability to differentiate is reduced and ultimately lost when they reach full maturity.
Mechanisms of regulation of these complex processes are not well understood. The promise of novel stem cells-based therapies, such as therapies for impaired organs, degenerative diseases (Ehnert et al. 2009;Gratwohl and Baldomero 2009;Macarthur et al. 2009) or reconstitution of blood structure after chemotherapy in treatment of leukemias (Marciniak-Czochra and Stiehl 2012), has led to increase of interest in this field. Several mathematical and numerical models were developed to help in understanding stem cell differentiation (Stiehl and Marciniak-Czochra 2010;Marciniak-Czochra et al. 2009;Till et al. 1964;Wu et al. 2009;Colijin and Mackey 2005).
The aforementioned models present different mathematical approaches to describe processes of differentiation and self renewal (Stiehl and Marciniak-Czochra 2010;Marciniak-Czochra et al. 2009;Doumic et al. 2011) involving stem cells proliferation (Till et al. 1964) and mechanisms of stem cell fate decision (Wu et al. 2009) or stem cells regulation system (Colijin and Mackey 2005). Each of these models describes different parts of the process of homeostasis in adult tissue. However, an important difference is how an environmental and intracellular perturbations during those processes are considered, and if they have a reflection in the mathematical form of the models.
The main purpose of this study is to investigate a stochastic stability of the system presented in Marciniak-Czochra et al. (2009) and study how the model responds to noise. We choose this particular model, for its simplicity and for its general macroscopic approach to the problem. An additional advantage of this model is a promising application in describing blood reconstitution after chemotherapy (Marciniak-Czochra and Stiehl 2012).
The paper is organised as follows. First, we present a short description of the model developed in Marciniak-Czochra et al. (2009) and also the biological and mathematical assumptions related to this model. Then, we propose a stochastic version of this model using a system of stochastic ordinary differential equations that are in the form of the original model. Further, we present some numerical simulations of the stochastic model. In the Appendix, we prove mathematical proof of the long-time behaviour of the distribution of the stochastic process described by the stochastic version of the model. We study the evolution of the distribution of the solution of stochastic version of the model which is related to some Markov semigroup. We also investigate the existence of the stationary distribution of the stochastic model. Finally, the asymptotic stability of the Markov semigroup related to the stochastic model is shown. The proof is based on the theory of Markov semigroups and the method of the Khasminskií function (Rudnicki et al. 2002;Pichór and Rudnicki 1997;Skwara 2010a,b).

Deterministic Model
This section is devoted to the deterministic model proposed in Marciniak-Czochra et al. (2009) and analysed in Stiehl and Marciniak-Czochra (2010), Nakata et al. (2012). We assume that there exists a discrete chain of differentiation stages and skipping of some stages during the process of differentiation is impossible. Cell behaviour at each stage of maturation is described by parameters of death rate, proliferation rate and a probability of differentiation. Additionally, it is assumed that the system is regulated by a single cytokine similarly as the red blood cells production is controlled by erythropoietin (Fried 2009;Metcalf 2008;Lasota et al. 1981;Adamson 1994;Ratajczak et al. 1997) or the process of specializing of granulocytes is by G-CSF (Metcalf 2008;Semerad et al. 2002;Price et al. 1996).
Let c 1 , . . . , c n describe the density of population of stem cells, cells at ith stage of differentiation and mature cells, respectively. According to the papers (Stiehl and Marciniak-Czochra 2010;Marciniak-Czochra et al. 2009;Marciniak-Czochra and Stiehl 2012) the phrase 'density of the population' corresponds to the size of the population. Let s be the concentration of the signalling molecules. Following Marciniak-Czochra et al. (2009) it is assumed that the concentration s depends only on the density of mature cells c n . Assuming that the dynamics of the cytokine is a fast process and where μ is the death rate of the mature cells population, and k is a positive constant. Using a quasi-steady state approximation yields to an algebraic dependence of the form described by . (2) Assuming that both proliferation and differentiation processes are regulated leads to the following system of differential equation describing dynamics of n cell subpopulations where p i for i = 1, . . . , n − 1 is the proliferation rate of the population at stage i, μ i for i = 1, . . . , n denotes the death rate at the stage i and a i for i = 1, . . . , n − 1 is the maximal fraction of self-renewal at the stage i. The real fraction of self-renewal at time t is then given by a i s k 1 (t) and defined as a fraction of the direct progeny of cells at stage i which are at the same differentiation stage as their progenitors. The following assumptions are made t ∈ [0, ∞), The death rate, proliferation rate and initial conditions are non-negative and the fraction of self-renewal is between 0 and 1, which corresponds to different types of differentiation: symmetric self-renewal, symmetric differentiation and asymmetric divisions. The model is based on the assumption that cells divide at the rate p i s k2 (t), which results into p i s k2 c i (t) of descendant cells in a unit time t and stage i = 1, . . . , n. The fraction a i of progeny cells remains at the same stage of differentiation as the parent cell, while 1 − a i fraction of the progeny cells differentiates, i.e. transfers to the higher differentiation stage. Additionally, cell death at the rate μ i is modelled.
We also note that when the population of mature cells c n reaches some value, then the term (2a i s k 1 (t) − 1) p i s k 2 (t)c i (t) becomes negative and the number of cells at stage i decreases. On the other hand when the density of the mature cells is low, then is positive and the number of cells at stage i increases providing that the death rates are not too high. It shows how the dynamics of each cell subpopulation depends on the level of mature cells.
The model (3) is well posed, the solution exists and it is unique for t ∈ [0, ∞) and for non-negative initial condition, the solution of system (3) remains non-negative, which is proved in Stiehl and Marciniak-Czochra (2010). Additionally, assuming that it is shown that system (3) has a unique positive steady-state (Stiehl and Marciniak-Czochra 2010). The first inequality provides that there exists a level of the density of mature cells such that c 1 > 0. In other words, the population of stem cells does not simply decrease and becomes extinct but replenishes itself. The second inequality of (5) says that the signal intensity for self-maintenance of the population of stem cells is lower than the concentration needed at some stage of maturation i to maintain the population without influx from the stage i − 1. This interpretation is straightforward if we set d 1 = · · · = d n−1 = 0 and d n > 0. Then, we obtain that condition (5) is equivalent to The latter assumption on the death rates might be applied to cell systems, such as granulopoietic system (Stiehl and Marciniak-Czochra 2010).
In previous works, different numbers of compartments n were chosen. As proposed in Marciniak-Czochra and Stiehl (2012) n = 8 is a direct reflection of the biological case of differentiation of white blood stem cells into Neutrophil Granulocyte (Jandl 1996). The number n = 6 considered in Marciniak-Czochra et al. (2009) corresponds to long term repopulating stem cells, short term repopulating stem cells, multipotent progenitor cells, committed progenitor cells, precursors and mature cells. We may also consider model simplifications with n = 2 or 3. In the first case we treat all immature cells jointly with stem cells as one population and mature cells as another one. The second case is related to the properties of dividing cells. The first compartment might be understood as stem cells and their direct progeny, which are characterised by a low frequency of divisions, the second one as the direct progenitors of mature cells, which are characterised by a high frequency of division. The third compartment is treated as a population of mature cells, which do not divide any more.
In the rest of this chapter we will consider the following 2-dimensional version of the model (7): where as was mentioned in the previous section, processes ξ 1 , may be interpreted as a density of population of cells which are not yet differentiated and ξ 2 as a density of mature cells population. The latter model is well posed, the solution exists and is unique, which results directly from the following consideration. System (8) is understood as a solution of the following stochastic integral equation One may note that the following formulas are bounded Using comparison theorem for one-dimensional Itó process (Ikeda and Watanabe 1981, p. 352) and (10) we may find two stochastic processes (N 1 (t), N 2 (t)), (M 1 (t), M 2 (t)) described by the following systems of Itô stochastic differential equations and satisfying the following inequalities Since a solution of the system (11)-(12) exists and is unique, then the process described in Eq. (9) exists and is unique. It is straightforward that the processes ξ 1 , ξ 2 are nonnegative if and only if the initial processes are non-negative. Since the solution of the Eqs. (11) is as follows and M 1 (0), M 2 (0) ≥ 0, the processes ξ 1 (t), ξ 2 (t) are non-negative.

Biological Interpretation of the Noise in the Model
There are different ways to introduce stochastic perturbations into a continuous model. In some cases, stochastic perturbations might be computed as a limit of the discrete random process. In population dynamics, a random birth-death process is usually the base of such consideration. On the other hand, a non-degenerated diffusion in a system of stochastic differential equations, which is proportional to the variables, may be understood as an environmental stochastic perturbation. Similarly to papers (Allen 2003;Arnold et al. 1979;Dalal et al. 2008) a stochastic noise in the model (8) is understood as an environmental perturbation which act independently and proportionally on each population. Noise in the model (8) may be introduced as a limit of birth-death process but in our model it would be difficult to justify it. It is related to the macroscopic response of the population for environmental perturbation. Additionally, the noise in the first equation of system (8) might be understood as an environmental stochastic noise that has an impact on the death and proliferation process. We also assume that this noise has no impact on the regulation of the differentiation process. Independence of the noise in system (8) is understood as a lack of correlation between Wiener processes W 1,t , W 2,t . The case that the processes W 1,t , W 2,t are correlated is much more complicated and needs different mathematical methods to be investigated hence, we omit this case in this paper.
Proposition 1 Let (ξ 1 (t), ξ 2 (t)) be a solution of the system (8) and assume that (8). If it is not satisfied, the process ξ 1 in the model (8) tends to zero almost sure. And from the fact that the growth of the process ξ 2 depends only of the level of the process ξ 1 , we may conclude that the process ξ 2 (t) tends to zero almost sure. Analytical consideration of this fact and the proof of the proposition 1 are presented in Appendix. The proposition 1 shows an important fact that under the condition (2a − 1) p − α 2 1 2 > 0, the distribution of the populations stabilises and tends to the stationary distribution.
The simulations are done using the Matlab R2007b, and SDE Toolbox created by Picchini (2007). The stochastic trajectories were approximated using Euler-Maruyama scheme (Picchini 2007) while the approximation of mean values is done using Monte-Carlo simulations.
All the simulations present mean trajectory and quartiles of 150 realisations of solutions of the model (7)  We present three numerical simulations of the stochastic model with n = 8 for three different sets of diffusion coefficients. Figure 1 depicts the simulation with small diffusion coefficients for all population. It is well known that the realisation of the process with small diffusion coefficients is similar to the solution of deterministic model (Freȏdlin and Wentzell 1998).  (7) with large diffusion coefficients for subpopulations where a 1 = 0.735, a 2 = 0.7298, a 3 = 0.7245, a 4 = 0.7140, a 5 = 0.5775, a 6 = 0.4725, a 7 = 0.3675, p 1 = 0.006, p 2 = 0.03, p 3 = 0.18, p 4 = 0.6, p 5 = 0.65, p 6 = 1, p 7 = 1.5, μ 1 = μ 2 = μ 3 = μ 4 = μ 5 = μ 6 = μ 7 = 0, μ 8 = 2.77, α 1 = 0.01, α 2 = 0.012, α 3 = 0.018, α 4 = 0.03, α 5 = 0.18, α 6 = 0.5, α 7 = 0.6, α 8 = 1, k = 1.28 × 10 −9 In simulations presented in Fig. 2, we increase the diffusion coefficients for the stem cells population and populations at first two stages of maturation. A large noise of first three populations leads to the extinction. In this description by 'extinction' we mean a convergence of the expected value to zero. However, the interesting observation is that the noise added to the first stages of differentiation does not have any influence on the population of mature cells and the populations at late stages of differentiation. On the other hand, Fig. 3 presents the simulations where we increased the noise for the population of mature cells and the three last stages of differentiation. One may notice that the effect seen in Fig. 2 is obtained in Fig. 3. While the addition of noise disturbs the processes at last stages of differentiation, the first stages of differentiation seem to be insensitive for such perturbation. We might form a conclusion of those three simulations that the effect of noise directed to some subpopulation is dumped at the sufficiently distant subpopulation. Figure 4 presents a realisation of the 2-dimensional deterministic and stochastic version of the model (7). It is seen that in the Fig. 4b we present the case where perturbations applied to the subpopulations of undifferentiated cells and mature cells are relatively small consequently, the mean value of 150 trajectories of the stochastic process is almost indifferent to the numerical solution of the equivalent deterministic model shown in the Fig. 4a. In the Fig. 4c one can notice that the level of the noise that is required to lead most paths to zero is relatively large and even in this figure there are still many of trajectories that reach a high level. One can observe that even in this case the mean value of the simulated trajectories stabilises at some time t, which is a consequence of the asymptotic stability of the related Markov process.

Conclusions
Simulations presented in Figs. 1, 2 and 3 show the behaviour of the stochastic model (7) and these need to be interpreted together with the form of the deterministic model. In Fig. 2, we see the extinction of the first three subpopulations. Provided that a 4 > 1 2 , the fourth subpopulation should take the role of the stem cells and maintain the whole process. On the other hand, since the third quartiles of the first three subpopulations reach a high level in comparison to the unperturbed case, some of the trajectories also reach large values. However, such a case is not observed in subsequent subpopulations. In Fig. 3 we may observe similar effects. The four last subpopulations are perturbed, although, mean values and quartiles of the trajectories of the first three subpopulations are similar to those presented in Fig. 1. Since some of the realisations of the last subpopulation reach very low levels, the negative feedback should lead to an increase of the stem cell population; however, this is not observed in the stochastic simulation. In Fig. 3, we note that some of the trajectories of the fourth subpopulation reach twice level observed in Fig. 1. Therefore, in the case when the last subpopulatons are perturbed, subpopulations which are the nearest to the maturation level are mobilized to stabilize the process.

Basic Definitions
Let (S, , m) be a σ -finite measure space and let D ⊂ L 1 = L 1 (S, , m) be the set of densities i.e.
for almost all y ∈ S and for every f ∈ D. We call function k a kernel of the operator P.
Definition 4 A density f * is called invariant if P(t) f * = f * for each t > 0.

Definition 5 The Markov semigroup {P(t)} t≥0 is called asymptotically stable if there is an invariant density
Remark 6 Let {P(t)} t≥0 be an integral Markov semigroup generated by the equation where A is some differential operator. If there exists a unique stationary solution u * ∈ D of (23), u * > 0, then {P(t)} t≥0 is asymptotically stable which is equivalent to

Definition 7 A Markov semigroup {P(t)} t≥0 is called sweeping with respect to a set
Theorem 11 (Rudnicki et al. 2002, p. 10 Assume that the semigroup {P(t)} t≥0 is generated by the following Fokker-Planck equation where matrix a = [a i j ] is symmetric and positive definite. It means that d i j=1 For the following theorem let A be a differential operator such as on the right-hand side of the Eq. (29).
Theorem 12 (Pichór and Rudnicki 1997, p. 59) Let {P(t)} t≥0 be the Markov semigroup and let A * be a adjoint operator of the infinitesimal generator A of this semigroup. Assume that there exist a C 2 function V : S → [0, ∞), r > 0 and ε > 0 such that

Then{P(t)} t≥0 is asymptotically stable and V is a Khasminskií function for the Markov semigroup {P(t)} t≥0 .
Approximation Theorem for One Dimensional Itô stochastic Differential Equation Let μ, σ be Lipschitz functions and assume that σ (x) > 0. Let X (t) be the solution of the following SDE with For every x let If I 1 (x 0 ) = ∞ and I 2 (x 0 ) = ∞, then If I 1 (x 0 ) < ∞ and I 1 (x 0 ) < ∞, then If An invariant density f * of the semigroup {P(t)} t≥0 related to the Eq. (32) should satisfy the following equation and invariant density f * exist if and only if Additionally, we may obtain that for C > 0. We consider the first equation from the system of stochastic differential Eq. (12) where a > 1 2 , p, α > 0 and N (0) = n 0 . An invariant density for the diffusion N (t) does not exist what is seen by To investigate the behaviour of trajectories, we need to notice that N (t) > 0 a.s. for n 0 > 0. Hence, we may write the following integrals One may notice that for (2a − 1) p − α 2 2 < 0 we get that I 1 (n 0 ) < ∞ and I 2 (n 0 ) = ∞. Therefore, from (35), we obtain that and for (2a − 1) p − α 2 2 > 0 we obtain that I 2 (x) < ∞ and I 1 (n 0 ) = ∞ so from (36) As is shown in Sect. 3 the stochastic process ξ 1 (t) is bounded from above by the process N (t) a.s. Therefore, for (2a − 1) p − α 2 2 < 0 we obtain that Proof of the Proposition 1 To show that the statement of the Proposition 1 is true, firstly, we show that the density of the solution of the system (9) is related to an integral Markov semigroup {P(t)} t≥0 . From the theorem 9 we know that the integral semigroup satisfies Foguel alternative. Therefore, semigroup {P(t)} t≥0 is asymptotically stable or sweeping from the family of all compact sets. In order to exclude sweeping we construct a function V such that the condition (31) from the proposition 1 is satisfied.
where m is a Lebesgue measure on R 2 and B(R 2 ) is a σ -field of Borel sets in R 2 . We consider the system of differential equations (8), where c 1 (t), c 2 (t) are stochastic processes adapted to the natural filtration. The solution of this system is non-negative providing the initial processes are non-negative. We may make the following technical substitution ξ 1 = e x , ξ 2 = e y and consider the following system of differential equations where the process (x(t), y(t)) ∈ R 2 . The probability density of the transition probability function of the process (x(t), y(t)) is described by the following Fokker-Planck equation where g 1 , g 2 : R 2 → R are defined as follows Let A be the differential operator defined as Then, the adjoint operator of operator A is given by Let {P(t)} t≥0 be the Markov semigroup related to the Eq. (47) in the following way where v 0 (x, y) is the initial probability density of the process (x(t), y(t)). Additionally, the operator A defined in (49) is an infinitesimal generator of the semigroup The following remark shows that the semigroup {P(t)} t≥0 generated by the model (8) is integral.
Remark 13 If the Wiener processes W 1,t , W 2,t are not correlated, then the operator A is elliptic and for every (x 0 , y 0 ) the transition probability function P(t, x 0 , y 0 , B) of model (46) is absolutely continuous with the respect to the Lebesgue measure and has a density ν(t, x 0 , y 0 , x, y) ∈ C 2 (0, ∞) × R 2 × R 2 . In addition, (52) Hence the semigroup {P(t)} t≥0 is integral.
Therefore, by Theorem 9 the Markov semigroup {P(t)} t≥0 satisfies the Foguel alternative. According to the Theorem 12, in order to exclude sweeping we need to construct the Khasminskií function for the semigroup {P(t)} t≥0 . To simplify the form of the proof, we firstly show the general way of finding Khasminskií function.
Level lines of a Khasminskií function V (x, y) for sufficiently large R and x 2 + y 2 > R should consist of line segments and circle segments. A construction should be done in the following way. The vector field (g 1 (x, y), g 2 (x, y)) at the level line V (x, y) = z is directed inside the domain = {(x, y): V (x, y) ≤ z}, see Fig. 5 for graphical interpretation. Since the vector field is directed inside the domain we may notice that at the line segments of the Khasminskií function, for some b 0 > 0 the following condition is true At the circle segments we may easily check that Therefore, because α 1 , α 2 are constant, then for sufficiently large r we may obtain that Hence, from (54), (55) we conclude that there exists such b 0 > ε > 0 that the following condition holds for some r > 0 and (x, y) / ∈ B(r ). Since ε > 0 the condtion from the Theorem 12 holds. To finish the proof of asymptotic stability of the Markov semigroup {P(t)} t≥0 we need function V to be C 2 -function. The following technique presents how to obtain such condition.
We convolute our function V with a function f η constructed in the following way. Let f η be a sufficiently smooth function with compact support contained in a closed ball B(η) ⊂ R 2 with centre in (0, 0) and radius η. Let K η be η-neighbourhood of the set K , K 1 , respectively. It means that K η = {x ∈ R 2 : inf u∈K x − u ≤ η}. Hence, let V η = V * f η and then for x / ∈ K η we obtain B(η). In this way we obtain C 2 -function V such that the condition from theorem 12 holds and the proof of the preposition 1 is finished. Therefore, the semigroup {P(t)} t≥0 is asymptotically stable.
Exact construction of Khasminskií function is very technical. We have omitted pure algebraic computation and left only important conditions and inequalities. The proof is divided into two parts. In the first one we assume that In the second one we assume that there exists n ∈ N such that At first we show seweral conditions which holds for model (9). Then we present the exact formula of the function V (x, y) and domains D 1 , . . . , D 12 for related line and circle segments. Consequently we show the boundaries of the derivatives of the function V which lead to estimation of A * V on each domain D 1 , . . . , D 12 . Similar computation is done for the second case. Ultimately we obtain function V for which A * V < −ε for some ε > 0. By virtue of (4), (58) the following inequalities hold for some sufficiently small We fix ε > 0 such that the Eq. (60) hold and define the following constants for r > 0 We may notice that there exists δ > 0 and sufficiently big M > 0 such that the following inequalities hold We note that for every y > x − γ 1 and for every y < x − γ 2 where domains D 1 , . . . , D 12 are defined as follows for some z 0 > 0 We can estimate the first and the second order derivative of the function z = V (x, y) in the domains D 1 , . . . , D 12 and obtain the following inequalities One can also notice that in the domains D 1 , D 7 the following equalities are obtained ∂z ∂ x + ∂z ∂ y = 1 for x, y ∈ D 1 , ∂z ∂ x + ∂z ∂ y = −1 for x, y ∈ D 7 .
In the first part of the proof, under the assumption (58) we may construct the Khasminskií function in the more approachable form, we may obtain the following function and show that for the operator A * defined in (50) there exists ε > 0 such that the following inequality holds for x 2 + y 2 > K where K > 0 is some sufficiently big but finite constant. The latter form of Khasminskií function is easier to generalise and use in the similar forms of the model (8).