|Home | About | Journals | Submit | Contact Us | Français|
Mathematical modeling is being increasingly recognized within the biomedical sciences as an important tool that can aid the understanding of biological systems. The heavily regulated cell renewal cycle in the colonic crypt provides a good example of how modeling can be used to find out key features of the system kinetics, and help to explain both the breakdown of homeostasis and the initiation of tumorigenesis.
We use the cell population model by Johnston et al.5 to illustrate the power of mathematical modeling by considering two key questions about the cell population dynamics in the colonic crypt. We ask: how can a model describe both homeostasis and unregulated growth in tumorigenesis; and to which parameters in the system is the model most sensitive? In order to address these questions, we discuss what type of modeling approach is most appropriate in the crypt.
We use the model to argue why tumorigenesis is observed to occur in stages with long lag phases between periods of rapid growth, and we identify the key parameters.
Mathematical modeling can be a powerful tool for understanding biologically observed phenomena which cannot be understood by verbal reasoning alone.1 One such example is that of homeostasis in the colonic crypt. The single layer of epithelial cells that line the crypt is renewed every two to three days by a number of long-living stem cells that remain at the bottom of the crypt.2 As the stem cells divide, their progeny migrate up the crypt wall, and once at the top they are shed into the lumen or undergo apoptosis. This general structure of stem, transit and differentiated cells is also found in many other biological systems, for example the hematopoietic system.
Many modeling studies have tried to capture this tightly regulated system, and three common approaches are to use compartmental, simulation or stochastic models (reviewed in ref. 3). Compartmental models include those by Tomlinson and Bodmer,4 Johnston et et al.,5 Boman et al.,6 Hardy and Stark,7 Wodarz8 and Paulus et al.,9 and simulation models have been presented by Loeffler et al.10,11 Models that capture the stochastic effects of mutational changes include those by Nowak et al.,12,13 Michor et al.,14-17 Komarova and Wang18 and Komarova.19,20
In this paper we focus on a compartmental model, and discuss mathematical techniques that can be used to address key questions relating to the population dynamics of the system. We will also show how this mathematical approach can suggest which parameters are most important to determine experimentally, and which it is less important to know accurately. We shall draw some examples from the recent paper by Johnston et al.5
We consider two biological questions: (1) how can we formulate a model that will allow for both the tight regulation of the number of cells in a healthy crypt, but also the breakdown of homeostasis when cancer occurs? (2) Which are the important parameters in this system? These questions are addressed in Examples 2 and 3 respectively. Firstly, we need to decide the most appropriate form of mathematical modeling for the cell populations, and this is discussed in Example 1. We begin, however, by summarizing our compartmental model5 which we will later use to address the key questions relating to colorectal cancer.
In our compartmental model,5 which is summarized in Figure 1, we followed Tomlinson and Bodmer4 who separated the cells in the crypt into populations of stem cells (denoted N0), fully-differentiated cells (denoted N2), and a third population for all the transit-amplifying cells that are in the process of differentiating, referred to as semi-differentiated cells (denoted N1). At the end of each cell cycle the stem cells can die, differentiate into transit cells or renew with probabilities a1, a2 and a3 respectively, whereas transit cells will die, differentiate into fully-differentiated cells or renew with probabilities b1, b2 and b3 respectively. The fully-differentiated cells will be removed from the system with probability c. The probabilities can also be thought of as the proportions of the cell populations that undergo each process. The cell cycle times for the stem and transit cell populations are denoted by t0 and t1 respectively and measured in hours.
We use this model to address the question of modeling homeostasis in Example 2. First, however, we discuss discrete, age-structured and continuous models and their relevance to crypt modeling.
What is the simplest and most appropriate technique for modeling the cell division process in the colonic crypt? In this section we first compare and contrast the discrete, age-structured and continuous modeling approaches, illustrating them by considering the stem cell population from Figure 1, and then explain when each approach is appropriate.
A discrete modeling approach follows the number of cells in a population at regular time intervals, and does not take into account where the individual cells are in their cell cycle. As in the model by Tomlinson and Bodmer,4 we count the number of stem cells N0(G) at each generation G measured in time units of the stem cell division time t0. The parameter a3 denotes the proportion of the stem cell population that renews at each division, assuming that all divisions are symmetric (asymmetric divisions can be thought of as half a symmetric renewal and half a symmetric differentiation). We relate the number of stem cells at any generation to the number in the previous generation by the difference equation
Assuming an initial number of stem cells, 0, this can be solved to find the population at each future Gth generation by the formula
An alternative approach is to use an age-structured model5 to find the number of stem cells N0(t,a) at a continuous time point t and age a. In this context, the “age” of a cell refers to the length of time since its last division. This model consists of partial differential equations (PDEs) that quantify the rate of change of N0 as both time t and age a vary [see equation (3)]. For example, assuming (as we did in ref. 5) that the cells only die, differentiate or renew at the end of their cell cycles, we obtain the PDE
which is valid for all ages 0<a<t0. At the end of each cell cycle we also have that
which is a more specific form of equation (1) relating the numbers of cells being born and those reaching the end of their cell cycle at time t. If an initial distribution of the ages of cells n0(a) is known, the solution is given by
where G is the number of times that cells of age a have been through the cell cycle at time t.
A third approach that can be adopted is that of continuum modeling which follows the number of cells N0(t) at a continuous time t. Whereas in the discrete and age-structured models the cell population renewal, differentiation and death are assumed to occur only at the end of cell cycles, in a continuous model these processes are assumed to occur continuously in time. In the discrete and age-structured models we consider the proportions ai, bi and c of the populations that undergo each process at the end of each cell cycle, but in the continuous model we must monitor the rates at which these processes occur in time, which we denote by their Greek equivalents αi, βi and γ (measured in hours−1), as shown in Figure 1. Thus, for example, over a small time τ we would expect a proportion α2 τ of stem cells to have differentiated, on average.
This continuous approach yields ordinary differential equations (ODEs) that follow the population as it changes in time. The corresponding ODE for the stem cell compartment is
which has a solution of the form
where α = α3 − α1 − α2 represents the net per-capita growth rate of stem cells.
Note that in order to compare formula (5) with (2) and (7) it is necessary to integrate the solution over all possible ages to obtain the total numbers for each cell population (as in ref. 5). These solutions are compared in Figure 2, with an initial age profile for (5) where all the cells have the same age. Each of the solutions has the same overall growth behavior, and any of these models could be used to describe the stem cell population. However, which approach is the most appropriate to use when modeling all the cells in the crypt, and which produces the model that is most mathematically tractable?
Let us consider both the stem cells and the transit-amplifying cells, which are assumed to divide on different timescales, t0 and t1 respectively. There is an implicit assumption of synchrony of cell division associated with the discrete model, where all the cells divide at the same time and are at the same stage in their cell cycle. Although not biologically realistic, in the case of only one cell type this assumption is at least self-consistent because every cell will divide once every generation. However, the discrete model cannot be used to describe two cell types with different generation times since, except in the special case where the stem cell cycle time is an integer multiple of the transit cell cycle time, t0 = pt1 (where p is an integer), the cells that are differentiating from the stem cell compartment will enter the transit cell population when the existing transit cells are already partly through their cell cycle, which breaks the synchrony condition.
Hence, in the case of two populations dividing at different rates, it is necessary to use either the age-structured or continuous models. If it is necessary to follow specific variations at times on the order of the cell cycle then the age-structured model is the best choice. The age-structured model is more accurate than the continuous model but produces more complicated solutions, as illustrated in ref. 5 where we solved the age-structured model for two different initial age profiles. Alternatively, when a population changes over a timescale of many cell cycles, the continuous model can capture the same overall behavior with more clarity: ODEs are much easier to solve mathematically than difference equations or PDEs. In ref. 5 we showed that the continuous and age-structured models produced solutions with the same behavior, and we derived conditions relating their parameters.
In conclusion, discrete modeling can be an informative technique when there is only one cell division timescale to consider, but if there is more than one timescale, as is the case with cells in the colonic crypt, then an age-structured or continuous approach must be adopted, and then the choice comes down to the timescale of interest. In Example 2 we investigate the global stability behavior of the system over long timescales, and so we use a continuous model.
Homeostasis in the crypt is observed to control cell numbers, but this regulation is broken in tumorigenesis where cell populations grow without bound. In this section we discuss the need to model homeostasis, how it can be overcome in tumorigenesis, and how a mathematical model can capture both these processes.
One feature of all of the solutions (2), (5) and (7) to the models in Example 1 is that the stability of the population is critically dependent on the choice of one particular parameter value. The age-structured and discrete model solutions vary in time according to functions of the form (2a3)G, where G is the number of divisions that the cells have undergone. The continuous model solutions vary with time-dependent functions eαt, where α represents the net per-capita growth rate of the stem cells. It is often assumed that stem cells normally divide asymmetrically, producing one daughter stem cell and one daughter transit cell each generation. However it is equally possible that stem cells, sometimes at least, divide symmetrically producing two daughter cells of the same type, either both stem cells or both transit cells. In that case, for the age-structured model, the stem cell population will remain stable if, and only if, exactly half of all symmetric divisions produce two stem cells while the other half produce two transit cells, a3 = 1/2. This is equivalent to the stem cells dividing asymmetrically on average. If a3 > 1/2 then the population will increase exponentially while if a3 < 1/2 then the population will decrease to zero. The same applies in the continuous model, with exponential growth (decay) occurring if α is greater (less) than zero, and a stable population occurs only if α = 0. When the stability of a system is dependent on a specific parameter value being taken a model is said to be structurally unstable. Such models are usually unrealistic since there is often some natural variation in parameter values.
Furthermore, when some stem cells are removed from the crypt, it has been observed that a homeostatic response is induced to create more stem cells.9 It is not known whether this is due to more stem cells dividing symmetrically, or to de-differentiation of first generation transit cells, but neither scenario can be modeled when these parameters are forced to take fixed values. A dynamic system requires regulation of cell numbers to maintain equilibrium.
One method that captures the homeostatic regulation of the stem cells is to introduce feedback into the system by making the parameter values dependent on the population sizes. In ref. 5 we used the continuous model to include these variable parameters by allowing the rate of stem cell differentiation to depend on the size of the stem cell population. Equally, we could have chosen the feedback to act through the rates of death or renewal instead.
Let us replace α2 by the variable differentiation rate R = R(N0), and note that the rates of stem cell death and renewal remain fixed throughout. If the stem cell number N0 becomes too large then R will increase so that more stem cells differentiate than renew. In practice this might be achieved by producing more transit cells and fewer stem cells, either by changing the proportion of symmetric divisions, or by changing the cell cycle time of the symmetric divisions. Alternatively, if the number of stem cells drops too low then R will decrease and symmetric divisions will produce more stem cells. In order to achieve both these effects, R(N0) must be an increasing function of N0. Many functions R(N0) could be chosen that would have this property, but the different functional forms could produce different qualitative effects on the system.
In our paper,5 we proposed feedback in the differentiation rate of the form
where R(N0) replaces α2 in (6), and k0 and m0 are nonnegative constants. The parameter k0 is measured in hours−1 and represents the speed of response of the feedback with large (small) k0 inducing a fast (slow) feedback response from the system. The parameter m0 is dimensionless and represents feedback saturation.
When m0 > 0, the differentiation rate cannot increase beyond a maximum size α2 + k0 / m0, and there is a nonzero stem-cell steady state of the form
as long as the net per-capita growth rate α = α3 - α1 - α2 is in the range 0 < α < k0 / m0. If the growth rate exceeds this value, then the stem cell population grows without bound. This is a saturating feedback because the population dependence in the differentiation rate will regulate the population size if the growth rate is low enough, but once the growth rate passes its saturation limit a steady state can no longer be achieved.
When m0 = 0, the feedback reduces to a linear form and the saturation is switched off. This means that any increase in the size of N0 will be matched by an equivalent increase in R, and there is a nonzero stable steady state
as long as the net per-capita growth rate α is positive, i.e., as long as cells are not removed faster than they are produced. This is the same form as the feedback chosen by Wodarz,8 who used a logistic form with the un-mutated rate of change of stem cells taken to be αN0(1 − N0/K) with a stem-cell carrying capacity K = α / k0. The linear feedback controls the population growth to the extent that unbounded growth is not permitted for any positive parameter values.
Which of these two forms is more appropriate for describing the cell population dynamics in the crypt? Both forms capture homeostatic regulation of stem cells and produce models that are structurally stable. However, in a system with the linear form, unbounded growth can only be achieved by removing the feedback completely (by setting the speed of feedback response k0 to zero), but with the saturating form unbounded growth can also be achieved by overwhelming the feedback. Therefore, only the saturating form can describe both homeostasis and unbounded growth of cell numbers without removing the feedback. There are two possibilities for how the saturating feedback might occur. The stem and transit populations might both be subject to saturating feedback, and consequently both populations would be able to initiate tumorigenesis. Alternatively, the stem cells may be subject to linear feedback which maintains equilibrium, with the transit cells subject to saturating feedback and possible tumorigenesis. We conclude this section by discussing both possibilities and giving examples of each.
Firstly, we consider the case of both the stem and transit cell populations being subject to saturating feedback. We choose the same form of feedback for the transit cells as for the stem cells, replacing the transit cell differentiation rate β2 by β2 + k1N1 / (1 + m1N1). Since the fully-differentiated cells do not divide and are only removed from the system at a rate γ, there is no need to include feedback in that population. The ODEs and cell population steady states for this system are shown in the Appendix. The stability of the replenishing cycle in the crypt is dependent on just two parameters: the net per-capita growth rates α = α3 - α1 - α2 and β = β3 - β1 - β2 for the stem and transit cell populations respectively. This stability is summarized in Figure 3 for a range of values of the parameters α and β. If α < 0 and β < k1 / m1 (Region I), the stem cells cannot sustain their number and consequently the crypt dies out. If 0 < α < k0 / m0 and β < k1 / m1 (Region II), the crypt will be in a normal equilibrium. If either α > k0 / m0 or β > k1 / m1 (or both) there is unregulated cell population growth in the system. For α > k0 / m0 and β < k1 / m1 (Region III), the cancer stem cells driving tumorigenesis are the tissue stem cells, whereas for α < k0 / m0 and β > k1 / m1 (Region IV), the cancer stem cells are derived from transit cells and the stem cell population either remains fixed (α > 0) or dies out (α < 0). Finally, if both α > k0 / m0 and β > k1 / m1 (Region V), unbounded growth is driven by cancer stem cells derived from both tissue stem cells and transit cells. So, in summary, the crypt will only maintain a healthy balance of proliferation, differentiation and death if the growth rates satisfy 0 < α < k0 / m0 and β < k1 / m1.
Secondly, we consider the case where the stem cells are now subject to linear feedback and the transit cells are again regulated by saturating feedback. The same form of feedback is chosen for the transit cells as in the last case, and the steady states for this system are shown in the Appendix. The stability of the population numbers is again dependent on the two net per-capita growth rates α = α3 - α1 - α2 and β = β3 - β1 - β2, which is summarized in Figure 4. If α < 0 and β < k1 / m1 (Region I), the stem cells cannot sustain their number and consequently the crypt dies out. If α > 0 and β < k1 / m1 (Region II), the crypt will be in a healthy equilibrium. If β > k1 / m1 (Region III), the tissue stem cells will either remain fixed (α > 0) or die out (α < 0), and unbounded growth is initiated from a cancer stem cell derived from the transit cell population. Here, unbounded growth can only be driven by transit cells and not tissue stem cells.
It is assumed that a mutation gives a cell a selective advantage,4 which would change the parameters associated with that cell and its feedback. This parameter change will increase the size of the steady states unless all the parameters change in combination such that the steady state is unaltered and the mutation is neutral. This increase in steady state could be interpreted as the first stage in the process transforming a normal crypt via an adenomatous polyp to a carcinoma.
It is observed that tumors grow to a certain size and then plateau until another growth spurt occurs. There are several long lag phases between periods of tumor growth, before the cancer finally sets in. Our feedback model can capture this behavior, with each mutation leading to a parameter change and consequently an increase in the steady state size. If any of these mutations result in α > k0 / m0 or β > k1 / m1, then the saturating feedback is overcome, there is no steady state, and exponential growth is initiated.
We illustrate this mutational process with two numerical examples. For this purpose, we make the simplifying assumption that once a mutation occurs its selective advantage is instantly conferred to all the other cells. Firstly, we consider the system with saturating feedback in both stem and transit cells, with the initial parameter set α1 = 0.1, α2 = 0.3, α3 = 0.69, β1 = 0.1, β2 = 0.3, β3 = 0.397, γ = 0.139, k0 = m0 = 0.1, k1 = 0.0003 and m1 = 0.0004, which gives α = 0.286, β = −0.0027 and produces critical saturation threshold values of α = k0 / m0 = 1 and β = k1 / m1 = 0.75. Therefore the population is stable with N0* = 4, N1* = 85 and N2* = 200 in Region II of Figure 3. If a mutation occurs (in either β1 or β3) that increases β to 0.1, then the stem cells remain at N0* = 4 and the other population steady states increase to N1* = 410 and N2* = 1,197. Now suppose a second mutation raises β to 0.2, which further raises the steady states to N1* = 925 and N2* = 3,344. Finally, a third mutation raises β to 0.8 which is larger than the critical threshold and lies in Region IV of Figure 3. Consequently there can no longer be a steady state, resulting in unregulated growth in the cell populations. Note that the initiating mutation for each of these changes may have occurred in stem cells but a selective advantage is only expressed in transit cells. The effect of this series of mutations is illustrated in Figure 5, where the different parameter regimes corresponding to each phase of growth are labeled (1), (2), (3) and (4) respectively, and these are also marked on the (α, β) parameter space in Figure 3.
For the second numerical example of a sequence of mutations, we consider the system with linear feedback in stem cells and saturating feedback in transit cells. The initial parameter set is α1 = 0.1, α2 = 0.3, α3 = 0.69, β1 = 0.1, β2 = 0.3, β3 = 0.388, γ = 0.1345, k0 = 0.07, k1 = 0.0002 and m1 = 0.0004, which gives α = 0.286, β = −0.0117 and produces the critical saturation threshold β = k1 / m1 = 0.75. Therefore the population is stable with N0* = 4, N1* = 85 and N2* = 200 in Region II of Figure 4. If a mutation occurs (in either β1 or β3) that increases β to 0.1, then the stem cells remain at N0* = 4 and the other population steady states increase to N1* = 654 and N2* = 1,962. Now suppose a second mutation raises α to 0.42, which further increases the steady states to N0* = 6, N1* = 676 and N2* = 2,042. Finally, a third mutation raises β to 1 which is larger than the critical threshold and lies in Region III of Figure 4. Consequently there can no longer be a steady state resulting in unregulated growth in the cell populations. These different regions of the parameter space are labeled (1), (2), (3) and (4) respectively in Figure 4.
Which parameters are the most influential in a system? This is an important question that mathematical modeling can help to address by using a sensitivity analysis. If it can be shown that some parameters have a strong influence while others have little or no effect, then experiments need only focus (at least initially) on the most influential processes. In this section, we discuss sensitivity coefficients, and how they can be used to ascertain key parameters.
By way of example, let us consider the effect that changes in the net per-capita growth rate α have on the stem cell steady state N0*. Suppose we change α by a small amount Δα, and assume that this induces a small change in N0* of ΔN0*. The sensitivity coefficient of N0* with respect to α is defined to be the relative change in N0* divided by the relative change in α
In the limit as Δα tends to zero, this can be expressed in terms of the derivative
If S is independent of α then we are effectively finding S such that N0* αS, but in general S will be a function of all the parameters in the system, so this relation is only true at the particular point in the parameter space at which S has been calculated. If S > 1, a percentage change in α will produce a larger percentage change in N0*; if S = 1, a percentage change in α causes equal percentage change in N0*; if 0 < S < 1, then a percentage change in α produces a smaller percentage change in N0*; if S = 0, a change in α has no effect on N0*; and if S < 0, then an increase (decrease) in α causes a decrease (increase) in N0*.
Applying formula (12) to the stem cell steady state with saturating feedback given in equation (9), the sensitivity coefficient of N0* with respect to α is SN0*,α = 1 + m0N0*. Since both N0* and m0 are positive, SN0*,α > 1 which means that an increase in α will produce a greater percentage change in N0*. Evaluating the sensitivity coefficient using parameter values α = 0.286 and k0 = m0 = 0.1 gives SN0*,α = 1.4. Alternatively, if α = 0.8, then SN0*,α = 5 and the dependence of N0* on α is much stronger. This illustrates how the dependence on α can vary at different points in the parameter space.
Now let us consider the sensitivity coefficient of the stem cell steady state to the parameter m0, denoted by SN0*,m0. It can be shown that SN0*,m0 = SN0*,α −1. Since the sensitivity coefficient SN0*,m0 will always be less than SN0*,α, we can say that the stem cell steady state is more sensitive to α than it is to m0.
By performing this analysis for each parameter to find the effect on each steady state, it can be shown that changes in the parameters α or k0 have the largest effect on the stem cell number, changes in β or k1 have the greatest effect on the semi-differentiated cell population, and changes in β, γ or k1 cause the biggest changes in the fully-differentiated cell population. Changes in m0 or m1 produce a less significant change in the steady states. However, we must remember that m0 and m1 are crucial in determining the stability boundaries in Figures Figures33 and and44.
The power of mathematical modeling is that it can provide insight into the behavior of complex interacting processes. In this Extra View we have outlined some of the relevant techniques that can be used to model the cell population dynamics in the colonic crypt, and we have illustrated these techniques with three examples.
Firstly, we discussed discrete, age-structured and continuous modeling approaches, why age-structured modeling is necessary to capture the dynamics of more than one population dividing on different timescales, and how continuum modeling is a more accessible tool.
Secondly, we explained how both homeostasis and tumorigenesis can be modeled by adopting a saturating form of feedback.5 We presented a scenario for the progression to cancer through a succession of mutations that lead to steady states with greater numbers of cells, until one mutation causes the net per-capita growth rate of stem or transit cells to exceed a critical value and unregulated cell population growth begins. In doing this, we made the simplifying assumption that the selective advantage gained by a mutation will be instantly conferred to the whole population, but a more detailed model might track the progress of healthy and mutant cells separately.
Finally, we discussed how a sensitivity analysis can be used to highlight the key parameters in a model. Modeling can prioritize the list of parameters that is to be measured. Results suggest that the key parameters are the net per-capita growth rates of stem and transit cells, the speed of response of stem and transit cells to changes in their number and the rate at which fully-differentiated cells are removed/die. We have shown how by incorporating the more stringent linear feedback into the stem cells, but not the transit cell compartment, it is then only the transit cells that can become the cancer stem cells. This may well be the best reflection of the actual situation in the large colorectal crypt.
As modeling and experiments become ever more intertwined, and increasing levels of system complexity are found, we believe that mathematical modeling will play a significant part in future research developments.
We acknowledge the support provided by the funders of the Integrative Biology project: the EPSRC (GR/S720231/01) and IBM. Matthew D. Johnston was supported by an EPSRC DTA graduate studentship (Award No. EP/P500397/1) and the Sarah and Nadine Pole Scholarship, which are gratefully acknowledged. Walter F. Bodmer was supported by a program grant from Cancer Research UK. Philip K. Maini was partially supported by a Royal Society-Wolfson Research Merit Award and NIH Grant U56CA113004 from the National Cancer Institute.
The system of ODEs describing the change of the cell populations in the crypt, with saturating feedback in the differentiation rates of stem and transit cells, is
which gives steady states
where D = α2N0* + k0N0*2/(1+m0N0*) = (α3 − α1)N0* is the stem-cell differentiation rate.
Linear feedback in the stem cells and saturating feedback in the transit cells. When the stem cells are subject to linear feedback and there is saturating feedback governing the transit cells, the system of ODEs (13)-(15) and the steady states (16)-(18) remain the same, except now the stem cell saturation is switched off (m0 = 0).