A MATHEMATICAL MODEL OF STEM CELL REGENERATION WITH EPIGENETIC STATE TRANSITIONS

. In this paper, we study a mathematical model of stem cell regeneration with epigenetic state transitions. In the model, the heterogeneity of stem cells is considered through the epigenetic state of each cell, and each epigenetic state deﬁnes a subpopulation of stem cells. The dynamics of the sub- populations are modeled by a set of ordinary diﬀerential equations in which epigenetic state transition in cell division is given by the transition probability. We present analysis for the existence and linear stability of the equilibrium state. As an example, we apply the model to study the dynamics of state transition in breast cancer stem cells.


1.
Introduction. Stem cells play crucial roles during development, tissue regeneration and healthy homeostasis in a whole life cycle. Stem cells provide regeneration in self-renewing tissues through proliferation, differentiation, and apoptosis. A well controlled population dynamics of stem cells is essential for healthy tissue physiological functions [18,27]. However, despite the long-running investigation of stem cell biology, the mechanisms by which stem cell numbers and activity are regulated are still not completely understood [22].
Stem cell biology is population biology. Many mathematical models of the population dynamics have been widely studied in understanding how stem cell regeneration is modulated in different context [6,10,15,17,24,25,26,27,33,39,43]. In most models, the dynamics of a homogeneous cell pool or the lineage of several homogeneous subpopulations are formulated through a set of differential equations. However, the heterogeneity of stem cells is highlighted in recent years due to novel experimental techniques at single cell level [7,11,12,14,36,41]. The heterogeneity is mostly originated from random changes of epigenetic state at each cell cycle, including DNA methylations, histone modifications, and transcriptions of genes and noncoding RNAs. For heterogeneous populations in which qualitatively different subpopulations of cells coexist and transit to each other, the validity of traditional population models is not clear [23].
Fifty years ago, Till et al. proposed a mathematical model of stem cell proliferation to consider the inherently random dynamics of individual cells based on a stochastic birth-death process [38]. This work opened up a perspective on stem cell 1380 QIAOJUN SITU AND JINZHI LEI biology of stochastic heterogenous dynamics of stem cell behavior [22]. In 2014, Lei et al. offered another approach to this problem [18]. The authors outline a general mathematical framework that applies tools from optimization theory to understand stem cell dynamics. In their model, stem cell numbers are regulated by rates of proliferation, differentiation, and apoptosis that are continually tuned by both genetic and epigenetic feedback mechanisms to maximize population performance. Key to the process is the classification of the stem cell population over a variety of different epigenetic states, and the association of different epigenetic states with proliferation, differentiation, and apoptosis. Both total cell numbers and the distribution of epigenetic states of the population are regulated by system-level feedbacks [22]. By adapting ideas from evolutionary theory and population biology, Lei et al. have investigated the cell population dynamics under various control strategies and the evolution of the optimal strategy through numerical simulations [18]. Nevertheless, many questions related to the dynamical properties of the model questions are not considered. In this paper, we present mathematical analysis on the existence and stability of the steady states, which are essential for the homeostasis of the long-term stem cell regeneration dynamics.

2.1.
Model description. The model studied in this paper is an extension from a discrete dynamical model of heterogeneity of stem-cell regeneration [18], which was established based on the G 0 cell cycle model with epigenetic transition between cell cycles. Stem cells at cell cycling are classified into resting (G 0 ) or proliferating (G 1 , S, G 2 and M ) phases ( Fig. 1) [6,24]. During each cell cycle, a cell in the proliferating phase either undergoes apoptosis with the probability µ or divides into two daughter cells at the end of M phase. During the proliferating phase, each mother cell duplicates its genome and regenerates the epigenetic states therein, including the patterns of DNA methylations and nucleosome histone modifications [29,40]. However, the epigenetic states are not perfectly inherited, which may lead to dynamical heterogeneity in daughter cells distinguished from the state of mother cell [35,36,37]. One mother cell divides into two daughter cells, each of which can be one of the two different types, either a differentiated cell (with the probability κ) or remains a stem cell (with the probability (1 − κ)) [17,28,34]. A cell in resting phase either returns to the proliferating phase with the probability β, or irreversibly removed in a rate γ due to aging, death, or differentiation.

Mathematical formulation.
To establish the formulation of the above model of stem cell regeneration, we introduce Ω = {X i } i≥1 for a space of all possible epigenetic states of resting phase stem cells. Here the epigenetic state of a cell can be any heritable quantities in gene activity whose changes are not caused by changes in DNA sequence, for example, the expression levels, the patten of DNA methylation, the state of histone modifications, or expression levels of cell surface marker genes that distinguish differentiated cells from stem cells. The dynamics of system state is measured by the number of stem cells at time t during the resting phase with different epigenetic state X i , denoted as N (t, X i ). The total number of stem cells in resting phase is then given by (1)  Figure 1. Model illustration. During stem cell regeneration, cells in the resting phase either enter the proliferating phase with a rate β, or be removed from the resting pool with a rate γ. The proliferating cells undergo apoptosis with a probability µ. Each daughter cell generated from mitosis is either a differentiated cell (with a probability κ) or a stem cell (with a probability (1 − κ)).
Here we always assume discrete epigenetic states. Mathematically, it is easy to extend the discussions to the continuous situation by replacing the summation with an integral over all X in Ω. Each cell in the resting phase reenters cell division at a rate of β that is dependent on the total cell population N (t) through various types of cytokines [5,17,27], as well as the epigenetic state X i of the cell. Thus, the proliferation rate has a form β(N, X i ). Heterogeneity in cell apoptosis, aging, and differentiation are considered, so that the probabilities µ, κ and the rate γ are dependent on the epigenetic state, and are denoted as µ(X i ), κ(X i ), and γ(X i ), respectively. Here µ(X i ) gives the probability that a cell with state X i (the state before cell division) goes to apoptosis in the proliferating phase (mainly at S phase), γ(X i )∆t gives the probability that a cell with state X i at resting phase moves out due to some reasons within a small time interval [t, t + ∆t] (for example, due to cell aging, death, or differentiation). Cell differentiation often happens at mitosis due to the asymmetric cell division [16,28]. A cell with state X i at M phase becomes a differentiated cell after cell division with the probability κ(X i ). Finally, the epigenetic state of a proliferating cell undergoes state transition from G 1 to M phase due to stochastic dynamics in DNA methylations and histone modifications [29,40]. Here, we do not model the detailed biochemical reactions in this process, instead we introduce the transition function p(X i , X j ) for the probability that a mother cell of state X j gives a daughter cell of state X i (at the end of mitosis phase). It is obviously At a small time interval [t, t + ∆t], a cell in the resting phase with state X i either reenters the proliferating phase with the probability β(N, X i )∆t, removed with the probability γ(X i )∆t, or remains at the resting phase. During the proliferating phase, a cell that enters G 1 phase with state X j either undergoes apoptosis with the probability µ(X j ), or produces two daughter cells. Each daughter cell gains epigenetic state X i with the probability p(X i , X j ), and undergoes differentiation with the probability κ(X i ). These processes lead to the following iteration equation for the resting-phase cell number (also referred to [18]) which gives the following differential equation for all X i ∈ Ω. Here the coefficient 2 means that each mother cell generates two daughter cells in cell division.
Equation (3) gives the basic dynamical equation of our model of stem cell regeneration with epigenetic transition. In the equation, β(N, X i ) and γ(X i ) are the proliferation and removing rate, respectively, and µ(X i ), κ(X i ), p(X i , Y j ) are probability functions, therefore, we always assume the following: Furthermore, the function β(N, X) is continuous with N ∈ R + .
Remark 1. The cell division process usually takes some time. Thus, a delay τ is often introduced in the modeling of stem cell regeneration [19,24]. In this case, equations in (3) are rewritten as the following delay differential equations

Remark 2. Adding equations (3) with respect to all
where If we omit the heterogeneity so that all functions are independent of X i , and let If the time delay in cell division is considered, we obtain a delay differential equation Thus, we again obtain the ordinary/delay differential equation model of stem cell regeneration [5,17,20,27].
Remark 3. If the epigenetic state space Ω is continuous, we can replace the summation in (3) with an integral, so that the equation is rewritten as (10) Here, an integral term is used when the epigenetic state space is continuous. This idea of using an integral to represent the continuous properties of stem cell dynamics has also been used in previous studies [2,3,4].
Upon the establishment of mathematical models, basic dynamical analysis of the model equations are important for our understanding of the system behavior and future extension of the model. In the next section, we give preliminary analysis results for the modeling equation.
3. Analytic results. In this section, we discuss the existence and stability of the equilibrium states of (3). It is easy to see that equation (3) always has a trivial equilibrium state N (t, X i ) ≡ 0. However, this solution is not biologically interesting.
Here we mainly discuss the positive equilibrium state.
In general, equation (3) includes a huge number (as many as possible epigenetic states) of differential equations that are coupled to each other. Thus, it is difficult to analyze the dynamics in general. Here, to make the analysis possible, we introduce a homogenous proliferation assumption: (A1) The proliferation rate β is independent of the epigenetic state of the cell, i.e., β(N, X) = β(N ). Biologically, cell division is mainly triggered by the cytokines secreted from other cells in the niche.
Theorem 3.1. Consider equation (11) and the general assumption (4), if the following conditions are satisfied then there exists at least one positive equilibrium state.
Proof. Let N * = (N * (X 1 ), · · · , N * (X n )) T and Multiplying the i'th equation in (13) Thus, introducing a coefficient matrix yields the following linear equation Now, equation (15) indicates that N * is an eigenvector of A corresponding to the eigenvalue 1 β(N0) . Thus, (15) has a positive solution if and only if A has an eigenvalue λ * ∈ [ 1 β max , 1 β min ], and this eigenvalue has a non-negative eigenvector z = (z 1 , · · · , z n ) T . In this case, the equilibrium state is given by Next, we show that conditions (1)-(3) ensure the existence of the eigenvalue λ * . First, when q ii ≥ 1 2 , A is a non-negative matrix. Moreover, since the transition matrix P is irreducible, A is irreducible as well 1 . Thus, we have the following results according to Perron-Frobenius theorem: 1. The spectral radius of A, denoted by λ * , is a positive real number and an eigenvalue of A. 2. λ * is bounded by the sum of row elements of A, i.e., 3. λ * is a simple eigenvalue and the corresponding unique eigenvector z = (z 1 , · · · , z n ) T is positive (i.e. z i > 0). 4. λ * is the only eigenvalue of A with a non-negative eigenvector.
Finally, letting it is easy to see that N * (X i ) (i = 1, · · · , n) gives a non-negative solution to equation (15).

Remark 4.
Biologically, q ij gives the probability that a mother cell of state X j gives a daughter cell of state X i after one cell cycle. Hence, the condition q ii ≥ 1 2 means that each cell has a probability not less than 1 2 to generate a daughter cell with the same epigenetic state. In particular, if we omit the heterogeneity, then q = (1 − µ)(1 − κ). In this case, the condition q ≥ 1 2 is necessary to have a positive steady state.
Thus, the two matrices A and P share the same indexes of non-zero non-diagonal elements. Hence, A is irreducible if and only if P is irreducible.
proliferation rate β(N ) is a monotone function (which is true in many cases), the corresponding solution(the positive equilibrium state) is unique.
Remark 6. If the transition matrix P is reducible, we can separate the system into several irreducible subsystems and rewrite the model equation in hierarchical structure(refer [21] for a computational method to determine the reducibility of the matrix). Now, to further simplify the equations, we introduce additional assumptions below: (A2) During cell cycle, cell apoptosis is mainly caused by serious DNA damage, which is mostly dependent on the number of initial damage sites and the activity of the DNA damage response pathways [30,31,42]. Here we assumed that the probability of apoptosis µ is independent of the epigenetic state under consideration. (A3) The transition probability p(X, Y ) is normally dependent on both epigenetic state of the mother cell and the daughter cell. Here, for simplicity, we assumed that the transition process approaches the statistical equilibrium state so that the probability p(X, Y ) is independent to the state of mother cell, i.e., p(X, Y ) ≡ p(X). Biologically, this condition means that the proliferating phase is long enough so that the epigenetic dynamics(DNA methylation and histone modifications) can reach the stationary state during proliferation. With assumptions (A1)-(A3), we further rewrite the equation (3) as The following theorem gives a necessary and sufficient condition for the existence of a positive equilibrium state of the equation (18).
Proof. At the equilibrium state, equation (18) gives Hence, we obtain Substituting (25) into (24), we have Thus, from (19), equation (26) has at least one positive solution if and only if (20) is satisfied, and when (20) is satisfied, the equilibrium is given by (25). If β(N ) is a monotone function, equation (26) has a unique solution, and hence the equilibrium state is unique. Remark 7. Theorem 3.2 indicates that the fraction of N * (X i ) at the equilibrium state is given by 3.2. Linear stability of the equilibrium state in two-state systems. Now, we study the linear stability of the equilibrium state under assumptions (A1)-(A3). Hereinafter, we always assume that β(N ) ∈ C 1 (R + , R + ). We first consider the system with two states, and the general situation of n ≥ 2 epigenetic states is discussed in the next section. Proof. When Ω = {X 1 , X 2 }, equation (18) becomes The linearization near the equilibrium state N * = (N * (X 1 ), N * (X 2 )) T is given by where A = (a ij ) is the Jacobi matrix Here we write β, β , N for β(N (t)), β (N (t)), N (t), respectively, for short. Thus, we obtain Next, we analyze the last term in (34), which determines the sign of Det(A). From Theorem 3.2, and (24)-(26), we have Hence, it is easy to see that Det(A) > 0 if and only if β (N 0 ) < 0. Next, we have It is easy to verify that tr are either all non-negative or all non-positive, then the equilibrium state N 0 is asymptotically stable.
Proof. First, we introduce functions f i of (N (t, X 1 ), · · · , N (t, X n )) as Equations in (18) are rewritten as Let N * = (N * (X 1 ), · · · , N * (X n )) and N 0 = N * 1 be the equilibrium state, i.e., Similar to the calculations in the proof of Theorem 3.3, we obtain where Hence, the coefficient matrix A of the linearization of (35) at the equilibrium state N = N * is given by Therefore, the eigenfunction of the coefficient matrix is given by − · · · − b n λ + a n )(λ + a 1 ) · · · (λ + a n ).
From (41), we have the following results for h(λ): From (25) and (40), First, we assume that b i = 0. Then (46) and (39) together yield the following expression (47) Thus, from (26), we obtain Now, we are ready to prove the main results.
Thus, in all situations, all n roots of the eigenfunction are negative, and the equilibrium is linearly stable.
Remark 8. Since γ(X i ) ≤ 1 and β (N 0 ) < 0, it is easy to verify that is enough to ensure the linear stability of the positive equilibrium. Moreover, if β(N ) is taken as a Hill type function (the function widely used in many stem cell population models [17,24,27]) the condition (49) is equivalent to Interestingly, we note that (51) is always satisfied when n = 1.

Stability of the zero solution.
Here, we study the stability of the zero solution 0 under assumptions (A1)-(A3), and prove that if equation (18) has a positive equilibrium state, the zero solution is unstable.
Remark 9. From the discussions in Theorem 3.2, equation (18) has at least one positive equilibrium state if and only if there exists N 0 > 0 so that F (β(N 0 )) = 1. Hence, Theorem 3.5 suggests that the zero solution is unstable if there is a positive equilibrium state, and is linearly stable if there is no positive equilibrium state (except the critical situation F (β(0)) = 1).

Simulations.
In this section, we present simulations of applying our model equation (3) to stem cell population dynamics.

4.1.
Transition dynamics of stem cell regeneration. First, we apply our model to a virtual tissue dynamics of stem cell regenerations. We assumed that in this virtual tissue the space of epigenetic states is given by Ω Here we write X i = i for simplicity. In simulations, the functions µ(X i ), κ(X i ) and γ(X i ) are taken as The proliferation rate function is given by The transition probability P (X i , X j ) is defined as Here the parameter ν is introduced to represent the effect of mother cells. When ν = 0, the transition probability is independent to the state of the mother cell. When ν is larger, however, only local transition is allowed, i.e., the states of the daughter cells are close to that of their mother cell. We take ν = 0, 20, 200 to examine the effects of mother cell control transitions. Simulation results are shown at Fig. 2. Both cell population and the percentage of different epigenetic states converge to a stable equilibrium state after a long time simulation. We note that when ν > 0, the transition matrix is different from the matrix for the case ν = 0 (independent to the mother cell) near the situation of identical division (X i = X j ). Fig. 2 show that there is only slight changes when ν increases from 0 to 20. This result suggests that the equilibrium state is robust with respect to a small perturbation to the state transition probability. However, when ν further increases to 200, there is an obvious change in the equilibrium state.  [13]. The proportion of distinct cell-state subpopulations remain stable in the cell line (proportion of B, S, and L = 97.3%, 1.9%, and 0.62%, respectively). When the pure subpopulations of cells are isolated and cultured, they can rapidly progress toward equilibrium proportions after 6 days growth in culture [13]. A quantitative Markov model of cell-state interconversion was developed to explain the state transition [13]. Here we applied our state transition model to describe the evolution dynamics based on state transition.  Table 1. Cells of these three states are assumed to have the same proliferation rate that is dependent only on the stem cell population N through function (54). Simulation results are shown at Fig. 3, which are in good agreement with experimental data and prediction based on a Markov chain model in [13]. Interestingly, starting from different initial conditions, all three state subpopulation proportions approach to the same equilibrium state. This result indicates the dynamic equilibrium of stem cell regeneration. Markers are data taken from [13]. Parameters are listed in Table 1.

5.
Discussions. Stem cell regeneration is essential during development and the maintaining of homeostasis. Epigenetic state transition is inherent to cell cycling due to the random dynamics of biochemical reactions involved in histone modification and DNA methylation during cell division. Here we have described a general form of continuous dynamical model of stem cell regeneration with epigenetic state transition. In the model, an individual cell is represented by its epigenetic state, and each cell can give two daughter cells with alterations in the epigenetic state. During cell division, the death rate and the probabilities of cell differentiation are assumed to be dependent on the epigenetic state. We give mathematical analysis to the model equations, basic dynamical properties, including existence and linear stability of the equilibrium state, are discussed under general assumptions. The current study is intended to bring a general consideration of the population dynamics of stem cell regeneration with epigenetic transition. Therefore, no molecular or mechanistic details were involved on our model. This leaves a wide space for the extension of our study. To investigate specific functions of a particular type of stem cells, one can add an additional layer of complexity into the model by incorporating corresponding genetic and molecular regulations into the equation of proliferation, differentiation, and apoptosis. In additional to the biological problems, further mathematical questions can be arisen when such details are added to the model, for example, the population dynamics with feedback regulations to stem cell differentiation [17], the signaling pathway to control to oscillatory dynamics in developing hair follicles [1,32], the complex dynamics in hematopoiesis and dynamical blood disease [8,9,44]. This study offers a new approach and mathematical framework to these issues of stem cell biology.