A general mathematical framework for understanding the behavior of heterogeneous stem cell regeneration

TL;DR

This paper introduces a comprehensive mathematical framework that models stem cell regeneration dynamics considering heterogeneity and epigenetic state transitions, extending classical models to include stochastic inheritance and potential tumor development insights.

Contribution

It develops a generalized model incorporating continuous epigenetic states and stochastic inheritance, advancing understanding of stem cell behavior and tissue dynamics beyond previous models.

Findings

The model captures tissue development, degeneration, and abnormal growth processes.

It provides a basis for studying gene mutation effects on tumor development.

The framework generalizes classical G0 cell cycle models.

Abstract

Stem cell heterogeneity is essential for the homeostasis in tissue development. This paper established a general formulation for understanding the dynamics of stem cell regeneration with cell heterogeneity and random transitions of epigenetic states. The model generalizes the classical G0 cell cycle model, and incorporates the epigenetic states of stem cells that are represented by a continuous multidimensional variable and the kinetic rates of cell behaviors, including proliferation, differentiation, and apoptosis, that are dependent on their epigenetic states. Moreover, the random transition of epigenetic states is represented by an inheritance probability that can be described as a conditional beta distribution. This model can be extended to investigate gene mutation-induced tumor development. The proposed formula is a generalized formula that helps us to understand various dynamic processes of stem cell regeneration, including tissue development, degeneration, and abnormal growth.

Equations

$$\frac{dQ}{dt}=-(\beta (Q)+\kappa )Q+2e^{-\mu \tau }\beta (Q_{\tau })Q_{\tau }.$$

$$\left\{\begin{array}[]{rcl}\dfrac{\partial Q(t,\mathbf{x})}{\partial t}&=&-Q(t,\mathbf{x})(\beta(\hat{Q},\mathbf{x})+\kappa(\mathbf{x}))\\ &&{}\displaystyle+2\int\beta(\hat{Q}_{\tau(\mathbf{y})},\mathbf{y})Q(t-\tau(\mathbf{y}),\mathbf{y})e^{-\mu(\mathbf{y})\tau(\mathbf{y})}p(\mathbf{x},\mathbf{y})d\mathbf{y}\\ \hat{Q}(t)&=&\displaystyle\int Q(t,\mathbf{x})\xi(\mathbf{x})d\mathbf{x}.\end{array}\right.$$

$$\frac{dQ}{dt}$$

$$\frac{\partial f(t,\mathbf{x})}{\partial t}$$

$$p_i(x_i,\mathbf{y})=\frac{x_i^{a_i(\mathbf{y})}(1-x_i{)}^{b_i(\mathbf{y})}}{B(a_i(\mathbf{y}),b_i(\mathbf{y}))},B(a,b)=\frac{\Gamma (a)\Gamma (b)}{\Gamma (a+b)},$$

$$\mathrm{E}(x_i|\mathbf{y})={\varphi }_i(\mathbf{y}),\mathrm{Var}(x_i|\mathbf{y})=\frac{1}{1+{\eta }_i(\mathbf{y})}{\varphi }_i(\mathbf{y})(1-{\varphi }_i(\mathbf{y})),$$

$$a_i(\mathbf{y})={\eta }_i(\mathbf{y}){\varphi }_i(\mathbf{y}),b_i(\mathbf{y})={\eta }_i(\mathbf{y})(1-{\varphi }_i(\mathbf{y})).$$

$$0<{\varphi }_i(\mathbf{y})<1,{\eta }_i(\mathbf{y})>0.$$

$$\frac{\partial Q_i(t,\mathbf{x})}{\partial t}$$

$$\hat{Q}(t)=\sum _{i=0}^3\int Q_i(t,\mathbf{x})\xi (\mathbf{x})d\mathbf{x}.$$

$$\left\{\begin{array}[]{rcl}\dfrac{\partial Q(t,\mathbf{x})}{\partial t}&=&-Q(t,\mathbf{x})(\beta(\hat{Q},\mathbf{x})+\kappa(\mathbf{x}))\\ &&{}\displaystyle+2\int\beta(\hat{Q}_{\tau(\mathbf{y})},\mathbf{y})Q(t-\tau(\mathbf{y}),\mathbf{y})e^{-\mu(\mathbf{y})\tau(\mathbf{y})}p(\mathbf{x},\mathbf{y})d\mathbf{y}\\ \dfrac{\partial P(t,\mathbf{x})}{\partial t}&=&\kappa(\mathbf{x})Q(t,\mathbf{x})-\nu(\mathbf{x})P(t,\mathbf{x}).\end{array}\right.$$

$$Q_{\mathrm{stem}}(t)={\int }_{x_0}^{1}Q(t,x)dx,Q_{\mathrm{progenitor}}(t)={\int }_0^{x_0}Q(t,x)dx,P(t)={\int }_0^{1}P(t,x)dx.$$

$$S(t)=Q_{\mathrm{stem}}(t)+Q_{\mathrm{progenitor}}(t)+P(t),$$

$$f(t,x)=\frac{Q(t,x)+P(t,x)}{S(t)}.$$

$$\begin{array}[]{lcl}\nabla s(t,a)=-\mu s(t,a),&(t>0,0<a<\tau)\\ \dfrac{dQ}{dt}=2s(t,\tau)-(\beta(Q)+\kappa)Q,&(t>0)\end{array}$$

$$s(t,0)=\beta (Q(t))Q(t),$$

$$s(0,a)=g(a),(0\le a\le \tau );Q(0)=Q_0.$$

$$\dfrac{dQ}{dt}=-(\beta(Q)+\kappa)Q+\left\{\begin{array}[]{lcl}2g(\tau-t)e^{-\mu t},&0<t\leq\tau\\ 2e^{-\mu\tau}\beta(Q_{\tau})Q_{\tau},&t>\tau\end{array}\right.$$

$$\frac{dQ}{dt}=-(\beta (Q)+\kappa )Q+2e^{-\mu \tau }\beta (Q_{\tau })Q_{\tau }.$$

$$[\mathrm{R}]+[{\mathrm{R}}^{\ast }]=mQ,$$

$$[\mathrm{R}]+n[\mathrm{L}]\leftrightarrows [{\mathrm{R}}^{\ast }].$$

$$[\mathrm{R}][\mathrm{L}{]}^n=K[{\mathrm{R}}^{\ast }]$$

$$\beta ={\beta }_0\frac{[\mathrm{R}]}{mQ}.$$

$$\frac{[\mathrm{R}]}{mQ}=\frac{1}{1+[\mathrm{L}{]}^n}.$$

$$\beta (Q)={\beta }_0\frac{{\theta }^n}{{\theta }^n+Q^n},$$

$$\beta (Q)={\beta }_0\frac{{\theta }^n}{{\theta }^n+Q^n}+{\beta }_1.$$

$$-(\beta (Q^{\ast })+\kappa )Q^{\ast }+2e^{-\mu \tau }\beta (Q^{\ast })Q^{\ast }=0,$$

$$\beta (Q^{\ast })=\frac{\kappa }{2e^{-\mu \tau }-1}.$$

$${\beta }_0>\frac{\kappa }{2e^{-\mu \tau }-1}-{\beta }_1>0.$$

$${\beta }_1\ge \frac{\kappa }{2e^{-\mu \tau }-1},$$

Full text

A general mathematical framework for understanding the behavior of heterogeneous stem cell regeneration

Jinzhi Lei

Zhou Pei-Yuan Center for Applied Mathematics, MOE Key Laboratory of Bioinformatics, Tsinghua University, Beijing 100084, China.

Abstract  Stem cell heterogeneity is essential for the homeostasis in tissue development. This paper established a general formulation for understanding the dynamics of stem cell regeneration with cell heterogeneity and random transitions of epigenetic states. The model generalizes the classical G0 cell cycle model, and incorporates the epigenetic states of stem cells that are represented by a continuous multidimensional variable and the kinetic rates of cell behaviors, including proliferation, differentiation, and apoptosis, that are dependent on their epigenetic states. Moreover, the random transition of epigenetic states is represented by an inheritance probability that can be described as a conditional beta distribution. This model can be extended to investigate gene mutation-induced tumor development. The proposed formula is a generalized formula that helps us to understand various dynamic processes of stem cell regeneration, including tissue development, degeneration, and abnormal growth.

Keywords: Heterogeneity, stem cell, cell cycle, epigenetic state, development, computational model

1 Introduction

Stem cell regeneration is an essential biological process in most self-renewing tissues during development and the maintenance of tissue homeostasis. Stem cells multiply by cell division, during which DNA is replicated and assigned to the two daughter cells along with the inheritance of epigenetic information and the partition of molecules. Unlike the accumulated process of DNA replication, inherited epigenetic information is often subjected to random perturbations; for example, the reconstruction of histone modifications and DNA methylation are intrinsically random processes of writing and erasing the modified markers [71, 91]. The stochastic inheritance of epigenetic changes during cell division can lead to stem cell heterogeneity which is important for the dynamic equilibrium of various phenotypic cells during tissue development. Accumulation of undesirable epigenetic changes may result in promoting or causing diseases [17, 18, 36, 43, 59, 60, 70, 72, 81, 92].

The heterogeneity of stem cells has been highlighted in recent years due to new technologies with single-cell resolution, which have led to the discovery of new cell types and changes in the understanding of differentiation landscapes [6, 9, 29, 50, 51, 69]. In early embryonic development, heterogeneous expression and histone modifications are correlated with correlated with cell fate and the dynamic equilibrium of pluripotent stem cells [33, 34, 68, 84]. Chromatin modifications in the human primary hematopoietic stem cell/progenitor cell (HSC/HPC) stage can lead to the dynamic equilibrium of heterogeneous and interconvertible HSCs [10, 89], as well as gene expression changes during differentiations [11]. Moreover, applications of single-cell RNA sequencing have revealed the continuous spectrum of differentiation in zebrafish [54], mice [65], and human HSCs [85]. These findings have challenged the demarcation between stem cells and progenitor cells and have led to the evolving understanding of the complex hematopoietic differentiation landscape [47, 66].

Heterogeneity plays an important role in the development of drug resistance. Cancer development is driven by evolutionary selection on somatic genetic alterations and epigenetic alterations, which result in the multistage tumorigenesis and heterogenous cancer cell phenotypes [22, 23, 39, 51, 53, 64, 90]. Tumors with different subtypes often differ in the treatment response and patient survival [13, 23, 73], and treatment stress can also induce cancer cell plasticity and drug resistance [24, 48, 57, 80, 82]. Cell plasticity is often associated with epigenetic modifications, and targeting the epigenetic regulators, such as the polycomb group protein EZH2, has been an attractive strategy in cancer treatment [20, 77, 83]. To better understand the progress of tumorigenesis and drug resistance, we need to develop predictive models of the evolutionary dynamics of cancer [3, 31].

Despite the central role of stem cell regeneration in tissue development, a quantitative investigation of the process is well beyond the ability of current technologies. Furthermore, in many fields of biological science, mathematical modeling tools have aided in improving the understanding of the principles of related processes [2, 45, 61, 62]. In 2007, Weinberg posed the following question [86]: “can algebraic formulae tell us more than reasoning about the behavior of complex biological systems?” Various computational models have been established in studies of tissue development and cancer systems biology under different circumstances [3, 4, 16, 25, 26, 27, 87, 88]. Nevertheless, a unified formulation that bypasses detailed assumptions is required to provide more basic logic of the biological behaviors of these complex systems. In this study, based on the general process of the cell cycle and heterogeneous stem cell regeneration, we established a general mathematical framework to formulate the dynamics of heterogeneous stem cell regeneration. The model framework includes essential cellular behaviors, including proliferation, apoptosis and differentiation/senescence; however, it bypasses the biological details of signaling pathways. The heterogeneity of stem cells and epigenetic inheritance during the cell cycle are key points in model development. Various formulas can be applied to different processes, such as embryonic development, tissue disease and degeneration, and tumor development.

The aim of this paper was to introduce a new general formula for the dynamics of stem cell regeneration with an emphasis on the effects of cell heterogeneity; therefore, a discussion of concrete conclusions based on the formula was not included. The simulation results below were included to demonstrate the potent application of the model and were not related to any actual biological processes.

2 Results

2.1 The G0 cell cycle model for homogeneous stem cell regeneration

A classical model that is used to describe the dynamics of stem cell regeneration is the G0 cell cycle model proposed in the 1970s [7, 55]. In this model, homogeneous stem cell cycles are classified into resting (G0) or proliferating (G1, S, and G2 phases and mitosis) phases (Figure 1A). During each cell cycle, a cell in the proliferating phase either undergoes apoptosis or divides into two daughter cells; however, a cell in the resting phase either irreversibly differentiates into a terminally differentiated cell or returns to the proliferating phase. This can be modeled by an age-structure model for cell numbers in the resting phase and proliferating phase. Integrating the age-structure model through the characteristic line method provides the following delay differential equation (Material and methods)

[TABLE]

Here, $\beta(Q)$ is the proliferation rate, $\mu$ is the apoptosis rate of cells in the proliferating phase, $\tau$ is the duration of the proliferating phase, and $\kappa$ is the rate of removing cells out of the resting phase, which includes terminally differentiation, cell death, and senescence (hereafter, we call $\kappa$ the differentiation rate for simplicity). Hereafter, the subscript indicates a time delay, i.e., $Q_{\tau}$ indicates $Q(t-\tau)$. The proliferation rate $\beta(Q)$ describes how cells regulate the self-renewal of stem cells through secreted cytokines and is often given by a decrease function and $\beta_{0}<\beta(Q)<\beta_{\infty}$ (Material and method). Typically, for normal individuals, we usually have $\beta_{\infty}=\lim_{Q\to+\infty}\beta(Q)=0$⁡ because of the inhibition of the cell cycle pathway.

The G0 cell cycle model and its extensions are widely used to investigate hematopoietic stem cell dynamics [1, 19, 49, 56]; dysregulation of the apoptosis rate or differentiation rate of hematopoietic stem cells can result in serious periodic hematopoietic diseases [12]. Moreover, from (1), the stem cell dynamics are mainly determined by pathways related to stem cell proliferation, apoptosis, differentiation, senescence, and growth. Major oncogenic signaling pathways obtained from an integrated analysis of genetic alterations in The Cancer Genome Atlas (TCGA) [74] show direction connections to the coefficients $\beta(Q),\mu$, and $\kappa$ in (1) (Figure 1B) (Material and methods). Equation (1) is capable of describing the population dynamics of stem cell regeneration. Nevertheless, cell heterogeneity is not included in the model and has been highlighted in recent years for the understanding of cancer development and drug resistance in cancer therapy.

2.2 The general framework of heterogeneous stem cell regeneration

To extend the abovementioned G0 cell cycle model to include cell heterogeneity, we introduce a quantity $\mathbf{x}$ (scalar or vector) for the epigenetic state of a cell and denote $Q(t,\mathbf{x})$ as the cell number at time $t$ with state $\mathbf{x}$ (Figure 2A). In general, $\mathbf{x}$ can refer to the expression levels of marker genes, histone modifications in nucleosomes, or DNA methylations associated with DNA segments and can be measured by single-cell sequencing techniques. Specifically, we often refer to $\mathbf{x}$ as quantities that affect signaling pathways that control cell cycle progression, apoptosis, and cell growth, so that the coefficients $\beta$, $\mu$, and $\kappa$ and the duration of the proliferating phase $\tau$ in (1) are cell specific and dependent on the state $\mathbf{x}$ in the cell. Moreover, cells in the niche can interfere with stem cell self-renewal through released cytokines. Let $\xi(\mathbf{x})$ denote the effective cytokine signal produced by a cell with state $\mathbf{x}$, and $\hat{Q}=\int Q(t,\mathbf{x})\xi(\mathbf{x})d\mathbf{x}$ denotes the total concentration of effective cytokines that regulate cell proliferation. The proliferation rate in (1) becomes $\beta(\hat{Q},\mathbf{x})$.

While the cell-to-cell variability is considered, the inheritance of epigenetic states of cells during cell division is essential to shape the distribution of cell heterogeneity. Many biological processes, such as the random partition of molecules [38], random inheritance of nucleosome modification [14, 71] and DNA methylation [91, 36], can be involved in the affecting the inheritance of epigenetic states from mother to daughter cells after cell division. Many efforts have been made to model epigenetic cell memory [14, 30, 37, 38, 79]; however, it remains challenging to develop precise models for this process that is not yet clear. Nevertheless, while we overlook the biological details and focus on the changes of epigenetic states, we introduce the inheritance probability $p(\mathbf{x},\mathbf{y})$, which represents the probability that a daughter cell of state $\mathbf{x}$ comes from a mother cell of state $\mathbf{y}$ after cell division. Therefore, $p(\mathbf{x},\mathbf{y})d\mathbf{x}=1$ for any $\mathbf{y}$. Based on the abovementioned assumptions and the similar argument to the abovementioned homogeneous model, the dynamical equation for $Q(t,\mathbf{x})$ is as follows (Material and methods):

[TABLE]

Here, the integrals are taken over all possible epigenetic states. Moreover, if we consider discrete states, such as gene mutations, we can extend the integrals to include the summation over all discrete states. Equation (2) extends the previous G0 cell cycle model and provides a general framework for heterogeneous stem cell regeneration.

Equation (2) is an autonomous system in which the rate functions and inheritance probability are not explicitly dependent on the time t. Nevertheless, while we apply the equation to situations with time dependence, such as embryo development, environmental changes, injury, and external stimuli, the time-dependent rate functions and the inheritance probability can be included in a straightforward manner. An example is shown in Figure 5 below.

Based on the (2), while we introduce appropriate definitions for the dependence of kinetic rates on epigenetic states and the inheritance function, we are able to model various processes of stem cell regeneration, such as tissue growth, degeneration, and abnormal growth (Figures 2B-G). For simplicity, in all simulations shown in this paper, we consider one epigenetic state $x$ ($0\leq x\leq 1$) that affects only cell proliferation and differentiation in a manner similar to the stemness so that a larger value of x indicates higher stemness (Material and methods). Figure 2B-C shows the dynamics of tissue growth starting from a small population cells with high levels of stemness toward a steady state. There is a temporary transition at the early stage characterized by a rapid increase in the cell number and a subpopulation of cells with low level stemness (Figures 2B-C, red arrows). Figures 2D-E shows the dynamics of degeneration with alterations to the inheritance function, and Figures 2F-G shows the abnormal growth due to a decreasing differentiation rate and an alteration to the inheritance function. Both processes include a short-term stage of biphenotypic cell populations with both high and low stemness cells (Figure 2E and G, red curve). Moreover, the simulations show that the steady state heterogeneity can be restored from cell subpopulation fractions (Figure 3), which is in agreement with experiments that were previously explained by transcriptome-wide noise [35, 52, 89].

Equation (2) describes the evolution of the cell numbers with various epigenetic states; however, the total cell number $Q(t)=\int Q(t,\mathbf{x})d\mathbf{x}$ and the density of cells with different epigenetic states $f(t,\mathbf{x})=Q(t,\mathbf{x})/Q(t)$ are relevant to the data. From equation (2), it is easy to obtain the equations for $Q(t)$ and $f(t,\mathbf{x})$ as follows (see Materials and methods):

[TABLE]

and

[TABLE]

Equations (2.2) and (2.2) provide the evolution dynamics of relative cell numbers that can be obtained from experiments by single-cell sequencing or flow cytometry. Here we note that when $\xi(\mathbf{x})=1$, we have $\hat{Q}(t)=Q(t)$, and hence, (2.2)-(2.2) provide a closed-form equation.

Here, the state variable $\mathbf{x}$ represents the epigenetic state, and $p(\mathbf{x},\mathbf{y})$ represents the inheritance function; hence, (2) describes the dynamics with epigenetic state transitions. Nevertheless, this equation can also describe the changes in genetic alternations if we consider $\mathbf{x}$ as the genetic state and $p(\mathbf{x},\mathbf{y})$ as the probability of point mutations. This is often the situation of genomic instability associated with cancer development [8, 32, 78, 93], and hence, the model can be used to study genetic heterogeneity in cancer development. In this paper, we focus on the equation with epigenetic state transitions and assumed that $\mathbf{x}$ always represents the epigenetic state.

2.3 Stochastic epigenetic state inheritance in the cell cycle

In Equations (2)-(2.2), the mathematical formula of epigenetic state-relevant coefficients should be expressed based on how the epigenetic states (or genes) affect the relative biological process. However, the inheritance function $p(\mathbf{x},\mathbf{y})$ cannot be determined from the biological process of cell division. Here, we derive a phenomenological inheritance function to represent the stochastic inheritance of the epigenetic states. More specifically, let $\mathbf{x}=(x_{1},x_{2},\cdots,x_{n})$ represent the expression level of n marker genes, and derive the inheritance function $p_{i}(x_{i},\mathbf{y})$ for each gene, and $p(\mathbf{x},\mathbf{y})=\prod_{i=1}^{n}p_{i}(x_{i},\mathbf{y})$.

We assume that the epigenetic state of a daughter cell is a random number with the distribution depending on that of the mother cell. In previous studies based on stochastic simulations of gene expression coupled with nucleosome modifications over multiple cell cycles [37, 40], we found that the nucleosome modification level of daughter cells, considering the nucleosome modification level of mother cells, which is normalized to the interval [0,1], can be well-described by a beta-distributed random number. Therefore, we generalized these findings and defined the inheritance function $p_{i}(x_{i},\mathbf{y})$ through the beta distribution density function as follows:

[TABLE]

where $\Gamma(z)$ is the gamma function and $a$ and $b$ are shape parameters depending on the state of the mother cell. We assumed that the mean and variance of $x_{i}$, considering $\mathbf{y}$, is as follows:

[TABLE]

and the shape parameters are (Materials and methods)

[TABLE]

Here, we note that $\phi_{i}(\mathbf{y})$ and $\eta_{i}(\mathbf{y})$ always satisfy

[TABLE]

Hence, the inheritance function can be determined through the predefined functions $\phi_{i}(\mathbf{y})$ and $\eta_{i}(\mathbf{y})$, often through data-driven modeling or assumptions, that satisfies (7).

2.4 Modeling tumor development with cell-to-cell variance

As shown in Figure 2F-G, to mimic the process of abnormal cell growth, we varied the differentiation rate and the inheritance probability. These variances to the model parameters can be a consequence of changes in the microenvironment that may affect all stem cells in the niche. Nevertheless, to model tumor development considering driver gene mutations to individual cells, we need to modify the model equations to include the mutants.

To show the framework for modeling tumor development induced by driver gene mutations, we consider the process with two types of mutations that increase the proliferation rate and decrease the differentiation rate (Figure 4A). Hence, let $Q_{i}(t,\mathbf{x})\ (i=0,1,2,3)$ represent the wild-type $(i=0)$ and the three mutant subpopulations $(i=1,2,3)$ cell counts, and $p_{i,j}$ represents the mutation rates. For the simplicity, we assume that gene mutations occur during cell division, and two daughter cells have the same mutant type. Therefore, equation (2) can be extended as follows:

[TABLE]

and

[TABLE]

Here, we consider only the driver mutation types, and at most one mutation occurs in each cell cycle, so that only the mutation rates $p_{0,1},p_{0,2}$ and $p_{1,3},p_{2,3}$ are nonzero value; and otherwise, the mutation rate $p_{i,j}$ is zero (Figure 4A).

Figure 4B-C shows the simulated dynamics. Single mutant cells occur prior to the obvious increase in the cell number, and the mutant cells eventually develop to double mutations that dominate the cell population (Figure 4B). Moreover, our simulation suggests that stemness increases with evolutionary processes when we limit the mutations to proliferation and differentiation (Figure 4C). Here, we consider only two types mutations that often occur in the precancerous stage [15, 27]. To simulate a more complicated process of cancer development, we must extend the simulation to include more mutations, such as apoptosis, DNA damage repair, and immune response pathways.

2.5 Modeling tissue growth with cell lineage dynamics

In the abovementioned models, we considered only cells capable of self-renewal, e.g., stem cells and progenitor cells. Nevertheless, to model tissue growth, we must include terminal differentiated cells that lose the ability to progress through the cell cycle. Therefore, let $Q(t,\mathbf{x})$ represent cells with self-renewal ability as previously mentioned, and let $P(t,\mathbf{x})$ represent the number of terminally differentiated cells. The terminally differentiated cells are produced from the stem and progenitor cells with the rate $\kappa(\mathbf{x})$ and cleared with the rate $\nu(\mathbf{x})$. Hence, equation 2 can be reformulated as follows:

[TABLE]

In the simulations shown in Figure 2, by considering the epigenetic state $0\leq x\leq 1$ as a stemness index and by distinguishing the stem cells from progenitor cells with the boundary $x=x_{0}$ (Figure S1), the numbers of stem cells, progenitor cells, and terminally differentiated cells can be determined as follows:

[TABLE]

This equation provides a model of multistage cell lineages shown in previous studies [1, 21, 46]. The tissue size is given by the total cell number as follows:

[TABLE]

and the distribution of stemness among all tissue cells is given by

[TABLE]

Figure 5 shows the simulated dynamics beginning with $100$ stem cells, which reveal the transition between stem and progenitor cells and the differentiation to terminally differentiated cells. Figure 5B shows the density of phenotypically different cell populations among the stem cells, progenitor cells and the terminally differentiated cells.

3 Discussion

Stem cell and progenitor cell regeneration is a basic cellular behavior associated with development, aging, and many complex diseases in multicellular organisms. In this study, to overlook the genetic details, we established a general mathematical framework to describe the process of stem cell and progenitor cell regeneration. This framework highlights cell heterogeneity and connects heterogeneity with cellular behaviors, e.g., proliferation, apoptosis, and differentiation/senescence. Cell heterogeneity is often associated with epigenomic markers that are subject to stochastic inheritance during cell division and is described by an inheritance probability function. Hence, the framework is a multiscale model that incorporates microscopic epigenetic state and gene expressions with macroscopic tissue growth through mesoscopic cell behaviors. We believe that this formula is helpful in answering the Weinberg question [86]. Despite the generality of this formula, different assumptions regarding the kinetic rate function and the inheritance probability can be applied to describe various biological processes related to stem cell regeneration (Figure 2).

In our framework, all stem and progenitor cells are described with a single compartment model, and different phenotypic cells are not distinguished explicitly. This approach differs from differentiation tree models that are widely used to describe the maintenance of hierarchically organized tissues. Recent experimental results have challenged the discrimination between stem and progenitor cell populations and have shown a continuous spectrum of results from cell differentiation [47, 65, 85]. Stochastic state transitions between different phenotypic cells lead to a dynamic equilibrium among a population of self-renewing cells [28, 52, 75]. Our model suggested that discrimination between cell types may not be necessary to describe tissue homeostasis. Different subtypes of cells can be characterized by their kinetic rates of proliferation, apoptosis, differentiation, senescence, etc. For convenience, these dynamic features are referred to as the kinotype, which is analogous to the genotype, epigenotype, and phenotype of a cell. The kinotype of cells is often associated with specific genes enriched in the related pathways. If the relationship between kinetic rates and the expression of these genes are known, we can extend the proposed framework to include the roles of specific genes. Therefore, in the future, we aim to develop a predictable model to investigate how variations in specific genes serve to alter the long-term dynamics of tissue growth.

Although the probabilistic epigenetic inheritance was considered, equation (2) is a deterministic equation that describes the dynamics of cell densities with different epigenetic states. This model often provides information regarding the average of multiple cells. To model a single cell, we must perform stochastic simulations that explicitly account for random events. Equation (2) suggests a numerical scheme of multiscale modeling for tissue growth where a multiple cell system is represented by a collection of epigenetic states in each cell as $\Omega_{t}=\left\{[C_{i}(\mathbf{x}_{i})]_{i=1}^{Q(t)}\right\}$. In each cell cycle, each cell undergoes proliferation, apoptosis, or terminal differentiation with a probability following the given kinetic rate so that both the system state $\Omega_{t}$ and the cell count $Q(t)$ change, and the epigenetic state of each cell undergoing cell division changes according to the predefined inheritance probability $p(\mathbf{x},\mathbf{y})$. In our previous study, this computational model was applied to model the process of inflammation-induced tumorigenesis and reproduced the two-stage tumorigenesis dynamics and revealed the competing oncogenic and onco-protective roles of inflammation. Based on the simulation results, which include the evolution of single-cell states, we were able to uncover the detailed process of cancer development.

Material and methods

Resource

Source MATLAB code for the study is available from

https://github.com/jzlei/StemCell.

Age-structured model and delay differential equation models

In the G0 cell cycle model, $Q(t)$ is the number of resting-phase stem cells, $s(t,a)$ is an age-structured quantify to represent the population of proliferating stem cells, and the age $a=0$ is their time of entry into the proliferative state. The resting-phase cells can either reenter the proliferative phase at a rate $\beta(Q)$ or differentiate into downstream cell lines at a rate $\kappa$. The proliferating stem cells are assumed to undergo mitosis at a fixed time $\tau$ after entry into the proliferating compartment and to be lost randomly at a rate $\mu$ during the proliferating phase. Each normal cell generates two resting-phase cells at the end of mitosis. Here, the units of cell population are often measured by the number of cells per unit body weight, e.g., $\mathrm{cells/kg}$, and the rates of proliferation, differentiation, and apoptosis are often united with $\mathrm{day}^{-1}$.

The above assumptions yield the following partial differential equations [49]:

[TABLE]

Here, $\nabla=\partial/\partial t+\partial/\partial a$. The boundary condition at $a=0$ is as follows:

[TABLE]

and the initial conditions are

[TABLE]

Equations (9)-(11) provide a general age-structured model of homogeneous stem cell regeneration.

By integrating (9)-(10) with the characteristic line method, we obtain the following close-form differential equation [49]:

[TABLE]

where $Q_{\tau}=Q(t-\tau)$. When we consider only the long-term behavior $(t>\tau)$ and shift the original time point to $\tau$, the delay differential equation model is as follows

[TABLE]

This equation describes the general population dynamics of homogeneous stem cell regeneration.

Formulation of the proliferation rate

The effect of feedback regulation from the cell population to the proliferation rate is given by the function $\beta(Q)$. Biologically, the self-renewal ability of a cell is determined by both microenvironmental conditions, e.g., growth factors and various types of cytokines, and intracellular signaling pathways, e.g., growth factor receptors and cell cycle checkpoints, such as fibroblast growth factors (FGFs) and the transforming growth factor beta (TGF-$\beta$) family [58, 63, 67]. The exact activation pathways that regulate the self-renewal of stem cells are poorly understood. Here, we derived a phenomenological formulation based on simple but general assumptions.

There are positive and negative signals for stem cell proliferation. We assume that positive growth factors are secreted by the niche, and growth factor inhibitors are released by the cells. Different types of cytokines bind to the cell surface receptors to regulate cell behavior. Let $[\mathrm{L}]$ denote the concentration of ligands for growth factor inhibitor; $[\mathrm{R}]$, the density of free receptor; $[\mathrm{R}^{*}$ ], the density of activated receptor; $Q$, the stem cell number. The total number of receptors is

[TABLE]

where m is the average number of receptors per cell. If n ligands are required to activate one receptor, we assume that ligands bind to the receptor following the law of mass action as follows:

[TABLE]

At equilibrium, we have the following equation:

[TABLE]

where $K$ is the equilibrium constant. We assume that the activated receptors inhibit cell proliferation so that the proliferation rate β is proportional to the fraction of free receptors on a cell as follows:

[TABLE]

From (14)-(15), we obtain the following expression:

[TABLE]

When ligands are secreted from stem cells and are cleared at a constant rate, the ligand concentration is proportional to the cell number, which gives $[\mathrm{L}]=\sigma Q$. Thus, we obtain the final form of the proliferation rate as follows:

[TABLE]

where $\theta=(1/\sigma)^{1/n}$ is the $50\%$ effective coefficient (EC50).

From (16), the proliferation rate, which is important for tissue homeostasis, approaches [math] due to the antiproliferation signals when the cell number $Q$ is sufficiently large. However, the capabilities of self-sufficiency in growth signals and insensitivity to antigrowth signals are the two characteristics of cancer that enable malignant tumor cells to escape antigrowth signals [31]. Hence, to model tumor development, the proliferation rate can be modified to include a nonzero constant $\beta_{1}$ for self-sustained growth signals as follow:

[TABLE]

Steady state of the G0 cell cycle model and oncogenic signaling pathways

From (13), the steady state $Q(t)\equiv Q^{*}$ is given by the equation

[TABLE]

which yields either $Q^{*}=0$, or

[TABLE]

When $\beta(Q)$ is given by (17), equation (13) has a unique positive steady state if and only if

[TABLE]

In particular, when

[TABLE]

(13) has only a zero steady state, and the zero solution is unstable. Hence, all positive solutions approach infinity, which corresponds to uncontrolled growth. Therefore, the inequality (19) summarizes a general condition to have uncontrolled growth, i.e., malignant tumors. Biologically, equation (19) is satisfied if there is self-sufficiency in growth signals and/or insensitivity to antigrowth signals (increasing $\beta_{1}$), evasion of apoptosis (decreasing $\mu$), and dysregulation in the differentiation and/or senescence pathways (decreasing $\kappa$), which are well known hallmarks of cancer [31].

Age-structured model of heterogeneous stem cell regeneration

When heterogeneity in stem cells is considered, and assumed that the apoptosis rates of cells during cell division are dependent on the epigenetic state of the cell before entering the proliferating phase, the age-structured model (9) becomes

[TABLE]

and

[TABLE]

While we considered the epigenetic state $\mathbf{x}$ in the first equation as a parameter, the characteristic line method remains valid, which gives the following equation (here we show only the result of long-term behavior):

[TABLE]

Substituting $s(t,\tau(\mathbf{x}),\mathbf{x})$ into the second equation, we obtained the following equation:

[TABLE]

which gives the equation (2) for heterogeneous stem cell regeneration.

Let

[TABLE]

which is the total cell number, and integrate (Age-structured model of heterogeneous stem cell regeneration) with $\mathbf{x}$. Notably, when $p(\mathbf{x},\mathbf{y})d\mathbf{x}=1$, we obtain the following equation:

[TABLE]

Define

[TABLE]

as the density of cells with a given epigenetic state $\mathbf{x}$, then

[TABLE]

Hence, we obtained the equation

[TABLE]

The inheritance probability $p(\mathbf{x},\mathbf{y})$

In (2), the inheritance probability $p(\mathbf{x},\mathbf{y})$ is essential to describe the heterogeneity of cells. The function $p(\mathbf{x},\mathbf{y})$ is associated with the process of cell division during which the epigenetic code and molecules distribute to daughter cells through complex regulation mechanisms that are not well understood. Hence, the general mathematical formula of the function $p(\mathbf{x},\mathbf{y})$ remains unknown. Here, we proposed an attempt to define the function based on the random inheritance of histone modifications.

In eukaryotic cells, most DNA sequences are enclosed in nucleosomes in which DNA sequences wrap around a histone octamer that is composed of one $\mbox{(H3-H4)}_{2}$ tetramer capped by two H2A-H2B dimers. These histones can undergo diverse posttranslational covalent modifications that lead to either active or repressive gene expression activities [5, 41, 44]. The patterns of histone modification dynamically change over time, and hence define a dynamic histone code for the transcription activity. The dynamics of histone modifications consist of complex process, including nucleosome assembly, writing and erasing of the modification markers, and random inheritance during DNA replication [71, 76]. Detailed computational models for the process of histone modification and random inheritance over the cell cycle remain a challenging issue in computational biology. While we consider the main process of writing and erasing the modification markers that are modulated by the related enzymes, the kinetics of histone modification can be modeled through stochastic simulations [37, 42].

In a proposed dynamic model of histone modification [37, 42], bivalent modifications of the histone H3, the trimethylation of H3 lysine 4 (H3K4me3) and the trimethylation of H3 lysine 27 (H3K27me3), were considered. Each H3 histone can be in one of the following states: unmodified (U), modified by the activating marker H3K4me3 (A), or modified by the repressing marker H3K27me3 (R). Each nucleosome can be in one of six physically nucleosome states, which include UU, AU, UR, AA, AR, and RR. The nucleosome states dynamically change according to methylation/demethylation, which are regulated by corresponding enzymes. During DNA replication, parental histones and newly synthesized histones are randomly distributed on daughter strands. To avoid the dilution of histone markers, maintenance modifications in the new histones can be achieved by using a neighboring histone as a template [71]. Hence, writing enzyme activities are dependent on the states of neighboring nucleosomes. Thus, changes in the nucleosome state over the cell cycle due to the random distribution of histone markers during DNA replication and kinetic methylation/demethylation can be tracked with a stochastic simulation [37].

Based on the abovementioned model simulation, we are able to study how the nucleosome state of daughter cells depends on that of the mother cells. For example, considering a DNA segment with $N$ nucleosomes, we counted the number ($N_{A}$) of nucleosomes with active markers (either AA or AU) at each cycle. The simulation results suggested that considering the state of the mother cell, the active nucleosome number is a random number with a binomial distribution with the parameter (success probability $p$) dependent on the state of the mother cell as follows [37]:

[TABLE]

Considering the nucleosome state through the fraction ($f_{A}=N_{A}/N$) of active nucleosome numbers, we can extend the binomial distribution to a continuous probability distribution defined on the interval $[0,1]$, which is given by the following beta distribution:

[TABLE]

Hence, for the specific situation of the random inheritance of histone modifications during the cell cycle, we can use the beta distribution probability as the inheritance probability function $p(x,y)$. Here, we extend this formulation to general cases and proposed (5) for the inheritance functions.

Beta distribution

The beta distribution is a family of continuous probability distributions defined on the interval $[0,1]$ and is parameterized by two positive shape parameters that appear as exponents of the random variable and that control the shape of the distribution. The probability density function (PDF) of the beta distribution, for $0\leq x\leq 1$, and the shape parameters $a,b>0$, is a power function of the variable $x$ and of its reflection $(1-x)$ as follows:

[TABLE]

where $\Gamma(z)$ is the gamma function.

For a random variable $X$ beta-distributed with parameters $a$ and $b$, which is denoted by $X\sim\mathrm{beta}(a,b)$, the mean and variance are as follows:

[TABLE]

Then, it is easy to obtain

[TABLE]

Hence, if we assume

[TABLE]

then

[TABLE]

which gives

[TABLE]

This gives equation (6) to determine the shape parameters from the functions $\phi_{i}(\mathbf{y})$ and $\eta_{i}(\mathbf{y})$.

Simulations for stem cell regeneration

Here, we present a simple example to show the numerical scheme to simulate stem cell regeneration based on the proposed model equations.

We consider a situation with one epigenetic state $x$ ($0\leq x\leq 1$) that affects only cell proliferation and differentiation so that only the rates $\beta$ and $\kappa$ are dependent on the epigenetic state $x$. Therefore, we have the following model equation:

[TABLE]

Here, $\xi(x)=1$ so that

[TABLE]

To specify the rate functions $\beta$ and $\kappa$, we assume that the state $x$ affects the proliferation and differentiation rates in a manner similar to the stemness so that large value of $x$ indicates the stem cells with a low proliferation rate, an intermediate value $x$ indicates progenitor cells with a high proliferation rate, and the terminated differentiation rate $\kappa$ is a decreasing function of $x$ that approaches zero when $x$ is large. This is mathematically expressed as follows:

[TABLE]

and

[TABLE]

The inheritance probability function $p(x,y)$ is defined from the beta distribution density function as with predefined function $\phi(y)$ and $\eta(y)$ as follows:

[TABLE]

and

[TABLE]

Figure 7 shows the functions $\beta_{0}(x)$, $\kappa(x)$, and $p(x,y)$.

In the simulations shown in Figure 2, the parameters are as follows: $\mu=2.0\times 10^{-4}\,\mathrm{h}^{-1},\tau=20\,\mathrm{h},\theta=10^{3}\mathrm{cells},a_{1}=5.8,a_{2}=2.2,a_{3}=3.75,b_{1}=4.0,\eta(x)=60$, and

• •

in (B)-(C): $\bar{\beta}=0.12\,\mathrm{h}^{-1},\kappa_{0}=0.02\,\mathrm{h}^{-1},\phi(x)=0.08+1.0\times\frac{(1.65x)^{1.8}}{1+(1.65x)^{1.8}}$;

• •

in (D)-(E): $\bar{\beta}=0.12\,\mathrm{h}^{-1},\kappa_{0}=0.02\,\mathrm{h}^{-1},\phi(x)=0.04+0.9\times\frac{(1.65x)^{1.8}}{1+(1.65x)^{1.8}}$;

• •

in (F)-(G): $\bar{\beta}=0.12\,\mathrm{h}^{-1},\kappa_{0}=0.0002\,\mathrm{h}^{-1},\phi(x)=0.08+1.2\times\frac{(1.65x)^{1.8}}{1+(1.65x)^{1.8}}$.

In the simulation show in Figure 3, $\bar{\beta}=4.8\mathrm{h}^{-1},\kappa_{0}=0.4\mathrm{h}^{-1}$, and other parameters are the same as those shown in Figure 2B.

In the simulation shown in Figure 4, we set the wild-type cell parameters according to those shown in Figure 2B, and $p_{0,1}=p_{0,2}=0.5\times 10^{-4},p_{1,3}=p_{2,3}=1\times 10^{-4}$. For mutation 1, the proliferation rate is twice of wild-type cells, and for mutation 2, the differentiation rate is $1/10$ of wild-type cells.

In the simulation shown in Figure 5, we set $x_{0}=0.7$, the differentiation rate increases from [math] to normal level following $\kappa(t)=0.04\times(\frac{1}{2}-\frac{1}{1+e^{t/10000}})$, and the other parameters are the same as those shown in Figure 2B.

Acknowledagement

This work was support by the National Natural Science Foundation of China (NSFC 91730301 and 11831015).