Figures
Abstract
Feedback mechanisms within cell lineages are thought to be important for maintaining tissue homeostasis. Mathematical models that assume well-mixed cell populations, together with experimental data, have suggested that negative feedback from differentiated cells on the stem cell self-renewal probability can maintain a stable equilibrium and hence homeostasis. Cell lineage dynamics, however, are characterized by spatial structure, which can lead to different properties. Here, we investigate these dynamics using spatially explicit computational models, including cell division, differentiation, death, and migration / diffusion processes. According to these models, the negative feedback loop on stem cell self-renewal fails to maintain homeostasis, both under the assumption of strong spatial restrictions and fast migration / diffusion. Although homeostasis cannot be maintained, this feedback can regulate cell density and promote the formation of spatial structures in the model. Tissue homeostasis, however, can be achieved if spatially restricted negative feedback on self-renewal is combined with an experimentally documented spatial feedforward loop, in which stem cells regulate the fate of transit amplifying cells. This indicates that the dynamics of feedback regulation in tissue cell lineages are more complex than previously thought, and that combinations of spatially explicit control mechanisms are likely instrumental.
Author summary
Tissues in higher organisms are maintained at a homeostatic setpoint level, which is determined by a complex set of regulatory mechanisms. Negative feedback loops from differentiated cells onto stem cell division patterns, which have been experimentally documented, are thought to contribute to homeostatic regulation. This is in part based on the analysis of mathematical models. Typically, however, such models are based on ordinary differential equations, which assume perfect mixing of cells and molecule, while in vivo, tissues are characterized by pronounced spatial structure. Here, we show that this negative feedback mechanism loses the ability to maintain tissue homeostasis in spatially explicit computational models. Negative feedback on stem cell division, however, can lead to the organization of stem and differentiated cells into distinct, spatially separated islands, where stem cells comprise a small fraction of the cell population. Reduction of negative feedback control in the model results in a breakdown of these spatial structures and stem cell enrichment. We further show that control loops can maintain homeostasis in this spatial setting if the negative feedback on stem cell divisions is combined with other, experimentally documented, spatially explicit control mechanisms, indicating that more complex combinations of control loops are instrumental.
Citation: Uhl P, Lowengrub J, Komarova N, Wodarz D (2022) Spatial dynamics of feedback and feedforward regulation in cell lineages. PLoS Comput Biol 18(5): e1010039. https://doi.org/10.1371/journal.pcbi.1010039
Editor: Feng Fu, Dartmouth College, UNITED STATES
Received: August 20, 2021; Accepted: March 18, 2022; Published: May 6, 2022
Copyright: © 2022 Uhl et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript and its Supporting Information files.
Funding: The work was funded by grant NSF DMS-1815406 (DW, NK), grant NIH/NCI U54-CA217378 (JL, DW, NK), and NSF-Simons Center for Multiscale Cell Fate Research, CMCF NSF DMS-1763272 (NK). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Summer Salary was obtained from grant DMS-1815406 (DW and NK), DMS-1763272 (NK), U54-CA217378 (DW and JL).
Competing interests: The authors have declared that no competing interests exist.
Introduction
Tissue homeostasis is central to the functioning and survival of higher organisms. Adult tissues are maintained by tissue stem cells that can both self-renew and differentiate to give rise to transit amplifying cells, which in turn give rise to terminally differentiated cells. Stem cells can divide asymmetrically [1], where one daughter is another stem cell and the other is a differentiating cell. In mammalian systems, however, data indicate that a stochastic symmetric division model also plays an important part, where a stem cell gives rise to either two daughter stem cells or to two daughter differentiated cells [2–4]. In such settings, it is thought that feedback mechanisms are required to prevent unbounded cell growth and to maintain homeostasis. Corresponding feedback signaling molecules have been much discussed in the literature in different tissues, and loss of feedback signals has been implicated in carcinogenesis [5,6]. Several feedback molecules have been shown to determine the fate of cell divisions, influencing whether self-renewal or differentiation can occur. Examples include growth differentiating factor 11 (GDF11) and Activin βB (ACTβB), which negatively regulate self-renewal rates in progenitor and stem cells in the olfactory epithelium of mice [7,8]; transforming growth factor beta (TGF-β) [9], which is mutated in a variety of tumors [10–12]; the bone morphogenetic protein 4 pathway (BMP4) that is inactivated in glioblastomas [13]; and the adenomatous polyposis coli (APC) tumor suppressor gene that is inactivated in colorectal cancer, with concomitant activation of the Wnt cascade [14].
A growing mathematical literature has emerged that investigates feedback control in relation to tissue homeostasis, and loss of feedback control in relation to carcinogenesis [6–8,15–26]. One particular approach focused on the notion that feedback factors produced by differentiated cells might play an important role for determining the fate of cell divisions. Specifically, in the olfactory epithelium, there appears to be negative feedback from differentiated cells both on the probability of stem cell self-renewal and on the rate of cell division, mediated by GDF11 and Activin β [7,8]. These observations motivated mathematical models showing that such negative feedback from differentiated cells onto stem cell division patterns can play an important role in the maintenance of tissue homeostasis [8,21–24,27]. Mathematical models predicted that in the absence of this feedback, unbounded growth occurs, while negative feedback from differentiated cells inhibiting stem cell self-renewal can result in a stable equilibrium.
These mathematical models of negative feedback regulation are typically based on ordinary differential equations (ODEs), which assume perfect mixing (mass action or meanfield assumption) of cells and molecules. In other words, no spatial structure was assumed. This applies to most mathematical models of stem cell regulation, with a few exceptions, e.g. [28]. Here, we re-formulate negative feedback models within cell lineages assuming spatially restricted dynamics, and specifically examine how this affects the ability of negative feedback loops to maintain tissue homoeostasis. We further construct models of other control mechanisms where spatial interactions are important, and determine their effect on tissue homeostasis. We start by briefly reviewing the previously published ODE models that were used to study negative feedback from differentiated cells on stem cell division patterns in relation to the ability to maintain tissue homeostasis. We then introduce the spatial modeling framework and explore how this affects the dynamics.
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+f1Dκ1), where constant p’ is the basic self-renewal probability of stem cells in the absence of any feedback, and parameters f1 and κ1 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+f2Dκ2), where r’ is the basic stem cell division rate in the absence of any feedback (and similarly, parameters f2 and κ2 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]
Results
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 pdiv. 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 pself (self-renewal), and a differentiated cell is placed there with a probability 1- pself. A stem cell can die with a probability pSdeath. If the sampled spot contains a differentiated cell, cell death occurs with a probability pDdeath. In addition to these processes, stem and differentiated cells can attempt a migration event with a probability pmig. 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, zi, is given by the following ODE: dzi/dt = c–bzi−mgzi + gZi. 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 Zi = Σjzj denotes the sum of all feedback factor populations, zj, 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 pself = p(0)self / (1+hzi), where p(0)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, pmig, and the rate of feedback diffusion, g. If pmig = 0 and the value of g is low, the system is characterized by strong spatial restriction. For large values of pmig 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 pmig 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) = 0, D(0) = 0, Z(0) = 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.
(A) Dynamics of the agent-based model, depending on the carrying capacity K under the assumption that cells migrate with a high rate (large pmig) 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: Pdiv = 4.17x10-2, p(0)self = 0.7, PSdeath = 0, PDdeath = 8.3x10-3, pmig = 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 (pmig = 0, low g). Parameters were chosen as follows. Pdiv = 4.17x10-2, p(0)self = 0.7, PSdeath = 8.3x10-5, PDdeath = 4.17x10-3, pmig = 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.
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 Idisp 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].
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 Idisp = variance / mean is greater than 1, the spatial pattern is clumped. If the ratio Idisp is less than one, the distribution is uniform. The graph shows the value of Idisp 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 Idisp, i.e. from a uniform to a clumped distribution of cells across the space. Baseline parameter values were chosen as follows. Pdiv = 4.17x10-2, p(0)self = 0.7, PSdeath = 8.3x10-5, PDdeath = .17x10-3, pmig = 0., h = 4x10-3, c = 8.33, b = 8.3x10-3, g = 0.83; n = 200. For (A) h = 4x10-3, (B) h = 10−3, and for (C) the value of h was varied, as indicated. Units of rate parameters are in hours-1.
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 Idisp. 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 Idisp >1.5 (clustering) is recorded in red, and a value of Idisp <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).
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 Idisp, 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 Idisp < 1.5 (mostly uniform distribution). Red indicates Idisp > 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. Pdiv = 4.17x10-2, p(0)self = 0.7, PSdeath = 8.3x10-4, PDdeath = 4.17x10-3, pmig = 0., c = 8.33, b = 8.3x10-3, g = 0.833; n = 200. Units of rate parameters are in hours-1.
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 (qdiv). With a probability qself, this is a self-renewing division, giving rise to two TA daughter cells; with a probability 1-qself, this is a terminally differentiating division, giving rise to two daughter differentiated cells. TA cells are assumed to die with a probability qTdeath. 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, wi, is given by the following ODE: dwi/dt = c2 –b2wi−mg2wi + g2Wi. The parameter c2 represents the production rate of the feedback factor. We set c2 = 0 if the spot does not contain a stem cell, otherwise c2>0. Feedforward factors are assumed to decay with a rate b2 in each patch. Feedforward factors move to the nearest neighboring patches with a rate g2, representing diffusion processes. The parameter m again denotes the number of neighboring patches. The variable Wi = Σjwj denotes the sum of all feedforward factor populations, wj, 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.
(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: Pdiv = 4.17x10-2, p(0)self = 0.8, qdiv = 5.83x10-2, PSdeath = 0, PTdeath = 0, PDdeath = 4.17x10-3, pmig = 0., c2 = 0.833, b2 = 8.3x10-3, g2 = 0.417, h2 = 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: Pdiv = 4.17x10-2, p(0)self = 0.8, qdiv = 5.83x10-2, PSdeath = 8.3x10-5, PTdeath = 10−4, PDdeath = 4.17x10-3, pmig = 0., c2 = 0.833, b2 = 8.3x10-3, g2 = 0.417, h2 = 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. Pdiv = 4.17x10-2, p(0)self = 0.8, qdiv = 5.83x10-2, PSdeath = 8.3x10-5, PTdeath = 10−4, PDdeath = 4.17x10-3, pmig = 0., c = 8.33, b = 8.33x10-2, g = 0, h = 6.0x10-2, c2 = 8.33, b2 = 8.33x10-2, g2 = 3.33, h2 = 2.5. Units of rate parameters are in hours-1.
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.
Discussion
Previous mathematical modeling approaches [6,8,24,27], based on the assumption that cells mix perfectly (mass action), suggested that negative feedback from differentiated cells on the self-renewal probability of stem cells is an important determinant of tissue homeostasis. Here, we extended this analysis into a spatially explicit scenario where cells grow in a finite space, characterized by a given carrying capacity (maximum possible number of cells the system can sustain). This suggests that the negative feedback from differentiated cells onto the stem cell self-renewal probability cannot by itself maintain tissue homeostasis. The number of cells at equilibrium always scales with the carrying capacity of the system, even under the assumption that cells and feedback mediators mix relatively well due to fast migration and diffusion processes. In these models, an infinitely large space available for cell growth will lead to infinitely large cell population sizes, despite the presence of the negative feedback loop.
The negative feedback on stem cell self-renewal can only maintain cell density, not total cell numbers in the models considered here. In the spatially explicit model versions (without significant migration of cells), this negative feedback can further lead to the formation of spatial patterns. In the presence of relatively strong negative feedback, clumps of cells (containing both stem and differentiated cells) form in the model, separated by empty space. This is interesting, because it suggests that this negative feedback loop might be involved in the formation of spatial structures in the cellular microenvironment. At the same time, however, the formation of spatial structures around stem cell populations, in particular the formation of stem cell niches, is highly complex and includes many (partially unknown) components not currently taken into account in our model [36]. Therefore, it is difficult to interpret details of these spatial structures, such as the location of stem cells in the center and the location of differentiated cells around the periphery in Fig 2A. While this could be quantified further in the model, the physiological relevance remains unclear because these predicted spatial arrangements can be easily altered when further biological complexities, characteristic of stem cell niches, are introduced into the model. This was therefore not explored further.
As the strength of feedback inhibition is reduced, there is a sharp transition in the model away from the clumped spatial structure towards a uniform distribution of cells, where the stem cell population is dominant. This means that loss of negative feedback on stem cell self-renewal might lead to a collapse of spatial organization, which might be a first step towards malignancy, even though unbound growth does not yet occur. Whether this is indeed the case requires further investigation, in particular using models that take into account a higher biological complexity that characterize stem cell niche morphology.
In this context, it is interesting to note that in the regime with relatively strong feedback, where pronounced spatial structures form, stem cells are in the minority and form a small fraction of the total cell population (Fig 2A). As feedback inhibition becomes weaker and spatial structures disappear, the fraction of stem cells markedly increases in the model, i.e. we observe stem cell enrichment (Fig 2B). Stem cell enrichment is a phenomenon that occurs during carcinogenesis [33], supporting our notion that the transition from spatial structures to a more uniform distribution of cells could correlate with disease development.
While the model suggests that negative feedback on stem cell self-renewal can lead to spatial structures, the model further indicates that stronger spatial restrictions of feedback mediator diffusion limits the parameter regime in which this behavior is observed. The pattern formation is observed over relatively wide parameter regions only if the rate of feedback mediator diffusion is fast, i.e. if feedback factors secreted by a given cell can affect cells in a location further away. If feedback mediators are assumed to act only locally (through limited diffusion), the parameter region where the spatial patterns form becomes extremely narrow. In this case, as the strength of the negative feedback is reduced from high to low, the model behavior more or less transitions from population extinction (stronger feedback) to a uniform distribution of cells (weaker feedback).
It was interesting to observe, however, that a combination of the negative feedback loop with a feedforward loop from stem cells to transit amplifying cells can lead to homeostasis that is maintained by the control loops and independent of the carrying capacity (e.g. a self-limiting solution); that is, cell numbers settle around an equilibrium that is independent of the amount of space available. The feedforward loop assumed in the model was based on data from the airway epithelium in mice, where stem cells were shown to secrete a notch ligand, which enabled the maintenance of the transit amplifying cells [34]. In the absence of this feedforward signal, transit amplifying cells were shown to undergo terminal differentiation. This is an inherently spatial process, since the concentration of the feedforward mediator decreases with distance from the originating stem cell. This leads to the formation of a structure, where stem cells are located in the center, surrounded by transit amplifying cells, while differentiated cells are located at the surface of this area. The presence of only the feedforward control loop cannot maintain homeostasis because the stem cell population can expand outward and replace transit amplifying cells. The concomitant presence of the negative feedback loop, however, results in increased differentiation (rather than self-renewal) of stem cells, limiting their ability to spread outwards in space. This contributes to a stable equilibrium that is independent of the carrying capacity of the system. This type of equilibrium is observed even if the diffusion rate of the negative feedback mediator is assumed to be slow, i.e. if the negative feedback loop acts largely on a local level. Importantly, an equivalent non-spatial model does not exhibit control-mediated cell homeostasis. These results suggest that the interplay between different feedback control mechanisms within a cell lineage in a spatial setting can lead to homeostasis that is maintained only by the control loops (independent of carrying capacity), which we could not observe in a model without spatial structure.
While the combination of these particular feedback and feedforward mechanisms can maintain control loop-mediated homeostasis in our model (independent of carrying capacity), it is likely that other feedback configurations exist that give rise to similar results. In addition, it is important to point out that this is a relatively simple model that was aimed to test how previously described homeostatic mechanisms [6,8,24] translate into spatial settings. Nevertheless, this model is based on solid experimental data that document the biological importance of the feedback and feedforward loops underlying the mathematical models. As mentioned above, however, additional biological complexities, especially those that characterize stem cell niches, need to be considered to gain a more detailed picture of cell dynamics under homeostatic conditions. Besides the control loops involved in cells of a given lineage, complex signaling mechanisms between the cell lineage and its microenvironmental components exist that together yield homeostatic properties. These processes will need to be decoded with a combination of experiments and dynamical models. The results described here form a basis for understanding feedback dynamics in such more complex and realistic settings.
Materials and methods
The core computational models underlying this study are given by stochastic agent-based models that describe the dynamics of cell populations, coupled with a deterministic patch model that describes the dynamics of soluble feedback mediators. The models are described in detail in the Results section. They were coded in Intel Fortran, and simulated on a computer cluster. The basic codes that simulate the different models are provided in S1–S4 Data. For the limit of perfect mixing, the dynamics can be given by ordinary differential equations, which are also described in detail in the Results section.
Supporting information
S1 Text. Additional mathematical and computational details about the modeling analysis.
https://doi.org/10.1371/journal.pcbi.1010039.s001
(PDF)
S1 Data. Computer code that generates data that were used for Fig 2.
Details about the code and how to run it are given within the file.
https://doi.org/10.1371/journal.pcbi.1010039.s002
(F90)
S2 Data. Computer code that generates data that were used for Fig 2.
Details about the code and how to run it are given within the file.
https://doi.org/10.1371/journal.pcbi.1010039.s003
(F90)
S3 Data. Computer code that generates data that were used for Fig 3.
Details about the code and how to run it are given within the file.
https://doi.org/10.1371/journal.pcbi.1010039.s004
(F90)
S4 Data. Computer code that generates data that were used for Fig 4.
Details about the code and how to run it are given within the file.
https://doi.org/10.1371/journal.pcbi.1010039.s005
(F90)
References
- 1. Inaba M, Yamashita YM. Asymmetric stem cell division: precision for robustness. Cell stem cell. 2012;11(4):461–9. pmid:23040475
- 2. Simons BD, Clevers H. Strategies for homeostatic stem cell self-renewal in adult tissues. Cell. 2011;145(6):851–62. pmid:21663791.
- 3. Snippert HJ, van der Flier LG, Sato T, van Es JH, van den Born M, Kroon-Veenboer C, et al. Intestinal crypt homeostasis results from neutral competition between symmetrically dividing Lgr5 stem cells. Cell. 2010;143(1):134–44. pmid:20887898.
- 4. Nicolas P, Kim KM, Shibata D, Tavare S. The stem cell population of the human colon crypt: analysis via methylation patterns. PLoS computational biology. 2007;3(3):e28. pmid:17335343; PubMed Central PMCID: PMC1808490.
- 5. Lee J, Son MJ, Woolard K, Donin NM, Li A, Cheng CH, et al. Epigenetic-mediated dysfunction of the bone morphogenetic protein pathway inhibits differentiation of glioblastoma-initiating cells. Cancer cell. 2008;13(1):69–80. pmid:18167341
- 6. Rodriguez-Brenes IA, Komarova NL, Wodarz D. Evolutionary dynamics of feedback escape and the development of stem-cell-driven cancers. Proceedings of the National Academy of Sciences of the United States of America. 2011;108(47):18983–8. Epub 2011/11/16. 1107621108 [pii] pmid:22084071; PubMed Central PMCID: PMC3223454.
- 7. Gokoffski KK, Wu HH, Beites CL, Kim J, Kim EJ, Matzuk MM, et al. Activin and GDF11 collaborate in feedback control of neuroepithelial stem cell proliferation and fate. Development. 2011;138(19):4131–42. pmid:21852401; PubMed Central PMCID: PMC3171217.
- 8. Lander AD, Gokoffski KK, Wan FY, Nie Q, Calof AL. Cell lineages and the logic of proliferative control. PLoS biology. 2009;7(1):e15. Epub 2009/01/27. [pii]. pmid:19166268; PubMed Central PMCID: PMC2628408.
- 9. Neuzillet C, Tijeras-Raballand A, Cohen R, Cros J, Faivre S, Raymond E, et al. Targeting the TGFbeta pathway for cancer therapy. Pharmacology & therapeutics. 2015;147:22–31. pmid:25444759.
- 10. Jakowlew SB. Transforming growth factor-beta in cancer and metastasis. Cancer metastasis reviews. 2006;25(3):435–57. pmid:16951986.
- 11. Tian M, Neil JR, Schiemann WP. Transforming growth factor-beta and the hallmarks of cancer. Cellular signalling. 2011;23(6):951–62. pmid:20940046; PubMed Central PMCID: PMC3076078.
- 12. Drabsch Y, ten Dijke P. TGF-beta signalling and its role in cancer progression and metastasis. Cancer metastasis reviews. 2012;31(3–4):553–68. pmid:22714591.
- 13. Lim DA, Tramontin AD, Trevejo JM, Herrera DG, García-Verdugo JM, Alvarez-Buylla A. Noggin antagonizes BMP signaling to create a niche for adult neurogenesis. Neuron. 2000;28(3):713–26. pmid:11163261
- 14. Krausova M, Korinek V. Wnt signaling in adult intestinal stem cells and cancer. Cellular signalling. 2014;26(3):570–9. pmid:24308963
- 15. Kunche S, Yan H, Calof AL, Lowengrub JS, Lander AD. Feedback, Lineages and Self-Organizing Morphogenesis. PLoS computational biology. 2016;12(3):e1004814. pmid:26989903; PubMed Central PMCID: PMC4798729.
- 16. Youssefpour H, Li X, Lander AD, Lowengrub JS. Multispecies model of cell lineages and feedback control in solid tumors. Journal of theoretical biology. 2012;304:39–59. pmid:22554945; PubMed Central PMCID: PMC3436435.
- 17. Rodriguez-Brenes IA, Kurtova AV, Lin C, Lee YC, Xiao J, Mims M, et al. Cellular Hierarchy as a Determinant of Tumor Sensitivity to Chemotherapy. Cancer research. 2017;77(9):2231–41. pmid:28235762; PubMed Central PMCID: PMC5487257.
- 18. Stiehl T, Ho AD, Marciniak-Czochra A. Mathematical modeling of the impact of cytokine response of acute myeloid leukemia cells on patient prognosis. Sci Rep. 2018;8(1):2809. pmid:29434256; PubMed Central PMCID: PMC5809606.
- 19. Arino O, Kimmel M. STABILITY ANALYSIS OF MODELS OF CELL PRODUCTION SYSTEMS Mathematical Modeling. 1986;7:1269–300.
- 20. Marciniak-Czochra A, Stiehl T, Ho AD, Jager W, Wagner W. Modeling of asymmetric cell division in hematopoietic stem cells—regulation of self-renewal is essential for efficient repopulation. Stem Cells Dev. 2009;18(3):377–85. pmid:18752377.
- 21. Komarova NL. Principles of regulation of self-renewing cell lineages. PloS one. 2013;8(9):e72847. pmid:24019882; PubMed Central PMCID: PMC3760876.
- 22. Komarova NL, van den Driessche P. Stability of Control Networks in Autonomous Homeostatic Regulation of Stem Cell Lineages. Bulletin of mathematical biology. 2017. pmid:28508298.
- 23. Yang J, Axelrod DE, Komarova NL. Determining the control networks regulating stem cell lineages in colonic crypts. Journal of theoretical biology. 2017;429:190–203. pmid:28669884; PubMed Central PMCID: PMC5689466.
- 24. Lo WC, Chou CS, Gokoffski KK, Wan FY, Lander AD, Calof AL, et al. Feedback regulation in multistage cell lineages. Math Biosci Eng. 2009;6(1):59–82. Epub 2009/03/19. pmid:19292508; PubMed Central PMCID: PMC2756546.
- 25. Nakata Y, Getto P, Marciniak-Czochra A, Alarcón T. Stability analysis of multi-compartment models for cell production systems. Journal of biological dynamics. 2012;6(sup1):2–18. pmid:22873671
- 26. Stiehl T, Marciniak-Czochra A. Characterization of stem cells using mathematical models of multistage cell lineages. Mathematical and Computer Modelling. 2011; pmid:21076663
- 27. Konstorum A, Hillen T, Lowengrub J. Feedback Regulation in a Cancer Stem Cell Model can Cause an Allee Effect. Bulletin of mathematical biology. 2016;78(4):754–85. pmid:27113934.
- 28. Chou C-S, Lo W-C, Gokoffski KK, Zhang Y-T, Wan FY, Lander AD, et al. Spatial dynamics of multistage cell lineages in tissue stratification. Biophysical journal. 2010;99(10):3145–54. pmid:21081061
- 29. Rodriguez-Brenes IA, Wodarz D, Komarova NL. Characterizing inhibited tumor growth in stem-cell-driven non-spatial cancers. Mathematical biosciences. 2015;270(Pt A):135–41. pmid:26344137.
- 30. Weiss LD, van den Driessche P, Lowengrub JS, Wodarz D, Komarova NL. Effect of feedback regulation on stem cell fractions in tissues and tumors: Understanding chemoresistance in cancer. Journal of theoretical biology. 2021;509:110499. pmid:33130064
- 31. Simons BD, Clevers H. Stem cell self-renewal in intestinal crypt. Experimental cell research. 2011;317(19):2719–24. pmid:21787769
- 32. Ritsma L, Ellenbroek SI, Zomer A, Snippert HJ, de Sauvage FJ, Simons BD, et al. Intestinal crypt homeostasis revealed at single-stem-cell level by in vivo live imaging. Nature. 2014;507(7492):362–5. pmid:24531760
- 33. Enderling H. Cancer stem cells: small subpopulation or evolving fraction? Integrative Biology. 2014;7(1):14–23. pmid:24840956
- 34. Pardo-Saganta A, Tata PR, Law BM, Saez B, Chow RD-W, Prabhu M, et al. Parent stem cells can serve as niches for their daughter cells. Nature. 2015;523(7562):597–601. pmid:26147083
- 35. Hillen T, Enderling H, Hahnfeldt P. The tumor growth paradox and immune system-mediated selection for cancer stem cells. Bulletin of mathematical biology. 2013;75(1):161–84. pmid:23196354.
- 36. Jones DL, Wagers AJ. No place like home: anatomy and function of the stem cell niche. Nature reviews Molecular cell biology. 2008;9(1):11–21. pmid:18097443