Previously published, basic models without spatial restriction

An ordinary differential equation model has been used to describe tissue hierarchy dynamics in a healthy tissue [8,15,24], and the models presented here build on these approaches. While cell lineages consist of stem cells, transit amplifying cells, and terminally differentiated cells, we can make a simplification and take into account only stem cells (which encompass all the proliferating cells) and differentiated cells. Hence, with this simplification, the stem cell compartment can be thought to contain both stem and transit amplifying cells, and such models have been shown to display qualitatively similar properties compared to three-compartment models with respect to the types of dynamics investigated here [6,24] (below we explore a more complex model that requires three compartments to be explicitly included). Denoting stem cells (SC) by S and differentiated cells (DC) by D, the model is given by: (1) Stem cells divide with a rate r. With a probability p, the division results in two daughter stem cells (self-renewal), and with a probability (1-p), the division results in two daughter differentiated cells (differentiating division). Differentiated cells are assumed to die with a rate α. These equations capture a probabilistic model of tissue control, where on the cell population level a fraction of the symmetric divisions result in two daughter stem cells and the remaining fraction results in two daughter differentiated cells. In addition to symmetric divisions, asymmetric divisions may play a role in tissue renewal. With asymmetric cell division, a stem cell gives rise to one stem cell and one differentiated cell, thus maintaining a constant population of stem cells. While not included here, it has been previously shown that the incorporation of asymmetric cell divisions into this modeling framework does not fundamentally alter the properties of the model [29], and in fact mathematically can be reduced to the same equation structure as (1)[30]. System (1) is only characterized by a neutrally stable family of nontrivial equilibria if p = 0.5 [6,8]. If p>0.5, infinite growth is observed. If p<0.5, the cell population goes extinct.

If we include the assumption that differentiated cells secrete negative feedback factors that influence stem cell division patterns, however, more realistic dynamics can be observed [6,8,24]. In particular, it has been shown that negative feedback from differentiated cells onto the probability of self-renewal, p, results in the existence of a stable equilibrium of cells, which might contribute to tissue homeostasis. This is because increased numbers of differentiated cells shift the division pattern in favor of differentiation, which limits overall cell growth. Mathematically, this has been expressed by p(D) = p’/(1+f₁D^κ1), where constant p’ is the basic self-renewal probability of stem cells in the absence of any feedback, and parameters f₁ and κ₁ determine the strength and power of the feedback. The existence of a stable equilibrium requires p’>0.5.

In addition, feedback loops have been proposed where differentiated cells reduce the rate of stem cell division [8,24], which can be expressed as r = r’/(1+f₂D^κ2), where r’ is the basic stem cell division rate in the absence of any feedback (and similarly, parameters f₂ and κ₂ determine the strength and power of the feedback). While this feedback can influence the dynamics of the system, it does not contribute to the existence of a stable equilibrium, and hence homeostasis.

Under these assumptions, the stable equilibrium population size is given by [29]

An agent-based model with diffusion of feedback mediators

We consider a two-dimensional spatial model of growing cell populations that are regulated by a population of diffusing negative feedback mediators. The model contains a cell layer and a layer that describes the dynamics of feedback factors.

The dynamics of the cell population are described by a stochastic agent-based model, which assumes two cell populations: stem cells and differentiated cells. The agent-based model is given by a 2-dimensonal grid that contains nxn spots. Each spot can be empty, contain a stem cell, or contain a differentiated cell. During one time step, the grid is sampled N times, where N is the total number of cells in the system. If the sampled spot contains a stem cell, a division event can occur with a probability p_div. A target spot is chosen randomly from the eight nearest neighboring spots, into which the offspring cell can be placed. If that spot is already occupied, the division event is aborted, introducing density dependence into the model (similar to logistic growth in ODEs). If the target spot is empty, a new stem cell is placed there with a probability p_self (self-renewal), and a differentiated cell is placed there with a probability 1- p_self. A stem cell can die with a probability p_Sdeath. If the sampled spot contains a differentiated cell, cell death occurs with a probability p_Ddeath. In addition to these processes, stem and differentiated cells can attempt a migration event with a probability p_mig. A spot is selected randomly from the eight nearest neighbors and if this spot is empty, the cell moves there. While the division events can in principle be formulated differently, the rules implemented here are meant to correspond to the previously published ODEs [6,8,24], with the addition of logistic growth rather than unlimited growth in the absence of feedback (see next section).

In contrast to previously published ODE models reviewed above [6,8,24], we here include the possibility of stem cell death for two reasons. First, tissue stem cells do die, although with a relatively small rate. In addition, when the model formulation contains a carrying capacity for the cell population, the absence of stem cell death leads to an artificial model outcome where the stem cells can exclude (or “outcompete”) the differentiated cells population, which is biologically unrealistic.

Note that in this algorithm, not all the cells will necessarily be sampled (and some might get sampled more than once). It is however necessary to tie a (biological) time-step with population size, to be able to compare different populations with each other. The following simple scaling argument can be made: in a simple birth process with a per-cell division rate L time^-1 and a population size N, during one time unit a cell will divide with probability L, so a total of NL cells will be expected to divide per time-unit. Similarly, in a spatial system, we need to sample N times per time-unit to achieve that average population division rate. In general, the larger the population size, the more events will happen during a time unit. Our approach is one of a number of ways to simulate this type of a stochastic process, another being the spatial Gillespie simulation, which implements a similar idea by setting the expectation of the time to the next reaction by the inverse of the sum of all the propensities (given by NL in the above example). Further details are provided in previous work [6].

The probability of self-renewal is influenced by the presence of feedback factors that can be secreted from differentiated cells. The dynamics of the feedback factors are described by a deterministic patch model, which includes production by differentiated cells, degradation, and diffusion of factors to the neighboring patches. Each spot on the cell grid is associated with a patch, in which the level of feedback factors is recorded. In each patch i, the level of the feedback factor, z_i, is given by the following ODE: dz_i/dt = c–bz_i−mgz_i + gZ_i. The parameter c represents the production rate of the feedback factor. We set c = 0 if the spot does not contain a differentiated cell, otherwise c>0. Feedback factors are assumed to decay with a rate b in each patch. Feedback factors move to the nearest neighboring patches with a rate g, representing diffusion processes. The parameter m denotes the number of neighboring patches. For our grid, m = 8, except for boundary patches, where m<8. The variable Z_i = Σ_jz_j denotes the sum of all feedback factor populations, z_j, among the patches adjacent to i. For each time step of the agent-based model, the ODEs in each patch were run for one time unit. The probability of self-renewal for a given cell in the agent-based part of the model is thus given by p_self = p⁽⁰⁾_self / (1+hz_i), where p⁽⁰⁾_self denotes the self-renewal probability in the absence of any feedback. The more feedback factors are locally present in the patch corresponding to the spot in the cellular grid, the lower the probability of stem cell self-renewal (and the larger the probability of differentiation). The parameter h describes the strength of feedback inhibition, with larger values of h corresponding to more potent inhibition. For simplicity, we do not include feedback on the stem cell division rate, because this has been shown in the ODEs to not contribute to the existence of a stable equilibrium. This feedback can, however, be easily incorporated into the model in the same way.

We note that the time scales of the agent-based portion of the model and the ODEs are the same, but that the processes occurring in these two parts of the model can occur at very different rates. As mentioned above, in the agent-based model, one time unit is defined by sampling N times, where N is the number of cells in the system. For each time unit in the agent-based model, the ODEs describing the dynamics of the feedback factor are also advanced by one time unit (as both cell dynamics and factor diffusion happen concurrently in the system). The amount of change that occurs during that one time unit can be very different for the agent-based portion and the ODEs, and is expressed in the rate constants and probabilities that determine the rate at which processes happen in the two parts of the model. We explore a large parameter space where time-scales of these different processes vary significantly.

The degree of spatial restriction in this system is given by two parameters: the migration rate of cells, p_mig, and the rate of feedback diffusion, g. If p_mig = 0 and the value of g is low, the system is characterized by strong spatial restriction. For large values of p_mig and g, the system approaches mass action dynamics (perfect mixing).

As initial conditions, a square of 7x7 spots in the center of the grid was filled with stem cells, and the simulation was run according to the rules described above. Parameter values are largely unknown, and hence were chosen for the purpose of demonstration. Parameters quoted in the figures are picked such that stem cells divide on average once a day [31,32], and terminally differentiated cells on average live for 10 days. The reported results, however, do not depend on these particular values, as was explored by repeated simulations and parameter sweeps. The basic codes that simulate the different models are provided in S1–S4 Data.

Dynamics assuming strong mixing of cells and feedback mediators

We start by assuming large values of p_mig and g, i.e. the migration rate of cells and the diffusion rate of feedback factors. In this limit, the average dynamics of the agent-based model converge to those predicted by the corresponding set of ordinary differential equations. Denoting the populations of stem cells, differentiated cells, and feedback factors by S, D, and Z, respectively, the ODEs are given as follows. (2) where negative feedback from differentiated cells onto the probability of self-renewal is given by , and p’ is the self-renewal probability of stem cells in the absence of feedback. In contrast to model (1) above, here we explicitly track the population of secreted feedback factors, Z. They are produced by differentiated cells with a rate ξ and decay with a rate β. Another change in the model (compared to model (1)) is the inclusion of a carrying capacity K, which corresponds to the maximum population size the system can sustain, independent of the negative feedback loop. This describes the finite grid size underlying the agent-based model. In the term describing the negative feedback on the stem cell self-renewal probability, the abundance of the secreted feedback mediator is divided by the carrying capacity K. The reason is that in the agent-based model, a larger system results in a lower amount of soluble feedback mediators available per cell due to diffusion, and this has to be captured in the corresponding ODE. Finally, we assumed that stem cells can die with a rate η, in addition to being removed through differentiation.

Fig 1A shows that this system of ODEs describes the average behavior of the agent-based model well for the limit of large cell migration and feedback diffusion rates. More generally, this ODE model is characterized by two equilibria. If p’<0.5, the cell population goes extinct and the system converges to the trivial equilibrium S⁽⁰⁾ = 0, D⁽⁰⁾ = 0, Z⁽⁰⁾ = 0. If p’>0.5, the system converges to a stable equilibrium at which all cell populations exist (not written down here due to complexity of expressions). An important property of this equilibrium is that the total cell population size is proportional to the carrying capacity, K (Fig 1B). As the value of the carrying capacity increases towards an infinitely large size, the number of cells also increase to an infinitely large size, despite the occurrence of negative feedback regulation. Therefore, in this model, negative feedback alone cannot maintain homeostasis of the cell population, defined as keeping the number of cells constant. The strength of negative feedback does, however, regulate the density of cells for a given carrying capacity, K. This behavior is in contrast to model (1) that does not take into account a carrying capacity. The reason is that model (1), as well as similarly structured models [24], assume that a given total amount of feedback mediators is equally effective regardless of the system’s size, while model (2) assumes that the same amount of feedback mediators becomes less effective at suppressing stem cell self-renewal for larger system sizes. In other words it assumes that locally, biological processes respond to concentrations rather than the total amounts of an factor, which corresponds to the formulation of the agent-based model and might be biologically more realistic.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Properties of the spatially explicit computational model with negative feedback.

(A) Dynamics of the agent-based model, depending on the carrying capacity K under the assumption that cells migrate with a high rate (large p_mig) and feedback mediators move across space at a high rate (large g). The horizontal lines represent equilibrium values derived from the corresponding ODE system (2). Parameters are given as follows. For agent-based model: P_div = 4.17x10^-2, p⁽⁰⁾_self = 0.7, P_Sdeath = 0, P_Ddeath = 8.3x10^-3, p_mig = 0.67, h = 1.6, c = 8.33, b = 4.17, g = 83.3; n = 100 and n = 200 for the small and large system, respectively. The average over 46 simulations are shown for each case; standard errors are too small to see. For ODEs: r = 4.17x10^-2, η = 0, α = 8.3x10^-3, ξ = 8.33, β = 4.17, p’ = 0.7, f = 1.6, K = 100x100 and 200x200 for the small and large systems, respectively. (B) Equilibrium properties of the corresponding ODE system (2) as a function of the carrying capacity, K. Parameters were the same as in (a). (C) Dynamics of the agent-based simulation, depending on the carrying capacity K under the assumption of spatial restriction (p_mig = 0, low g). Parameters were chosen as follows. P_div = 4.17x10^-2, p⁽⁰⁾_self = 0.7, P_Sdeath = 8.3x10^-5, P_Ddeath = 4.17x10^-3, p_mig = 0., h = 4.0x10^-3, c = 8.33, b = 8.3x10^-3, g = 0.83; n = 100 and n = 200 for the small and large system, respectively. The average over 46 simulations are shown for each case; standard errors are too small to see. Units of parameters are in hours.

https://doi.org/10.1371/journal.pcbi.1010039.g001

Dynamics assuming spatial restrictions

Here we assume that the processes of cell migration and the diffusion of feedback mediators are spatially restricted. For simplicity, we assume the strongest form of spatial restriction of cell migration, i.e. no migration occurs, and cells are assumed to place their offspring in a nearest neighboring spot. For feedback mediators, the rate of diffusion will be varied from low to high. In general, similar properties are observed compared to the well-mixed system studied in the previous section. That is, negative feedback can only regulate the density of the cells, but not total numbers: the larger the size of the grid, the larger the total cell population size (Fig 1C). Beyond the properties of the well-mixed model, however, we find that negative feedback inhibition can result in pronounced spatial patterns, in which a number of islands of stem and differentiated cells exist, separated by space that is not occupied by cells (Fig 2A). Within these islands, stem cells cluster in the inside, while differentiated cells are located around the periphery of the islands. Biological relevance of these clusters is discussed further below (Discussion Section). As the strength of feedback inhibition is reduced (lower value of parameter h), the distribution of cells within the space becomes uniform (Fig 2B). We quantified the degree of clustering by calculating the ratio of the variance to mean of the number of cells (which we denote by I_disp for index of dispersion) assuming the grid is subdivided into a number of relatively small squares. A ratio that is significantly greater than one indicates that the cell population is clustered across space. Fig 2C shows that as the strength of negative feedback inhibition is reduced, there is a remarkably sharp transition from a clustered distribution of cells towards a uniform distribution. This indicates the existence of two distinct states in the system, depending on whether the negative feedback can control density (leading to the patchy configuration) or whether this breaks down, resulting in a significantly higher density. If the clustered population structure corresponds to a normal state, then it is possible that the transition to a uniform distribution due to reduced negative feedback corresponds to a first step towards abnormal growth. Interestingly, in the presence of strong negative feedback and clustered population structure, the stem cell population is in the minority (Fig 2). For weak negative feedback and uniform cell distribution, however, the stem cell fractions are significantly larger (Fig 2), corresponding to stem cell enrichment [33].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Spatial patterns observed in the agent-based model with negative feedback.

Dark blue is empty space, light blue represents stem cells, and yellow represents differentiated cells. (A) With a stronger degree of negative feedback on the stem cell self-renewal rate, clumped spatial patterns are observed. Islands of stem and differentiated cells form, separated by empty space. Stem cells are in the minority. (B) For weaker negative feedback, these spatial patterns break down, and a uniform distribution of cells across space is observed. Also, stem cells become the dominant population. For A and B: The spatial picture is a snapshot in time at equilibrium, and the time series represents the average over 46 simulations; standard errors are too small to see. (C) The degree of clustering of cells across space can be quantified by dividing the space up into relatively small squares (10x10), and recording the number of cells per square. If the ratio of I_disp = variance / mean is greater than 1, the spatial pattern is clumped. If the ratio I_disp is less than one, the distribution is uniform. The graph shows the value of I_disp at the end of the simulation (at equilibrium). As the rate of negative feedback inhibition is increased from low to high, we observe a relatively sharp transition in the ratio I_disp, i.e. from a uniform to a clumped distribution of cells across the space. Baseline parameter values were chosen as follows. P_div = 4.17x10^-2, p⁽⁰⁾_self = 0.7, P_Sdeath = 8.3x10^-5, P_Ddeath = .17x10^-3, p_mig = 0., h = 4x10^-3, c = 8.33, b = 8.3x10^-3, g = 0.83; n = 200. For (A) h = 4x10^-3, (B) h = 10⁻³, and for (C) the value of h was varied, as indicated. Units of rate parameters are in hours^-1.

https://doi.org/10.1371/journal.pcbi.1010039.g002

Finally, we note that if the strength of negative feedback inhibition is large and crosses a threshold in this model, extinction of the cell population is observed.

Next, we investigated how the different outcomes depend on model parameters. This was done by randomly drawing the logarithm of parameter values from a uniform distribution and recording the ratio I_disp. In Fig 3, two parameters were varied simultaneously for any given plot: The rate of negative feedback inhibition on stem cell self-renewal, h, was always varied, together with a second parameter. The outcome of the simulation is color-coded in Fig 3, where each dot in the parameter space represents the outcome of an individual simulation. A value of I_disp >1.5 (clustering) is recorded in red, and a value of I_disp <1.5 (uniform distribution) is recorded in blue. Runs in which population extinction occurred are shown in yellow. As mentioned above, clustered cell population persistence occurs for intermediate rates of negative feedback inhibition, while larger and smaller rates of feedback inhibition result in population extinction, or a uniform distribution of cells across space, respectively. The width of the parameter region in which cell clusters are observed depends on parameters (Fig 3). In particular, it depends on the diffusion rate of feedback mediators, and on the decay rate of the feedback mediators. Larger diffusion rates, g, and slower decay rates of feedback mediators, b, lead to a broader region of feedback inhibition rates for which cell clusters are observed (Fig 3). If the diffusion rate of feedback mediators is low and their decay rate is high, the mediators secreted by a given cell mostly act locally, and the clustering regime is narrow. In the limit, the clustering regime is so narrow that the behavior of the system essentially transitions from extinction to a uniform invasion of the space by cells (Fig 3). If the diffusion rate of feedback mediators is faster and they decay slower, then the mediators secreted from a given cell affect cells in a larger area of the space, and the clustering behavior becomes pronounced (Fig 3). Therefore, the model suggests that the clustered persistence of cells requires feedback mediators to act beyond the immediate neighborhood of the cell from which they are secreted. Other model parameters do not have a significant influence on the range of feedback inhibition values across which clustered cell persistence is observed (Fig 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Parameter dependencies of outcome.

Each dot in the graph represents the long-term outcome of an individual simulation. Each simulation was run up to a time threshold, and the spatial distribution was determined by calculating I_disp, the index of dispersion. The time threshlold was determined by waiting until the temporal average of the stem cells did not change by more than 0.1% for 1000 consecutive time steps, after which the simulation was allowed to run for an amount of time that corresponds to 10 years in the simulation, to ensure that the dynamics are well in the steady state pase. Yellow indicates the extinction of the cells. Blue indicates a distribution that is characterized by I_disp < 1.5 (mostly uniform distribution). Red indicates I_disp > 1.5 (clumped distribution of cells). For each graph, two parameters were varied: the strength of negative feedback on the stem cell self-renewal probability, h, and a second parameter, as indicated in the individual graphs. The ranges of the parameters were chosen such that all three outcomes are seen, to illustrate how the different outcomes depend on parameters. Baseline parameters were chosen as follows. P_div = 4.17x10^-2, p⁽⁰⁾_self = 0.7, P_Sdeath = 8.3x10^-4, P_Ddeath = 4.17x10^-3, p_mig = 0., c = 8.33, b = 8.3x10^-3, g = 0.833; n = 200. Units of rate parameters are in hours^-1.

https://doi.org/10.1371/journal.pcbi.1010039.g003

Feedforward loop

So far, we have considered regulatory loops where differentiated cells influence the behavior of stem cells. Data from the airway epithelial tissue from mice [34], however, indicate that stem cells can also send a signal forward to their progeny and influence their behavior. In fact, this “feedforward” regulation is inherently spatial. The stem cells secrete a notch ligand to their daughter transit amplifying cells (secretory cells), and this signal is necessary to maintain the transit amplifying cell population. In the absence of this signal, the transit amplifying cells undergo terminal differentiation to become ciliated cells. Hence, the further the transit amplifying cells are located away from stem cells, the weaker this signal, and the more likely terminal differentiation occurs. Transit amplifying cells that are located closer to stem cells receive a stronger signal and are more likely to self-renew.

We adapt our agent-based model to describe this scenario. To do so, we expand the model complexity to include a population of transit amplifying cells (TA) in addition to stem cells. Now, stem cell differentiation results in the generation of two TA cells. Similarly to stem cells, TA cells have a probability to divide (q_div). With a probability q_self, this is a self-renewing division, giving rise to two TA daughter cells; with a probability 1-q_self, this is a terminally differentiating division, giving rise to two daughter differentiated cells. TA cells are assumed to die with a probability q_Tdeath. Feedforward mediators are secreted from stem cells. The dynamics of the feedforward factors are again described by a deterministic patch model. In each patch i (associated with a spot on the grid), the level of the feedforward factor, w_i, is given by the following ODE: dw_i/dt = c₂ –b₂w_i−mg₂w_i + g₂W_i. The parameter c₂ represents the production rate of the feedback factor. We set c₂ = 0 if the spot does not contain a stem cell, otherwise c₂>0. Feedforward factors are assumed to decay with a rate b₂ in each patch. Feedforward factors move to the nearest neighboring patches with a rate g₂, representing diffusion processes. The parameter m again denotes the number of neighboring patches. The variable W_i = Σ_jw_j denotes the sum of all feedforward factor populations, w_j, among the patches adjacent to i. Importantly, the probability of TA cell self-renewal (and thus maintenance) is determined by the amount of the feedforward factor according to , consistent with the experimental observations [34]. For large feedforward factor levels, the probability of TA self-renewal converges to one. In the absence of feedforward factors, the probability of TA self-renewal is zero. In other words, the feedforward factors secreted by stem cells are responsible for the maintenance of the TA cell population, and in the absence of these mediators, terminal differentiation is the only fate of the division.

It is instructive to start examining the dynamics of this system under the assumptions that stem cells and TA cells do not die (and only disappear through differentiation). In this case, a structure forms where stem cells are located in the center, surrounded by an area of TA cells, with differentiated cells at the periphery (Fig 4Ai). Interestingly, this structure is self-contained and stops growing after an initial phase of expansion. The size of this structure is independent of the grid size, and determined by the diffusion rate of the feedforward factor. The further the feedforward factors can diffuse from the stem cells, the larger the area of the TA cell population. In other words, in this scenario, the feedforward regulation can maintain homeostasis of cell numbers, and the same cell numbers are maintained no matter how large the available space is. Hence, homeostasis is maintained by the control loops and not by the carrying capacity of the system. The reason for this behavior is that there is competition between stem and TA cells for space, similar to the dynamics described by Hillen et al [35]. Stem cells enable self-renewal of TA cells in their immediate vicinity, and the TA cells block the stem cells from dividing further. Hence, the growth of the overall cell population is limited.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4.

(A) Properties of the agent-based models with the feedforward control loop only. Dark blue represents empty space, light blue stem cells, green transit amplifying cells, and yellow terminally differentiated cells. (i) In the absence of stem and transit amplifying cell death, the cell population stops growing and converges to an equilibrium that is independent of the carrying capacity, K (not shown). The transit amplifying cells block stem cell divisions due to lack of available space, and this prevents the cluster of cells from expanding. The picture shown corresponds to the population at steady state. Parameters were as follows: P_div = 4.17x10^-2, p⁽⁰⁾_self = 0.8, q_div = 5.83x10^-2, P_Sdeath = 0, P_Tdeath = 0, P_Ddeath = 4.17x10^-3, p_mig = 0., c₂ = 0.833, b₂ = 8.3x10^-3, g₂ = 0.417, h₂ = 2, n = 200. (ii) In the presence of cell death, however, this mechanism breaks down and the stem cells can continuously expand into empty space, provided by the death of transit amplifying cells. The picture represents a snap-shot during this cell expansion. Parameters were chosen as follows: P_div = 4.17x10^-2, p⁽⁰⁾_self = 0.8, q_div = 5.83x10^-2, P_Sdeath = 8.3x10^-5, P_Tdeath = 10⁻⁴, P_Ddeath = 4.17x10^-3, p_mig = 0., c₂ = 0.833, b₂ = 8.3x10^-3, g₂ = 0.417, h₂ = 2, n = 200. (B) Properties of the agent-based model that contains both the feedforward and the feedback loop, and assumes the occurrence of death for all cell populations. (i) Time series, in which the system / grid size was varied, n = 100 vs n = 150. Blue and purple show stem cells for the smaller and larger grid size, respectively. Light green and dark green show TA cells, for the smaller and larger grid size, respectively. Yellow and orange show differentiated cells, for the smaller and larger grid size, respectively. The lines present the average time series over 46 iterations of the simulation, and the dashed lines represent the average plus minus standard errors. The simulation was run for 100 years to show that this mechanism can in principle maintain tissue homeostasis for long human life-spans (in the absence of any mutations that might allow cell growth) (ii) A snapshot of the spatial configuration of cells at a specific time point during steady state. Parameters were chosen as follows. P_div = 4.17x10^-2, p⁽⁰⁾_self = 0.8, q_div = 5.83x10^-2, P_Sdeath = 8.3x10^-5, P_Tdeath = 10⁻⁴, P_Ddeath = 4.17x10^-3, p_mig = 0., c = 8.33, b = 8.33x10^-2, g = 0, h = 6.0x10^-2, c₂ = 8.33, b₂ = 8.33x10^-2, g₂ = 3.33, h₂ = 2.5. Units of rate parameters are in hours^-1.

https://doi.org/10.1371/journal.pcbi.1010039.g004

This homeostasis, however, represents a somewhat artificial situation due to the assumption that stem and TA cells can only be eliminated through differentiation processes, and no explicit cell death is assumed to occur. If cell death is assumed to occur in stem and TA cells, the spatial competition dynamics either go in favor of the TA cells, resulting in the exclusion of stem cells and thus in the extinction of the whole cell population, or in favor of stem cells, in which case the stem cell population can expand outward over time, leading to cell population growth limited only by the available space (Fig 4Aii), similar to the simulations with negative feedback on stem cell self-renewal. Hence, the control-loop-mediated homeostasis observed in this system is lost in the presence of stem cell death.

Combination of negative feedback and feedforward regulation

The above sections considered negative feedback on the stem cell self-renewal probability, and feedforward regulation separately. For each model, these mechanisms failed to maintain control loop-mediated homeostasis of cell numbers. Here, we consider a modified version of the agent-based model that includes both of these control processes at the same time. The two control mechanisms are implemented in the same way as before. We now observe parameter regions in which the control loops maintain a constant number of cells independent of the carrying capacity of the system. In other words, it is the control loops that maintain homeostasis and not the carrying capacity (Fig 4B), i.e. the number of cells converges to an equilibrium value that is independent of the size of the grid or carrying capacity (including parameter regions where stem and transit amplifying cells are assumed to die). The reason for this behavior is as follows. With the feedforward loop, the death of transit amplifying cells can make space for the invasion of self-renewing stem cells. If only the feedforward loop exists, this results in a continuous expansion of the stem cell population through space, which results in a concomitant expansion of transit amplifying and differentiated cells, and hence in uncontrolled growth. The addition of negative feedback on stem cell self-renewal, however, counters this process. The negative feedback mediators are assumed to be secreted from differentiated cells, which are located predominantly at the surface of the cell mass (Fig 4B). As the stem cells expand in space and reach locations closer to the differentiated cells, the negative feedback loop becomes important and counters this outward stem cell expansion by forcing terminal differentiation to occur. This prevents the unbounded expansion of the stem cells and ensures the existence of an equilibrium that is independent of the grid size (carrying capacity). This effect is observed even if the negative feedback factors act predominantly on a local level, including with zero diffusion rates (as long as the feedback mediators remain present for long enough to stop an expanding front of stem cells).

Exploration of the parameter space is computationally not feasible on a larger scale because each parameter combination needs to be run many times to obtain the average trajectories and thus to determine whether grid size determines the equilibrium number of cells. Additional simulations, however, are described in S1 Text (ii), demonstrating that these dynamics are observed over large parameter ranges.

It is important to point out that spatial dynamics are essential to observe homeostasis that is maintained by the control loops and independent of the carrying capacity of the system. It is possible to describe a non-spatial version of the feedforward control mechanism, where stem cells promote the self-renewal of TA cells regardless of spatial location. A corresponding ODE model that takes into account both the feedforward and the negative feedback mechanisms is presented in S1 Text (i). In the absence of space, the ability of these controls to maintain homeostasis is lost, and the equilibrium cell population sizes are directly proportional to the carrying capacity of the system.