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 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. Model

---

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 $\mu$ 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 $\kappa$) or remains a stem cell (with the probability $(1-\kappa)$) [17,28,34]. A cell in resting phase either returns to the proliferating phase with the probability $\beta$, or irreversibly removed in a rate $\gamma$ due to aging, death, or differentiation.

---

2.2. Mathematical formulation

To establish the formulation of the above model of stem cell regeneration, we introduce $\Omega = \{X_i\}_{i\geq 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

$\label{eq:Nt} N(t) = \sum\limits_{X_i\in \Omega} N(t, X_i).$

(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 $\Omega$.

Each cell in the resting phase reenters cell division at a rate of $\beta$ 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 $\beta(N, X_i)$. Heterogeneity in cell apoptosis, aging, and differentiation are considered, so that the probabilities $\mu, \kappa$ and the rate $\gamma$ are dependent on the epigenetic state, and are denoted as $\mu(X_i), \kappa(X_i)$, and $\gamma(X_i)$, respectively. Here $\mu(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), $\gamma(X_i) \Delta 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 + \Delta 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 $\kappa(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

$\label{eq:np} \sum\limits_{X_i\in \Omega} p(X_i, X_j) = 1, \forall X_j \in \Omega.$

(2)

At a small time interval $[t, t+ \Delta t]$, a cell in the resting phase with state $X_i$ either reenters the proliferating phase with the probability $\beta(N, X_i)\Delta t$, removed with the probability $\gamma(X_i)\Delta 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 $\mu(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 $\kappa(X_i)$. These processes lead to the following iteration equation for the resting-phase cell number (also referred to [18])

$\begin{align*} N(t+\Delta t, X_i) =& N(t, X_i) - N(t, X_i) (\beta(N, X_i) + \gamma(X_i)) \Delta t\\ & + 2 \sum\limits_{X_j \in \Omega} \beta(N, X_j) N(t, X_j) (1-\mu(X_j)) p(X_i, X_j) (1-\kappa(X_i)) \Delta t, \end{align*}$

which gives the following differential equation

$\label{eq:1} \left\{ \begin{array}{rcl} \displaystyle\dfrac{\partial N(t, X_i)}{\partial t} &=& \displaystyle-N(t, X_i) (\beta(N, X_i) + \gamma(X_i))\\ &&\displaystyle{} + 2 \sum\limits_{X_j\in \Omega} \beta(N, X_j) N(t, X_j) (1-\mu(X_j)) p(X_i, X_j) (1 - \kappa(X_i)), \\ N&=&\displaystyle\sum\limits_{X_i\in\Omega} N(t, X_i), \end{array} \right.$

(3)

for all $X_i\in \Omega$. 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, $\beta(N, X_i)$ and $\gamma(X_i)$ are the proliferation and removing rate, respectively, and $\mu(X_i), \kappa(X_i), p(X_i, Y_j)$ are probability functions, therefore, we always assume the following:

$\label{eq:ge} \beta(N, X_i), \gamma(X_i) > 0, 0\leq \mu(X_i), \kappa(X_i), p(X_i, X_j)\leq 1.$

(4)

Furthermore, the function $\beta(N, X)$ is continuous with $N\in \mathbb{R}^+$.

Remark 1. The cell division process usually takes some time. Thus, a delay $\tau$ 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

$\label{eq:1d} \left\{ \begin{array}{rcl} \displaystyle\dfrac{\partial N(t, X_i)}{\partial t} &=& \displaystyle-N(t, X_i) (\beta(N, X_i) + \gamma(X_i))\\ &&\displaystyle{} + 2 \sum\limits_{X_j\in \Omega} \beta(N_\tau, X_j) N(t-\tau, X_j) (1-\mu(X_j)) p(X_i, X_j) (1 - \kappa(X_i)), \\ N&=&\displaystyle\sum\limits_{X_i\in\Omega} N(t, X_i), \end{array} \right.$

(5)

here $N_\tau = N(t-\tau)$.

Remark 2. Adding equations (3) with respect to all $X_i$ in $\Omega$ gives the equation for $N(t)$:

$\dfrac{d N}{d t} = \sum\limits_{X_i \in \Omega} N(t, X_i) \left[\beta(N, X_i) (2 \sigma(X_i) (1-\mu(X_i)) - 1) - \gamma(X_i)\right],$

(6)

where

$\sigma(X_i) = \sum\limits_{X_j \in \Omega} (1-\kappa(X_j)) p(X_j, X_i).$

(7)

If we omit the heterogeneity so that all functions are independent of $X_i$, and let $\tilde{\beta}(N) = (2 \sigma (1-\mu) - 1) \beta(N)$, $N(t)$ satisfies

$\label{eq:n0} \dfrac{d N}{d t} = \tilde{\beta}(N) N - \gamma N.$

(8)

If the time delay in cell division is considered, we obtain a delay differential equation

$\label{eq:n1} \dfrac{d N}{d t} = \tilde{\beta}(N_\tau) N_\tau - \gamma N.$

(9)

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 $\Omega$ is continuous, we can replace the summation in (3) with an integral, so that the equation is rewritten as

$\label{eq:1i} \left\{ \begin{array}{rcl} \displaystyle\dfrac{\partial N(t, X)}{\partial t} &=& \displaystyle-N(t, X) (\beta(N, X) + \gamma(X))\\ &&\displaystyle{} + 2 \int_{\Omega} \beta(N, Y) N(t, Y) (1-\mu(Y)) p(X, Y) (1 - \kappa(X)) d Y\\ N&=&\displaystyle \int_{\Omega} N(t, X) dX. \end{array} \right.$

(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) \equiv 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 $\beta$ is independent of the epigenetic state of the cell, i.e., $\beta(N, X) = \beta(N)$. Biologically, cell division is mainly triggered by the cytokines secreted from other cells in the niche.

---

3.1. Existence and uniqueness of the positive equilibrium state

Assuming the condition (A1), equation (3) can be rewritten as

$\label{eq:trivial} \left\{ \begin{array}{rcl} \displaystyle\dfrac{\partial N(t, X_i)}{\partial t} &=& \displaystyle-N(t, X_i) (\beta(N) + \gamma(X_i))\\ &&\displaystyle{} + \beta(N)\sum\limits_{X_j\in \Omega}2 (1-\mu(X_j)) p(X_i, X_j)(1 - \kappa(X_i))N(t, X_j), \\ N&=&\displaystyle\sum\limits_{X_i\in\Omega} N(t, X_i). \end{array} \right.$

(11)

In the following discussions, we always assume that there are finite epigenetic states $\Omega = \{X_1, \cdots, X_n\}$, and introduce the following notations:

$\begin{eqnarray*} q_{ij}&=&(1-\mu(X_j)) p(X_i, X_j)(1 - \kappa(X_i)), \\ s_i &=& \sum\limits_{j = 1}^n \frac{2q_{ij}}{\gamma ({X_i})} - \frac{1}{\gamma ({X_i})}, \end{eqnarray*}$

$s = \mathop {\min }\limits_{1 \le i \le n} {s_i}, S = \mathop {\max }\limits_{1 \le i \le n} {s_i},$

$\beta_{\min} = \inf\limits_{N\in \mathbb{R}^+} \beta(N), \beta_{\max} = \sup\limits_{N\in \mathbb{R}^+} \beta(N).$

Here we note that $\beta(N), \gamma(X_i) > 0$ and $0\leq \mu(X_i), \kappa(X_i) \leq 1$ from our general assumption (4).

Theorem 3.1. Consider equation (11) and the general assumption (4), if the following conditions are satisfied

(1) ${q_{ii}}\geq \frac{1}{2}$ for $i=1, \cdots, n$,

(2) $\beta(N)\in C(\mathbb{R}, \mathbb{R}^+)$, and

$\label{eq:sS} \dfrac{1}{\beta_{\max}} < s \leq S < \dfrac{1}{\beta_{\min}},$

(12)

(3) the transition matrix defined with $P=[P(X_i, X_j)]_{i, j=1, \cdots, n}$ is irreducible, then there exists at least one positive equilibrium state.

Proof. Let $\mathbf{N}^*=(N^*(X_1), \cdots, N^*(X_n))^T$ and

$N_0 = \sum\limits_{i=1}^n N^*(X_i)$

denote the equilibrium state. From (11), we have

$-N^*(X_i) (\beta(N_0) + \gamma(X_i)) + \beta(N_0)\sum\limits_{j = 1}^n 2q_{ij}N^*(X_j)=0, \forall i = 1, \cdots, n. \label{etr1}$

(13)

Multiplying the i'th equation in (13) with $\dfrac{1}{\beta(N_0)\gamma(X_i)}$, we obtain

$-N^*(X_i) (\frac{1}{\gamma(X_i)} + \frac{1}{\beta(N_0)}) + \frac{1}{\gamma(X_i)}\sum\limits_{j = 1}^n 2q_{ij}N^*(X_j)=0,$

which gives

$- \frac{1}{\gamma(X_i)} N^*(X_i) + \frac{1}{\gamma(X_i)}\sum\limits_{j = 1}^n 2q_{ij}N^*(X_j)=\frac{1}{\beta(N_0)}N^*(X_i), i=1, \cdots, n. \label{etr2}$

(14)

Thus, introducing a coefficient matrix

$A = \left[ \begin{array}{cccc} \dfrac{2q_{11} - 1}{\gamma ({X_1})}&\dfrac{2q_{12}}{\gamma ({X_1})}& \cdots &\dfrac{2q_{1n}}{\gamma ({X_1})}\\ \dfrac{2q_{21}}{\gamma ({X_2})}&\dfrac{2q_{22} - 1}{\gamma ({X_2})}& \cdots &\dfrac{2q_{2n}}{\gamma ({X_2})}\\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{2q_{n1}}{\gamma ({X_n})}&\dfrac{2q_{n2}}{\gamma ({X_n})}& \cdots &\dfrac{2q_{nn}-1}{\gamma ({X_n})} \end{array} \right],$

(14) yields the following linear equation

$\label{etr3} A \mathbf{N}^* = \frac{1}{\beta ({N_0})} \mathbf{N}^*.$

(15)

Now, equation (15) indicates that $\mathbf{N}^*$ is an eigenvector of $A$ corresponding to the eigenvalue $\frac{1}{\beta(N_0)}$. Thus, (15) has a positive solution if and only if $A$ has an eigenvalue $\lambda^* \in [\dfrac{1}{\beta_{\max}}, \dfrac{1}{\beta_{\min}}]$, and this eigenvalue has a non-negative eigenvector $\mathbf{z}=(z_1, \cdots, z_n)^T$. In this case, the equilibrium state is given by

$\label{eq:N0s} \beta(N_0) = \dfrac{1}{\lambda^*}, N^*(X_i) = \frac{N_0 z_i}{\|\mathbf{z}\|_1}, \|\mathbf{z}\|_1=\sum\limits_{i=1}^n|z_i|.$

(16)

Next, we show that conditions (1)-(3) ensure the existence of the eigenvalue $\lambda^*$.

First, when $q_{ii}\geq\frac{1}{2}$, $A$ is a non-negative matrix. Moreover, since the transition matrix $P$ is irreducible, $A$ is irreducible as well¹. Thus, we have the following results according to Perron-Frobenius theorem:

¹The matrix $A$ can be rewritten as

$A=\Lambda_1P\Lambda_2-\Lambda_3,$

(17)

where $\Lambda_1, \Lambda_2, \Lambda_3$ are positive diagonal matrices

$\begin{eqnarray*} \Lambda_1&=&\mathrm{diag}(\dfrac{2(1-\kappa(X_1))}{\gamma(X_1)}, \cdots, \dfrac{2(1-\kappa(X_n))}{\gamma(X_n)}), \\ \Lambda_2&=&\mathrm{diag}(1-\mu(X_1), \cdots, 1-\mu(X_n)), \\ \Lambda_3&=&\mathrm{diag}(\dfrac{1}{\gamma(X_1)}, \cdots, \dfrac{1}{\gamma(X_n)}). \end{eqnarray*}$

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.

1. The spectral radius of $A$, denoted by $\lambda^*$, is a positive real number and an eigenvalue of $A$.

2. $\lambda^*$ is bounded by the sum of row elements of $A$, i.e.,

$s\leq \lambda^* \leq S.$

3. $\lambda^*$ is a simple eigenvalue and the corresponding unique eigenvector $\mathbf{z} = (z_1, \cdots, z_n)^T$ is positive (i.e. $z_i> 0$).

4. $\lambda^*$ is the only eigenvalue of $A$ with a non-negative eigenvector.

Next, since $\beta(N)$ is a continuous function, $\beta_{\min}\leq \beta(N)\leq \beta_{\max}\ (\forall N>0)$, and

$\dfrac{1}{\beta_{\max}} < s\leq \lambda^* \leq S < \dfrac{1}{\beta_{\min}},$

there exists at least one $N_0>0$ so that

$\dfrac{1}{\beta(N_0)} = \lambda^*.$

Finally, letting

$N^*(X_i) = \dfrac{N_0 z_i}{\|\mathbf{z}\|_1},$

it is easy to see that $N^*(X_i)\ (i=1, \cdots, 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}\geq \frac{1}{2}$ means that each cell has a probability not less than $\frac{1}{2}$ to generate a daughter cell with the same epigenetic state. In particular, if we omit the heterogeneity, then $q = (1-\mu)(1-\kappa)$. In this case, the condition $q \geq \frac{1}{2}$ is necessary to have a positive steady state.

Remark 5. From the proof, each eigenvalue of $A$ gives an equilibrium state if the eigenvalue satisfies $\lambda \in [\frac{1}{\beta_{\max}}, \frac{1}{\beta_{\min}}]$ and the corresponding eigenvector is non-negative. Furthermore, according to property 4 of Perron-Frobenius theorem, if the proliferation rate $\beta(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 $\mu$ 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) \equiv 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

$\label{mo_final1} \left\{ \begin{array}{rcl} \displaystyle\dfrac{\partial N(t, X_i)}{\partial t} &=& \displaystyle-N(t, X_i) (\beta(N) + \gamma(X_i))\\ &&\displaystyle {}+ 2 p(X_i) (1-\mu) (1 - \kappa(X_i)) \beta(N) N, \\ N&=&\displaystyle\sum\limits_{X_i\in\Omega} N(t, X_i). \end{array} \right.$

(18)

The following theorem gives a necessary and sufficient condition for the existence of a positive equilibrium state of the equation (18).

Theorem 3.2. Assume that $\beta(N)\in C(\mathbb{R}^+, \mathbb{R}^+)$, and

$\beta(N) > 0, \forall N\in \mathbb{R}^+, \mathrm{and}\ \gamma(X)>0, \ 0\leq \mu, \kappa(X), p(X) \leq 1, \forall X\in \Omega.$

Let

$\label{eq:bm} \beta_{\min} = \inf\limits_{N\in \mathbb{R}^+} \beta(N), \beta_{\max} = \sup\limits_{N\in \mathbb{R}^+} \beta(N),$

(19)

and

$\label{eq:hb} F(\beta) = 2 \sum\limits_{X_i\in \Omega} \dfrac{\beta}{\beta + \gamma(X_i)}p(X_i) (1 - \kappa(X_i)) ( 1 -\mu).$

(20)

Equation (18) has at least one positive equilibrium state if and only if

$F(\beta_{\min}) < 1 < F(\beta_{\max}),$

(21)

and the equilibrium state is given by

$N(t, X_i)\equiv N_0(X_i)=\dfrac{2\beta(N_0)N_0}{\beta(N_0)+\gamma(X_i)}p(X_i)(1-\kappa(X_i))(1-\mu),$

(22)

where $N_0$ satisfies $F(\beta(N_0)) = 1$. In particularly, if $\beta(N)$ is a monotone function, the positive equilibrium state is unique.

Proof. At the equilibrium state, equation (18) gives

$-N^*(X_i)(\beta(N_0)+\gamma(X_i))+2(1-\mu)p(X_i)(1-\kappa(X_i))N_0\beta(N_0)=0$

(23)

for all $X_i \in \Omega$, where

$\label{eq:n0s} N_0 = \sum\limits_{X_i\in \Omega} N^*(X_i).$

(24)

Hence, we obtain

$N^*(X_i) =\frac{2\beta(N_0) N_0}{\beta(N_0)+\gamma(X_i)} p(X_i) (1-\kappa(X_i)) (1-\mu). \label{e1}$

(25)

Substituting (25) into (24), we have

$\label{eq:hb2} F(\beta(N_0)) = 1.$

(26)

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 $\beta(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

$f(X_i) = \dfrac{2(1-\mu)p(X_i)(1-\kappa(X_i))\beta(N_0)}{\beta(N_0)+\gamma(X_i)}.$

(27)

---

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 $\beta(N)\in C^1(\mathbb{R}^+, \mathbb{R}^+)$. We first consider the system with two states, and the general situation of $n \geq 2$ epigenetic states is discussed in the next section.

Theorem 3.3. Consider equation (18) with two states, i.e., $\Omega=\{X_1, X_2\}$, and there is a positive equilibrium state with total cell population $N_0$, then the equilibrium state is asymptotically stable if and only if $\beta'(N_0)<0$.

Proof. When $\Omega = \{X_1, X_2\}$, equation (18) becomes

$\left\{ \begin{array}{rcl} \dfrac{\partial N(t, {X_1})}{\partial t}&=& - N(t, {X_1})(\beta (N(t)) + \gamma ({X_1}))\\ &&{} + 2(1 - \mu )p({X_1})(1 - \kappa ({X_1}))N(t)\beta (N(t))\\ &:=&f_1(N(t, X_1), N(t, X_2))\\ \dfrac{\partial N(t, {X_2})}{\partial t}&=& - N(t, {X_2})(\beta (N(t)) + \gamma ({X_2}))\\ &&{} + 2(1 - \mu )p({X_2})(1 - \kappa ({X_2}))N(t)\beta (N(t))\\ &:=&f_2(N(t, X_1), N(t, X_2))\\ N(t)&=&N(t, {X_1}) + N(t, {X_2}). \end{array} \right. \label{mo_2states}$

(28)

The linearization near the equilibrium state $\mathbf{N}^* = (N^*(X_1), N^*(X_2))^T$ is given by

$\dfrac{d \mathbf{x}}{d t} = A \mathbf{x},$

(29)

where $A = (a_{ij})$ is the Jacobi matrix

$A = \left.\left[ \begin{array}{*{20}{c}} \dfrac{\partial f_1}{\partial N(t, {X_1})}&\dfrac{\partial f_1}{\partial N(t, X_2)}\\ \dfrac{\partial f_2}{\partial N(t, {X_1})}&\dfrac{\partial f_2}{\partial N(t, X_2)} \end{array} \right]\right|_{\mathbf{N}=\mathbf{N}^*},$

where

$\frac{\partial f_1}{\partial N(t, X_2)} = -N(t, X_1)\beta'+2(1-\mu)p(X_1)(1-\kappa(X_1))(\beta+\beta'N),$

(30)

$\frac{\partial f_2}{\partial N(t, X_1)} = -N(t, X_2)\beta'+2(1-\mu)p(X_2)(1-\kappa(X_2))(\beta+\beta'N),$

(31)

$\frac{\partial f_1}{\partial N(t, X_1)} = \frac{\partial f_1}{\partial N(t, X_2)} -(\beta+\gamma(X_1)),$

(32)

$\frac{\partial f_2}{\partial N(t, X_2)} = \frac{\partial f_2}{\partial N(t, X_1)} -(\beta+\gamma(X_2)).$

(33)

Here we write $\beta, \beta', N$ for $\beta(N(t)), \beta'(N(t)), N(t)$, respectively, for short.

Thus, we obtain

$\begin{eqnarray*} &&\mathrm{Det}(A)\nonumber\\ &=&\left.\left(\frac{\partial f_1}{\partial N(t, X_1)}\frac{\partial f_2}{\partial N(t, X_2)}-\frac{\partial f_1}{\partial N(t, X_2)}\frac{\partial f_2}{\partial N(t, X_1)}\right)\right|_{\mathbf{N}=\mathbf{N}^*}\nonumber\\ &=&(\beta(N_0)+\gamma(X_1))(\beta(N_0)+\gamma(X_2))\nonumber\\ &&{}-(\beta(N_0)+\gamma(X_1))\left.\frac{\partial f_2}{\partial N(t, X_1)}\right|_{\mathbf{N}=\mathbf{N}^*}-(\beta(N_0)+\gamma(X_2))\left.\frac{\partial f_1}{\partial N(t, X_2)}\right|_{\mathbf{N}=\mathbf{N}^*}\nonumber\\ &=&(\beta(N_0)+\gamma(X_1))(\beta(N_0)+\gamma(X_2))\nonumber\\ &&\times\left.\left(1-\frac{1}{\beta(N_0)+\gamma(X_2)}\frac{\partial f_2}{\partial N(t, X_1)}-\frac{1}{\beta(N_0)+\gamma(X_1)}\frac{\partial f_1}{\partial N(t, X_2)}\right)\right|_{\mathbf{N}=\mathbf{N}^*}. \end{eqnarray*}$

(34)

Next, we analyze the last term in (34), which determines the sign of $\mathrm{Det}(A)$. From Theorem 3.2, and (24)-(26), we have

$\begin{eqnarray*} &&\left.\left(1-\frac{1}{\beta(N_0)+\gamma(X_2)}\frac{\partial f_2}{\partial N(t, X_1)}-\frac{1}{\beta(N_0)+\gamma(X_1)}\frac{\partial f_1}{\partial N(t, X_2)}\right)\right|_{\mathbf{N}=\mathbf{N}^*}\\ &=&1 - \dfrac{-N^*(X_2)\beta'(N_0)+2(1-\mu)p(X_2)(1-\kappa(X_2))(\beta(N_0)+\beta'(N_0) N_0)}{\beta(N_0) + \gamma(X_2)}\\ &&{} - \dfrac{-N^*(X_1)\beta'(N_0)+2(1-\mu)p(X_1)(1-\kappa(X_1))(\beta(N_0)+\beta'(N_0) N_0)}{\beta(N_0) + \gamma(X_1)}\\ &=&\frac{N^*(X_2)\beta'(N_0)-2(1-\mu)p(X_2)(1-\kappa(X_2))N_0\beta'(N_0)}{\beta(N_0)+\gamma(X_2)}\\ &&+\frac{N^*(X_1)\beta'(N_0)-2(1-\mu)p(X_1)(1-\kappa(X_1))N_0\beta'(N_0)}{\beta(N_0)+\gamma(X_1)}\\ &=&-\frac{\beta'(N_0)}{\beta(N_0)}\left(\frac{\gamma(X_2)N^*(X_2)}{\beta(N_0)+\gamma(X_2)}+\frac{\gamma(X_1)N^*(X_1)}{\beta(N_0)+\gamma(X_1)}\right). \end{eqnarray*}$

Hence, it is easy to see that $\mathrm{Det}(A)>0$ if and only if $\beta'(N_0)<0$.

Next, we have

$\begin{eqnarray*} &&\mathrm{tr}(A)\\ &=& -(2\beta(N_0)+\gamma(X_1)+\gamma(X_2))-N_0\beta'(N_0)\\ &&{}+\beta(N_0)(2(1-\mu)p(X_1)(1-\kappa(X_1))+2(1-\mu)p(X_2)(1-\kappa(X_2)))\\ &&{}+N_0\beta'(N_0)(2(1-\mu)p(X_1)(1-\kappa(X_1))+2(1-\mu)p(X_2)(1-\kappa(X_2)))\\ &=& -(\beta(N_0)+\gamma(X_1)+\gamma(X_2))\\ &&{}+\beta(N_0)\left(\frac{N^*(X_1)(\beta(N_0)+\gamma(X_1))}{N_0\beta(N_0)}+\frac{N^*(X_2)(\beta(N_0)+\gamma(X_2))}{N_0\beta(N_0)}-1\right)\\ &&{}+N_0\beta'(N_0)\left(\frac{N^*(X_1)(\beta(N_0)+\gamma(X_1))}{N_0\beta(N_0)}+\frac{N^*(X_2)(\beta(N_0)+\gamma(X_2))}{N_0\beta(N_0)}-1\right)\\ &=&-(\beta(N_0)+\gamma(X_1)+\gamma(X_2))+\left(\frac{N^*(X_1)}{N_0}\gamma(X_1)+\frac{N^*(X_2)}{N_0}\gamma(X_2)\right)\\ &&{}+\frac{\beta'(N_0)}{\beta(N_0)}(N^*(X_1)\gamma(X_1)+N^*(X_2)\gamma(X_2))\\ &=&-\beta(N_0)-(\gamma(X_1)+\gamma(X_2)-\frac{N^*(X_1)}{N_0}\gamma(X_1)-\frac{N^*(X_2)}{N_0}\gamma(X_2))\\ &&{}+\frac{\beta'(N_0)}{\beta(N_0)}(N^*(X_1)\gamma(X_1)+N^*(X_2)\gamma(X_2)). \end{eqnarray*}$

It is easy to verify that $\mathrm{tr}(A) < 0$ if $\beta'(N_0)<0$.

Now, if $\beta'(N_0) < 0$, we have $\mathrm{Det}(A) > 0$ and $\mathrm{tr}(A) < 0$, hence all eigenvalues of $A$ have negative real parts, which implies the stability of the positive equilibrium state. This theorem is proved.

---

3.3. Linear stability of the equilibrium state in multi-state systems

Now, we consider the linear stability of general cases with multiple epigenetic states under assumptions (A1)-(A3).

Theorem 3.4. Consider equation (18) with $n$ epigenetic states, i.e., $\Omega=\{X_1, \cdots, X_n\}$, and there is a positive equilibrium state with total cell population $N_0$.

(ⅰ) (The condition $\beta'(N_0) \leq0$ is necessary for the linear stability of the equilibrium state.

(ⅱ) If $\beta'(N_0)<0$, and the $n$ terms

$({\beta'}({N_0}){N_0}\dfrac{{\gamma ({X_i})}}{{\beta ({N_0}) + \gamma ({X_i})}} + \beta ({N_0})), (n=1, \cdots, n)$

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), \cdots, N(t, X_n))$ as

$f_i=-N(t, X_i)(\beta(N(t))+\gamma(X_i))+2(1-\mu)p(X_i)(1-\kappa(X_i))N(t)\beta(N(t)).$

Equations in (18) are rewritten as

$\dfrac{\partial N(t, X_i)}{\partial t} = f_i, (i=1, 2, \cdots, n)$

(35)

$N=\sum\limits_{i = 1}^n {N(t, {X_i})}.$

(36)

Let $\mathbf{N}^* = (N^*(X_1), \cdots, N^*(X_n))$ and $N_0 = \|\mathbf{N}^*\|_1$ be the equilibrium state, i.e.,

$\left. f_i \right|_{\mathbf{N}=\mathbf{N}^*} = 0.$

(37)

Similar to the calculations in the proof of Theorem 3.3, we obtain

$\left. {\frac{\partial {f_i}}{\partial N(t, {X_j})}} \right|_{\mathbf{N}=\mathbf{N}^*} = \left\{ \begin{array}{l} b_i - {a_i} & i = j\\ {b_i} & i \ne j \end{array} \right.$

(38)

where

$a_i=(\beta(N_0)+\gamma(X_i)),$

(39)

$b_i=2(1-\mu)p(X_i)(1-\kappa(X_i))(\beta(N_0)+\beta'(N_0)N_0) - N^*(X_i)\beta'(N_0).$

(40)

Hence, the coefficient matrix $A$ of the linearization of (35) at the equilibrium state $\mathbf{N} = \mathbf{N}^*$ is given by

$\begin{eqnarray*} A &=& \left.\left[ \begin{array}{*{20}{c}} \dfrac{\partial {f_1}}{\partial N(t, {X_1})}&\dfrac{\partial {f_1}}{\partial N(t, {X_2})}& \cdots &\dfrac{\partial {f_1}}{\partial N(t, {X_n})}\\ \dfrac{\partial {f_2}}{\partial N(t, {X_1})}&\dfrac{\partial {f_2}}{\partial N(t, {X_2})}& \cdots &\dfrac{\partial {f_2}}{\partial N(t, {X_n})}\\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial {f_n}}{\partial N(t, {X_1})}&\dfrac{\partial {f_n}}{\partial N(t, {X_2})}& \cdots &\dfrac{\partial {f_n}}{\partial N(t, {X_n})} \end{array} \right]\right|_{\mathbf{N}=\mathbf{N}^*}\\ &=&\left[ \begin{array}{*{20}{c}} b_1 - a_1 &b_1& \cdots &b_1\\ b_2& b_2 - a_2 & \cdots &b_2\\ \vdots & \vdots & \ddots & \vdots \\ b_n&b_n& \cdots & b_n - a_n \end{array} \right]. \end{eqnarray*}$

Therefore, the eigenfunction of the coefficient matrix is given by

$\begin{eqnarray*} h(\lambda) &=& \left| {\lambda I - A} \right| \nonumber\\ &=& \left| {\begin{array}{*{20}{c}} {\lambda + {a_1} - {b_1}}&{ - {b_1}}& \cdots &{ - {b_1}}\\ { - {b_2}}&{\lambda + {a_2} - {b_2}}& \cdots &{ - {b_2}}\\ \vdots & \vdots & \ddots & \vdots \\ { - {b_n}}&{ - {b_n}}& \cdots &{\lambda + {a_n} - {b_n}} \end{array}} \right|\nonumber\\ &=& (1 - \frac{b_1}{\lambda + a_1} - \cdots - \frac{b_n}{\lambda + a_n})(\lambda + a_1) \cdots (\lambda + a_n).\label{th_3_ja} \end{eqnarray*}$

(41)

From (41), we have the following results for $h(\lambda)$:

$h(+\infty) = +\infty>0,$

(42)

$h(0) = (1-\frac{b_1}{a_1}-\frac{b_2}{a_2}-\cdots-\frac{b_n}{a_n})\prod\limits_{i=1}^n a_i,$

(43)

$h(-a_i) = -b_i\prod\limits_{j \ne i} {({a_j} - {a_i})},$

(44)

$4h(-\infty) = {\mathop{\rm sgn}} ({( - 1)^{n}})\cdot\infty.$

(45)

From (25) and (40),

$\begin{eqnarray*} b_i&=&-N^*(X_i)\beta'(N_0)+2(1-\mu)p(X_i)(1-\kappa(X_i))(\beta(N_0)+\beta'(N_0)N_0)\nonumber\\ &=&2(1-\mu)p(X_i)(1-\kappa(X_i))\beta'(N_0)N_0(-\frac{\beta(N_0)}{\beta(N_0)+\gamma(X_i)}+1)\nonumber\\ &&{}+2(1-\mu)p(X_i)(1-\kappa(X_i))\beta(N_0)\nonumber\\ \label{eq:bi} &=&2(1-\mu)p(X_i)(1-\kappa(X_i))(\beta'(N_0)N_0\frac{\gamma(X_i)}{\beta(N_0)+\gamma(X_i)}+\beta(N_0)). \end{eqnarray*}$

(46)

First, we assume that $b_i\not= 0$. Then (46) and (39) together yield the following expression

$\dfrac{b_i}{a_i}=\dfrac{2(1-\mu)p(X_i)(1-\kappa(X_i))\beta(N_0)}{\beta(N_0)+\gamma(X_i)} + N_0 \beta'(N_0) \dfrac{2 \gamma(X_i) (1-\mu) p(X_i) (1-\kappa(X_i))}{(\beta(N_0)+\gamma(X_i))^2}. \label{th_3_biai}$

(47)

Thus, from (26), we obtain

$\begin{eqnarray*} &&1 - \sum\limits_{i = 1}^n \dfrac{b_i}{a_i}\nonumber\\ &=&1 - \sum\limits_{i = 1}^n \frac{2(1-\mu)p(X_i)(1-\kappa(X_i))\beta(N_0)}{\beta(N_0)+\gamma(X_i)}\nonumber\\ &&{} - N_0 \beta'(N_0) \sum\limits_{i = 1}^n \dfrac{2 \gamma(X_i) (1-\mu) p(X_i) (1-\kappa(X_i))}{(\beta(N_0)+\gamma(X_i))^2}\nonumber\\ \label{eq:1ab} &=&- N_0 \beta'(N_0) \sum\limits_{i = 1}^n \dfrac{2 \gamma(X_i) (1-\mu) p(X_i) (1-\kappa(X_i))}{(\beta(N_0)+\gamma(X_i))^2}. \end{eqnarray*}$

(48)

Now, we are ready to prove the main results.

(ⅰ). If $\beta'(N_0)>0$, since $a_i > 0 (i=1, \cdots, n)$, (48) implies $h(0) < 0$. However, $h(\lambda) > 0$ for large enough $\lambda$, hence $h(\lambda) = 0$ has a positive solution, and the equilibrium state is unstable. Thus, $\beta'(N_0)\leq 0$ is necessary for the linear stability of $\mathbf{N}^*$.

(ⅱ). First, we assume that $a_i$ are not equal to each other, without loss of generality, we assume $0 < a_1 < a_2 < \cdots < a_n$.

    (a) If $b_i\ (i=1, \cdots, n)$ are positive, from (44), we have

$h(0)>0, h(-a_1)<0, h(-a_2)>0, \cdots, \mathrm{sgn}(h(-a_n))= \mathrm{sgn}{( - 1)^{n}},$

    and there are $n$ negative roots between $(a_n, 0)$.

    (b) If otherwise, $b_i\ (i=1, \cdots, n)$ are negative, from (45), we have

$h(-a_1)>0, h(-a_2)<0, \cdots, \mathrm{sgn}(h(-a_n))= \mathrm{sgn}(( - 1)^{n + 1}),$

$\mathrm{sgn}(h(-\infty))=\mathrm{sgn}(( - 1)^{n}),$

    and there are $n$ negative roots between $(-\infty, a_1)$.

Next, if $b_i = 0$, it is easy to verify that $\lambda=-a_i < 0$ is a root of $h(\lambda)$. And if some $a_i = a_{i+1}=\cdots=a_{i+k}$, $\lambda = -a_i<0$ is a root of $h(\lambda)$ of multiplicity $k$. Under either situations, we can write $h(\lambda)$ as $(\lambda+a_i)^k \tilde{h}(\lambda)$ so that $\tilde{h}(-a_i)\not= 0$. Then, all roots of $\tilde{h}(\lambda)$ are negative following the similar discussions above.

Thus, in all situations, all $n$ roots of the eigenfunction are negative, and the equilibrium is linearly stable.

Remark 8. Since $\gamma(X_i)\leq1$ and $\beta'(N_0)<0$, it is easy to verify that

${\beta'}({N_0}){N_0}\dfrac{{1}}{{\beta ({N_0}) + 1}} + \beta ({N_0})>0\label{th_3_re}$

(49)

is enough to ensure the linear stability of the positive equilibrium. Moreover, if $\beta(N)$ is taken as a Hill type function (the function widely used in many stem cell population models [17,24,27])

$\beta(N)=\beta_0\frac{\theta^n}{\theta^n+N^n}, \label{hill}$

(50)

the condition (49) is equivalent to

$(\beta_0+1) > (n-1) (N_0/\theta)^n. \label{th_3_re2}$

(51)

Interestingly, we note that (51) is always satisfied when $n=1$.

---

3.4. Stability of the zero solution

Here, we study the stability of the zero solution $\mathbf{0}$ under assumptions (A1)-(A3), and prove that if equation (18) has a positive equilibrium state, the zero solution is unstable.

Theorem 3.5. Consider the equation (18), and $\beta(N)$ is a decreasing function. Let

$F(\beta) = \sum\limits_{i=1}^n \dfrac{2 (1-\mu) p(X_i) (1-\kappa(X_i)) \beta}{\beta + \gamma(X_i)}.$

(52)

(1) If there exists $N_0 > 0$ so that $F(\beta(N_0)) = 1$, then the zero solution $\mathbf{N}(t)\equiv \mathbf{0}$ is unstable.

(2) If $F(\beta(0)) < 1$, then the zero solution is linearly stable.

Proof. Consider the linearization of (18) near the zero solution, the coefficient matrix is given by

$\bar{A} = \left[ \begin{array}{cccc} \bar{b}_1 - \bar{a}_1 & \bar{b}_1 & \cdots & \bar{b}_1\\ \bar{b}_2 & \bar{b}_2 - \bar{a}_2 & \cdots & \bar{b}_2\\ \vdots & \vdots & \ddots & \vdots\\ \bar{b}_n & \bar{b}_n & \cdots & \bar{b}_n - \bar{a}_n \end{array} \right],$

where

$\begin{eqnarray*} \bar{a}_i &=&(\beta(0)+\gamma(X_i))>0\\ \bar{b}_i &=&2(1-\mu)p(X_i)(1-\kappa(X_i))\beta(0)>0. \end{eqnarray*}$

Thus, similar to the above argument, let $\bar{h}(\lambda)$ the eigenfunction of $\bar{A}$, we have

$\begin{eqnarray*} \bar{h}(+\infty)&=&+\infty>0, \\ \bar{h}(0)&=&(1- \frac{\bar{b}_1}{\bar{a}_1} - \frac{\bar{b}_2}{\bar{a}_2} - \cdots - \frac{\bar{b}_n}{\bar{a}_n})\prod\limits_{i=1}^n \bar{a}_i, \\ \bar{h}(-\bar{a}_i) &=& -\bar{b}_i \prod\limits_{j\not=i} (\bar{a}_j - \bar{a}_i). \end{eqnarray*}$

(1) If there exists $N_0>0$ so that $F(\beta(N_0)) = 1$, then $\beta(N_0) < \beta(0)$, and

$\begin{eqnarray*} &&1-\sum\limits_{i=1}^n \frac{\bar{b}_i}{\bar{a}_i}\\ & = &1 - \sum\limits_{i = 1}^n {\frac{2(1-\mu)p(X_i)(1-\kappa(X_i))\beta(0)}{\beta(0)+\gamma(X_i)}}\\ &=& F(\beta(N_0)) - F(\beta(0)). \end{eqnarray*}$

From the proof of Theorem 3.2, $F(\beta)$ is a monotone increasing function of $\beta$. Hence, $1-\sum_{i=1}^n \frac{\bar{b}_i}{\bar{a}_i}<0$, which implies $\bar{h}(0) < 0$ by $\bar{a}_i > 0\ (i=1\cdots, n)$. Thus, there exists an eigenvalue of $\bar{h}(\lambda)$ in $(0, +\infty)$, and the zero solution is unstable.

(2) If $F(\beta(0)) < 1$, then $\bar{h}(0) > 0$. Moreover, since $b_i > 0\ (i=1, \cdots, n)$, similar to the discussions in Theorem 3.4, assuming $0<\bar{a}_1<\bar{a}_2<\cdots<\bar{a}_n$, we have

$\bar{h}(0)>0, \bar{h}(-\bar{a}_1)<0, \bar{h}(-\bar{a}_2)>0, \cdots$

and all roots of $\bar{h}(\lambda)$ have negative real parts, hence the zero solution is linearly stable.

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(\beta(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(\beta(0)) = 1$).

---

4. 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 $\Omega = \{X_i\}_{i=1}^{300}$. Here we write $X_i = i$ for simplicity. In simulations, the functions $\mu(X_i)$, $\kappa(X_i)$ and $\gamma(X_i)$ are taken as

$\mu(X_i ) = 0.1, \kappa(X_i) = 0.1 + \dfrac{0.2}{1 + (X_i/30)}, \kappa(X_i) = 0.05 + \dfrac{0.35}{1 + (X_i/50)}.$

(53)

The proliferation rate function is given by

$\label{eq:beta} \beta(N) = \beta_0 \dfrac{\theta^n}{\theta^n + N^n}, \beta_0 = 0.7, \theta = 500, n= 3.$

(54)

The transition probability $P(X_i, X_j)$ is defined as

$p(X_i, X_j) = \dfrac{(300 - X_i) (1 + \nu e^{-0.2 (X_i - X_j)^2})}{\sum_{i=1}^{300} (300 - X_i) (1 + \nu e^{-0.2 (X_i - X_j)^2})}.$

(55)

Here the parameter $\nu$ is introduced to represent the effect of mother cells. When $\nu = 0$, the transition probability is independent to the state of the mother cell. When $\nu$ 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 $\nu = 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 $\nu > 0$, the transition matrix is different from the matrix for the case $\nu=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 $\nu$ 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 $\nu$ further increases to $200$, there is an obvious change in the equilibrium state.

---

4.2. State transitions of cancer cells

Experiments have found that in human breast cancer cell lines there are distinct phenotypic cells (stem-like (S), basal (B), or luminal (L)) [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.

The model includes three subpopulations of cells ($\Omega = \{S, B, L\}$), and cells can transit between different states. The parameters are listed at 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.

---

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.

---