Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Math Model Nat Phenom. Author manuscript; available in PMC 2010 July 9.
Published in final edited form as:
Math Model Nat Phenom. 2009 June 1; 4(4): 103–117.
doi:  10.1051/mmnp/20094403
PMCID: PMC2900806

Multi-scale Modelling for Threshold Dependent Differentiation

A. Q. Cai,1,2,* Y. Peng,1,2 J. Wells,3 X. Dai,3 and Q. Nie1,2


The maintenance of a stable stem cell population in the epidermis is important for robust regeneration of the stratified epithelium. The population size is usually regulated by cell secreted extracellular signalling molecules as well as intracellular molecules. In this paper, a simple model incorporating both levels of regulation is developed to examine the balance between growth and differentiation for the stem cell population. In particular, the dynamics of a known differentiation regulator c-Myc, its threshold dependent differentiation, and feedback regulation on maintaining a stable stem cell population are investigated.

Keywords: cell division, cell differentiation, stem cells, threshold, cell population balance model, c-Myc, epidermal development

1. Introduction

The fate of a cell, whether it dies, divides or differentiates, is usually governed by its gene regulatory network interpreting external stimuli and internal regulations. The robust maintenance of a cell type population size, i.e. homeostasis, requires complex intracellular and extracellular regulations. When homeostasis is lost, a cell population either grows unboundedly or loses its ability to reach optimal size. Consequently, birth defects or diseases such as cancer, may occur. Epithelial development in the skin epidermis is an excellent model system for studying homeostasis maintained by dynamic interactions of feedback and adaptation occurring at scales involving both the cell population and individual cells.

The epidermis that undergoes constant regeneration has been traditionally proposed to consist of three types of cells: stem, transit amplifying (TA), and terminally differentiated cells. The stem cells divide into stem cells and/or differentiate into TA cells. The TA cells are distinguished from the self-renewing stem cells by a limited proliferative potential, as they can only undergo limited rounds of division before becoming terminally differentiated cells [8]. The proto-oncogene c-Myc, an intracellular transcription factor, plays a paradoxal role in epidermis by promoting both proliferation and terminal differentiation [23]. When c-Myc levels are elevated in post-mitotic cells, these cells re-enter cell cycle to become proliferation-competent TA cells [18]. However, when c-Myc levels are aberrantly elevated in stem/TA cells using genetic means, this leads to increased terminal differentiation and depletion of the stem cell pool [2, 21]. A solution to the paradox seems to reside in the ability of c-Myc to drive stem cells into a TA phase, thereby initiating the irreversible process of terminal differentiation. Together, these studies have led to the hypothesis that c-Myc levels are normally higher in TA cells compared to stem cells [18]. In addition, cells occupying different lineage stages produce secreted factors (extracellular molecules), such as TGF-β. Such secreted factors are known to regulate intracellular molecules including c-Myc [9]. Thus, the cell populations are clearly determined by both intracellular and extracellular regulation and their interactions. As the populations of different cell lineage stages have different levels of c-Myc and the fate of individual cells are c-Myc dependent, an understanding of c-Myc dynamics and distribution is important for studying epidermal stem cell dynamics. Such knowledge may have general implications for other multi-stage cell lineages, for example that in the hematopoietic system where c-Myc is also found to promote differentiation [25].

To study the multi-scale interaction on stem cell populations, we present in this paper a simple cell population model incorporating the level of c-Myc, its potential extracellular regulations, and c-Myc threshold dependent cell differentiation. Unlike a classic logistic model [16] in which the steady state size of the cell population is limited by the carrying capacity, a generic parameter used to model the external constraints on growth, such as cell-cell competition for nutrients, spatial effects, and growth inhibition via contact [1], the new model reaches steady state through intracellular and extracellular interactions and regulations. This model has many similar features of population balance models [10], also referred to as structured models. Such models have been used for biological populations in terms of time and an intrinsic variable. The intrinsic variable of the population may correspond to cell maturity or age [4, 22], cell size [3], the concentration of a cell division labelling dye [13], or gene expression [15]. Our model in this paper emphasizes the study of the cell population as a function of c-Myc, a key intracellular regulator of cell proliferation and differentiation. In this regard, it is different from some other discrete (e.g. [7, 20]), continuous (e.g. [6, 11, 12]), and stochastic (e.g. [14]) cell-lineage models.

The paper is organized as follows: in Section 2, we present the mathematical model; in Section 3, we explore and discuss the limitations of a linear case; in Section 4, we study a nonlinear model with regulated threshold dependent differentiation; and in Section 5, we summarize our results and discuss possible extensions of the model.

2. Mathematical Model

We model the stem cell population with c-Myc number threshold dependent differentiation. Let us define C(t) as the number (or concentration) of c-Myc protein within individual stem cells. The maximum number (or concentration) of c-Myc within a cell is cmax. The exact number of c-Myc molecules in a cell is unknown, however, a transit amplifying cell is likely to have more c-Myc than a stem cell. In the model, cmax, the maximum number of c-Myc proteins in a cell, is chosen with an arbitrary unit. The quantities in the model related to c-Myc are all scaled by cmax, and C = cmaxc. Let N (c, t) be the stem cell population distribution in c. The typical epidermal cell division time is approximately 24 hours [24]. We chose the time unit as one day.

Then the model takes the form,


where H() is the Heaviside function and cSC is the threshold value for cell differentiation. The natural boundary conditions for this system are N = 0 at c = 0 and c = 1. The parameters for cell division and differentiation are τdiv and τdiff, representing the mean times to divide or differentiate, respectively. The second term on the left-hand side can be interpreted as population flux with respect to c-Myc. The change of the cell population distribution is dependent on the change of c-Myc, as well as cell division and differentiation. The number of cells with c between a0 and a1 is clearly a0a1Ndc, so the total stem cell number at time t is


The 22N(2c, t) term models symmetric cell divisions where the c-Myc level in the parent cell is evenly distributed to the daughter cells after each division. A factor of 2 is the doubling in cell number and the other factor of 2 is due to the different interval sizes occupied by the daughter cells in (c, c + Δc) and the mother cell in (2c, 2(c + Δc)).

The model (2.1) is solved numerically with a Matlab hyperbolic solver [19] using a Lax-Wendroff scheme. This numerical scheme is second order accurate in time and space, which corresponds to the c variable in this model. An independent implementation of the second order Lax-Wendroff scheme is also carried out for testing the Matlab solver.

3. Non-regulated differentiation: a linear case

The classical Malthusian growth model shows that an unregulated population exhibits either exponential growth or decay [16]. When the net growth rate of the population is zero, the total population size should be fixed for all time. However, the distribution of a population may change over time, even when the total population size is fixed. In this section, we consider a case of population growth without any regulation. Feedback regulation by the cell population on c-Myc dynamics, i.e. dcdt=f(N,c), is studied in Section 4.

Consider the growth of a population where differentiation can occur for all c without any threshold, i.e. cSC = 0. We assume that protein c is produced at a constant rate without significant degradation, dcdt=c¯>0. The model (2.1) then reduces to


Analytical approaches for population balance equations (3.1) can be found in [10]. The total cell number calculated by integrating Eq. (3.1) over C is


and it is Malthusian growth.

For such a system, only balanced cell division and differentiation, τdiv = τdiff, can lead to a steady state population distribution. We calculate and use the analytical steady state distribution to test the accuracy of the numerical solution. While it appears c-Myc could be unbounded as it is increasing at a constant rate, the c-Myc level inside the cell is in fact halved after each cell division and c-Myc within the cell population remains bounded.

When the cell division rate and the differentiation rate are different, the cell population grows or decays exponentially depending on the relative strength of the two rates. When cell differentiation occurs more frequently than division, i.e. τdiff < τdiv, the population eventually goes extinct even though the cell population distribution centered at c=12 divides and differentiates over time, and successive rounds of cell division create transient peaks in the population distribution (Figure 1). The total cell number of the numerical solution is consistent with the analytical solution given by Eq. (3.2) (Figure 1(c)). When cell division is more frequent than differentiation τdiv < τdiff, the population blows up over time (Figure 2). In both cases, the initial cell distribution is

Figure 1
Extinction of a cell population, c = 0.1, τdiv = 1, τdiff = 0.5. (a) The cell population distribution over time. (b) The cell population distribution at different times. (c) The total cell number over time; Solid line: ...
Figure 2
Unbounded growth of a cell population. c = 0.1, τdiv = 1, τdiff = 2. (a) The cell population distribution over time. (b) The cell population distribution at different times. (c) The total cell number over time. Solid line: ...

When cell division and differentiation are equal (τdiv = τdiff [equivalent] τ), the total population size remains unchanged over time while the cell population distribution varies over time and reaches a non-uniform steady state. Such detailed dynamical behavior could not be captured by classical homogeneous population growth models. In particular, for this case, the governing equation for N(c) becomes a linear system,


One may seek an analytical solution to the steady state with a Dirichlet series [5],


Then the recurrence relation dndn1=2(12n) is obtained, with a=2c¯τ, and the steady state solution becomes


where the constant d0 can be expressed in terms of the total cell number N(0),


Figures 3 and and44 show the approach towards the same steady state solution when two different initial cell distributions with the same total cell number are used. In Figure 3, Eq. (3.3) is used as the initial condition while in Figure 4 we use N(c,0)=10erf(210)erf(10/2)exp(40(c1/2)2) as the initial condition. The exact steady state distributions calculated using (3.6) are shown in Figure 3(b) and Figure 4(b) in red. It is interesting to observe that both populations starting with one peak exhibit multiple transient peaks before reaching the one-peak steady state.

Figure 3
Balanced cell division and differentiation; c = 0.1, τ = 1. (a) The population distribution over time. (b) The analytical solution of the steady state distribution is in red. The cell population at different times. (c) The absolute ...
Figure 4
Balanced cell division and differentiation; c = 0.1, τ = 1. (a) The cell population distribution over time. (b) The analytical solution of the steady state distribution is in red. (c) The absolute error of the numerical result ...

As seen in Figure 3(c) and Figure 4(c), the maximum absolute error between the exact analytical steady state solution and the numerical solution at t = 10 for successive halving of ΔC is roughly 0.25, 0.07, and 0.025, and 0.25, 0.06, 0.02, respectively. The maximum error is reduced by a factor of four when the spatial resolution increases by a factor of two. This suggests the numerical calculation is second-order accurate.

These calculations also indicate that without any regulation on cell population growth, a steady state population can only be maintained with equal cell division and differentiation rates for all time. The same steady state population distribution can be reached from different initial population distributions (i.e., cells expressing different levels of c-Myc) as long as initial total cell numbers are the same. In addition, (3.6–3.8) show the total initial cell number only affects the height of the steady state, not the distribution of the steady state, and the steady state is only controlled by cτ, a product of the mean growth rate for c-Myc and the differentiation and division rate.

4. Regulated threshold dependent differentiation

In this section we consider the full model, where cell differentiation is possible only when c-Myc inside a cell is above the threshold cSC.

After integrating the full model (2.1), the total cell number at steady state satisfies the following relationship:


Eq. (4.1) suggests that τdiff < τdiv is a necessary condition for existence of a steady state. This is in contrast to the case without differentiation threshold in Section 3. where the division rate must equal the differentiation rate for maintaining steady state population size.

With a differentiation threshold, there is only cell division in the region [0, cSC]. To counter this extra growth for achieving steady state, it is necessary for the cells with c above the threshold to be more likely to differentiate than to divide. In order to maintain stable population size, intuitively, we expect a regulation mechanism to bias cell differentiation when population size is large, and bias division when the population is small. Such regulation could be achieved via diffusible secreted factors and may occur directly by regulating c-Myc, or indirectly through other genes such as Ovol1 which regulates c-Myc [17]. It is natural to assume that the amounts of the diffusive cell-secreted factors are proportional to the total cell numbers. Here, we consider a simple positive regulation on the synthesis of c through a Hill function:


where c0 < c1. The positive regulation is to bias differentiation as the population size increases. Increasing population size leads to increased c-Myc synthesis via putative mechanisms mentioned above, thus intracellular c-Myc will reach the differentiation threshold faster. The parameter γ0 is a reference value for which the regulation becomes significant. The c-Myc synthesis rate is c0 when the stem cell population is relatively small N [double less-than sign] γ0 and c1 when N [dbl greater-than sign] γ0. The Hill exponent, n0, measures the slope of transition between the two states: c0 and c1.

There are two asympotic limits for (4.2). When the population size is relatively small, i.e. N [double less-than sign] γ0, then c approaches the value c0/cdeg. As the population size increases, the regulatory effect of the population becomes more significant when N [dbl greater-than sign] γ0 and c-Myc approaches the value c1/cdeg. It should be expected that when c1/cdeg < cSC the cell population will grow out of control because the c-Myc level within cells never increases above the differentiation threshold. However, the condition c1/cdeg > cSC can not guarantee existence of a steady state. As seen in Figure 5, both cases have c1/cdeg > cSC but when regulation strength is not strong enough the population grows exponentially within the time frame under examination. When the regulation becomes stronger, the population reaches a steady state.

Figure 5
The population may grow exponentially when the regulation strength is not strong enough. Solid line: c1 = 0.6, a case of exponential growth; “+”: c1 = 0.8, a case of existence of steady state. Other parameters are τdiv = 1, τ ...

Unlike the non-regulated population growth model presented in Section 3. the steady state population distribution through regulation is independent of the initial cell population size. Figures 6 and and77 are simulations of populations with different initial distributions and initial cell numbers, and all the populations reach the same steady state (Figures 7(b)). In the six cases of Figure 6, the total numbers of stem cells at the steady state are smaller than the initial total cell numbers for the first three cases, and it is the opposite for the second three cases.

Figure 6
Temporal dynamics for six different initial populations: (a) N(c, 0) = 20 exp(−40(c − 1)2), (b) N(c, 0) = 40 exp(−40c2), (c) N(c, 0) = 40 exp(−40(c − 1)2), (d) N(c, 0) = 2.5 exp(−40(c − 1)2), (e) ...
Figure 7
The systems with different initial population distributions all approach the same steady state. (a) The total cell number over time for the six different initial conditions. (b) The population distributions at t = 12 for the six cases shown in Figure ...

In these calculations, c0/cdeg < cSC < c1/cdeg. The regulation causes c to approach a value within the exclusive cell division region when the population size is small relative to γ0, and c to approach a value greater than the differentiation threshold otherwise. Transient proliferation can be seen where the initial distribution is concentrated below the differentiation threshold of c < cSC. Most notably in Figure 6(b), the cell population reaches a maximum before declining towards the steady state. The parameter values used for Figures 6 and and77 are τdiv = 1, τdiff = 0.2, c1 = 1, c0 = 0.2, γ0 = 2, cdeg = 1, n0 = 1, cSC = 0.4.

Finally, we study the effect of the different parameter values on the steady state distribution in Figure 8. When differentiation rate τdiff becomes larger, the total steady state population size becomes larger as well (Figure 8(a)) because the time to differentiation is longer. Smaller c0 values usually lead to a larger population at the steady state (Figure 8(b)). This is because slower c-Myc synthesis results in a longer time to reach the differentiation threshold, hence, more proliferation occurs before reaching steady state. The parameter γ0 measures the population size for which the feedback regulation becomes significant (i.e., large population size resulting in higher c-Myc expression). For larger γ0, the steady state population size is greater (Figure 8(c)) because the population grows to a larger size before being constrained by the regulation.

Figure 8
Steady state profiles at different parameter values. (a) τdiff = 0.1, 0.2, 0.4; (b) c0 = 0.05, 0.2, 0.6; (c) γ0 = 1, 2, 4. The other parameters are τdiv = 1, c1 = 1, cdeg = 1, n0 = 1, cSC = 0.4.

5. Discussion

We have studied a simple model of cell population growth with threshold dependent differentiation. One of the important features in the model is a multi-scale regulation between the putative extracellular secreted factors and intracellular molecules. Such regulation along with the threshold dependent differentiation is found to be critical for maintaining homeostasis of the cell population. A necessary condition for the existence of a steady state is the mean time to differentiation is less than the mean time for division. In addition, various conditions for creating a steady state cell population and dependence of the cell population on parameters have been discussed.

It would be of interest to further study the nature of the steady states systematically. Computationally, we have observed that different initial cell numbers (Figure 6 and and7)7) with different mean numbers of c-Myc approach the same steady state in some regions of parameter space. However, we also observed non-existence of steady state solutions for some other parameters (Figure 5). One question at hand is to investigate whether the steady state solution is a global attractor using an analytical approach.

In the current model, we have only considered a threshold for stem cell differentiation; there may very well be a threshold for cell division. Because quiescent stem cells tend to have very low levels of c-Myc [23], we might use a lower threshold on c for cell division. Also, we have observed that the stem cell population in our current model tends to either reach homeostasis or grow without control for most of the parameter values, with very few cases of vanishing population. This suggests that the simple model has a simple form that prevents the stem cell population from going to extinction. It would be interesting to study the growth of stem cells when the transit amplifying cells, whose cell cycle exit and differentiation might also be threshold dependent, are included in the model. Such extensions to multi-stage cell lineage should be straightforward. As more detailed and realistic intracellular controls through gene regulatory networks are built into such models, their study will provide better insights for the role of each regulatory component and the system behavior of the multi-stage cell lineage system.


The authors would like to thank Dr Ching-Shan Chou for discussions of numerical schemes. This work was supported by the NSF/NIH initiative on Mathematical Biology through R01GM75309 and R01GM67247 from the National Institute of General Medical Sciences, NIH P50GM76516, R01AR47320, and R01AR47320-08S1 grants.


1. Abercrombie M. The crawling movement of metazoan cells. Proc R Soc Lond B Biol Sci. 1980;207(1167):129–147.
2. Arnold I, Watt FM. c-Myc activation in transgenic mouse epidermis results in mobilization of stem cells and differentiation of their progeny. Curr Biol. 2001;11(8):558–568. [PubMed]
3. Basse B, Baguley BC, Marshall ES, Joseph WR, van Brunt B, Wake GC, Wall DJN. A mathematical model for analysis of the cell cycle in cell lines derived from human tumors. J Math Biol. 2003;47(4):295–312. [PubMed]
4. Bernard S, Pujo-Menjouet L, Mackey MC. Analysis of cell kinetics using a cell division marker: Mathematical modeling of experimental data. Biophys J. 2003;84:3414–3424. [PubMed]
5. van Brunt B, Wake GC, Kim HK. On a singular Sturm-Liouville problem involving an advanced functional differential equation. European J Appl Math. 2001;12:625–644.
6. Cai AQ, Landman KA, Hughes BD, Witt CM. T cell development in the thymus: From periodic seeding to constant output. J Theor Biol. 2007;249(2):384–394. [PubMed]
7. Clayton E, Doupé DP, Klein AM, Winton DJ, Simons BD, Jones PH. A single type of progenitor cell maintains normal epidermis. Nature. 2007;446:185–189. [PubMed]
8. Fuchs E, Raghavan S. Getting under the skin of epidermal morphogenesis. Nat Rev Genet. 2002;3:199–209. [PubMed]
9. Glick AB, Kulkarni AB, Tennenbaum T, Hennings H, Flanders KC, O’Reily M, Sporn MB, Karlsson S, Yuspa SH. Loss of expression of transforming growth factor β in skin and skin tumors is associated with hyperproliferation and a high risk for malignant conversion. Proc Natl Acad Sci USA. 1993;90:6076–6080. [PubMed]
10. Hjortsø MA. Population balances in biomedical engineering: Segregation through the distribution of the cell states. McGraw-Hill; 2006.
11. Johnston MD, Edwards CM, Bodmer WF, Maini PK, Chapman SJ. Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. Proc Natl Acad Sci USA. 2008;104(10):4008–4013. [PubMed]
12. Lo WC, Chou CS, Gokoffski KK, Wan FYM, Lander AD, Calof AL, Nie Q. Feedback regulation in multistage cell lineages. Math Biosci Eng. 2008;6(1):59–82. [PMC free article] [PubMed]
13. Luzyanina T, Roose D, Schenkel T, Sester M, Ehl S, Meyerhans A, Bocharov G. Numerical modelling of label-structured cell population growth using CFSE distribution data. Theor Biol Med Model. 2007;4(26) [PMC free article] [PubMed]
14. Mangel M, Bonsall MB. Phenotypic evolutionary models in stem cell biology: replacement, quiescence, and variability. PLoS one. 2008;3(2):e1591. [PMC free article] [PubMed]
15. Mantzaris NV. Single-cell gene-switching networks and heterogeneous cell population phenotypes. Comput Chem Eng. 2005;29:631–643.
16. Murray JD. Mathematical biology. Vol. 1. New York: Springer; 2002.
17. Nair M, Teng A, Bilanchone V, Agrawal A, Li B, Dai X. Ovol1 regulates the growth arrest of embryonic epidermal progenitor cells and represses c-Myc transcription. J Cell Biol. 2006;173(2):253–264. [PMC free article] [PubMed]
18. Pelengaris S, Littlewood T, Khan M, Elia G, Evan G. Reversible activation of c-Myc in skin: induction of a complex neoplastic phenotype by a single oncogenic lesion. Mol Cell. 1999;3(5):565–577. [PubMed]
19. Shampine LF. Solving hyperbolic PDEs in Matlab. Appl Num Anal Comp Math. 2005;2(3):346–358.
20. Tomlinson IP, Bodmer WF. Failure of programmed cell death and differentiation as causes of tumors: some simple mathematical models. Proc Nat Acad Sci USA. 1995;92:11130–11134. [PubMed]
21. Waikel RL, Kawachi Y, Waikel PA, Wang XJ, Roop DR. Deregulated expression of c-Myc depletes epidermal stem cells. Nat Genet. 2001;28(2):165–168. [PubMed]
22. Wang G. Estimation of the proliferation and maturation functions in a physiologically structured model of thymocyte development. J Math Biol. 2007;54:761–786. [PubMed]
23. Watt FM, Frye M, Benitah SA. Myc in mammalian epidermis: how can an oncogene stimulate differentiation? Nat Rev Cancer. 2008;8:234–242. [PMC free article] [PubMed]
24. Willie JJ, Jr, Pittelkow MR, Shipley GD, Scott RE. Integrated control of growth and differentiation of normal human prokeratinocytes cultures in serum-free medium: clonal analyses, growth kinetics and cell cycle studies. J Cell Physiol. 1984;121:31–44. [PubMed]
25. Wilson A, Murphy MJ, Oskarsson T, Kaloulis K, Bettess MD, Oser GM, Pasche AC, Knabenhans C, MacDonald HR, Trumpp A. c-Myc controls the balance between hematopoietic stem cell self-renewal and differentiation. Genes Dev. 2004;18:2747–2763. [PubMed]