article: Feedback regulation of cell lineage progression plays an important role in tissue size homeostasis, but whether such feedback also plays an important role in tissue morphogenesis has yet to be explored. Here we use mathematical modeling to show that a particular feedback architecture in which both positive and negative diffusible signals act on stem and/or progenitor cells leads to the appearance of bistable or bi-modal growth behaviors, ultrasensitivity to external growth cues, local growth-driven budding, self-sustaining elongation, and the triggering of self-organization in the form of lamellar fingers. Such behaviors arise not through regulation of cell cycle speeds, but through the control of stem or progenitor self-renewal. Even though the spatial patterns that arise in this setting are the result of interactions between diffusible factors with antagonistic effects, morphogenesis is not the consequence of Turing-type instabilities.
  Diffusible molecules, such as morphogens, growth factors, cytokines and chemokines, play many important roles in solid tissues, among which are conveying positional information, coordinating proliferation with differentiation, and enabling collective responses to disturbances [1–6]. Less clear, however, is the importance of diffusible signals in orchestrating the morphogenesis of solid tissues—i. e. the formation of three-dimensional structures such as sacs, tubes, sheets, buds and branches. Typically, the events of morphogenesis are described in terms of intrinsic cellular behaviors, such as cell migration [7], localized changes in cell shape [8–12], differential proliferation [13,14], regulated cleavage plane orientation [15–18], differential cell adhesion [19], and mechanical force production [20]. How much of a role long-range diffusible cues play in organizing these processes is poorly understood. Recent mathematical modeling has suggested that the ability of diffusible molecules to spontaneously create large-scale spatial patterns (“Turing patterns” [21,22]) might drive branching morphogenesis in the lung [23]. Turing patterns have also been implicated in digit morphogenesis in vertebrate limbs [24]. As the necessary conditions for Turing patterning are relatively restrictive—for example, pairs of factors that differ greatly in their diffusivity are required—how commonly this mechanism is exploited in morphogenesis remains to be seen. Here we identify a different mechanism by which diffusible signals can drive morphogenesis—not one based on creating patterns that act upstream of tissue growth or cell movement, but on regulating cell lineage progression, so that patterns arise directly out of tissue growth and cell displacement. As we describe below, such interactions can create large-scale spatial structure, and generate sharply-demarcated tissue forms, both spontaneously and in response to graded and/or transient extrinsic cues. Underlying such behavior is a negative feedback loop that has primarily been studied for its role in size homeostasis of continually self-renewing tissues [2,25–27]. The basis of this mechanism is the production of diffusible molecules by differentiated cells that act upon stem or progenitor cells to decrease the probability of self-renewal, and increase the probability of differentiation. In muscle, for example, such feedback comes from GDF8/myostatin [28,29], whereas in the olfactory epithelium the related molecules GDF11 and activin play an analogous role [30,31]. Such feedback causes stem cell populations to automatically coordinate their behaviors so that net renewal is maintained at exactly 50%—the value required for homeostasis [25,32]. Feedback of this type also produces steady states that are “perfectly” robust, insofar as tissue size is maintained at a set-point independent of initial conditions, cell cycle speed, or rates of cell turnover [25]. During development, the regulation of tissue and organ growth is less a matter of achieving homeostasis, than of efficiently arriving at a correct final size and shape (i. e. , morphogenesis). Indeed, for the many organs in which cell turnover is non-existent or negligible in the adult, a controlled growth trajectory during development is the only possible function of feedback regulation. Fortunately, mathematical modeling shows that negative feedback control of self-renewal, the same strategy that achieves robust steady states in tissues that turn over, can also achieve nearly perfectly robust final states for tissues that do not turn over [25]. Indeed, feedback molecules such as GDF8, GDF11, myostatin and other renewal-inhibiting factors are often highly expressed during development. Yet, interestingly, such molecules are often co-expressed with other diffusible molecules that have exactly the opposite effect on stem and progenitors cells—promoting renewal at the expense of differentiation. In the developing olfactory epithelium, for example, the same cells that produce GDF11, which inhibits the renewal of progenitor cells, also produce FGF ligands that stimulate renewal of the same cells [33]. In muscle, production of renewal-inhibiting GDF8 is similarly accompanied by production of renewal-stimulating FGFs [34]. In other tissues, pairs of co-expressed renewal-promoting and -inhibiting factors include Wnts and BMPs [35]; HGF and TGFβ [36]; and PTHrP and FGF [37]. Mixing positive and negative feedback might at first seem counterproductive, but it is well known that countervailing feedback loops—if they are non-linear and operate on sufficiently well-separated temporal or spatial scales—can give rise to interesting dynamic or spatially-dynamic behaviors. Indeed, Turing patterns may be seen as a special case of this general phenomenon. Here we show that the simultaneous positive and negative feedback regulation of cell renewal by diffusible signals tends naturally to lead to the spontaneous emergence of spatially dynamic behaviors, which manifest as sharply localized differences in net growth. We mathematically define the conditions required for underlying phenomena such as bistability, ultrasensitivity and spatial self-organization to arise in such systems, and show how they can be used to drive the formation of induced or spontaneous bulges, buds and fingers (the basic building blocks of many types of morphogenesis). We further show that even though such phenomena arise out of the opposing interactions of diffusible factors, the patterns that arise are not Turing patterns. Rather, pattern formation is driven by the coupling between tissue growth and the inherent ultrasensitivity of the underlying system, and not by intrinsic differences in the spatial ranges of action of diffusible positive and negative signals.
  Cell lineages may be described as pathways of differentiation that begin with a stem cell (SC), progress through some number of self-renewing progenitor stages, and end with one or more terminally differentiated (TD) cells. Here we consider only the simplest of lineages—ones that are un-branched uni-directional, and produce TD cells that are non-dividing (postmitotic) —because, as we will see, the spontaneous emergence of morphogenetic behaviors requires no more complexity than that. Following [25], the dynamics of such lineages may be written as a system of ordinary differential equations (ODEs), shown in section 1 of S1 Text, in which cell number (or concentration) at any lineage stage i is represented by the variable χi, cell cycle speed at stage i by rate parameter vi, and the probability of self-renewal at stage i by parameter pi (0 ≤ pi ≤ 1). In tissues undergoing cell turnover, a rate constant, d, for cell loss or death must also be introduced (minimally for the TD cell). Within this framework, feedback or feedforward effects may be represented as dependencies of parameters upon the χi. When seeking to capture the effects of diffusible molecules, modeling with ODEs is equivalent to assuming that cells are homogenously distributed (“well-stirred”) in space, so that diffusion gradients may be ignored. Although this is not generally a valid assumption for solid tissues, we may still understand the ODE representation of lineage dynamics as applying to sufficiently small spatial domains within tissues. In addition, because solutions to deterministic ODEs do not reveal the effects of stochastic variation that one might expect to encounter at low cell numbers, their analysis must often be supplemented with independent stochastic simulations, . In spite of these limitations, ODE formulations can serve as a good starting point for developing general insights, which can then be extended to regimes in which stochasticity, space and spatial organization matter explicitly. This is the approach we take here. Fig 1A depicts the simplest possible, continually-renewing unidirectional lineage: a single type of SC that gives rise directly to a single type of TD cell that then turns over with a constant half-life. The ODEs for this situation may be written as dχ0 (τ) dτ= (2P (τ) −1) χ0 (τ) dχ1 (τ) dτ=2 (1−P (τ) ) χ0 (τ) −δχ1 (τ) (1) in which χ0 and χ1 represent SC (cell type “0”) and TD (cell type “1”) cells, respectively; P stands for the SC renewal probability; τ refers to time in units of cell cycles; and δ is the death rate of TD cells relative to the cell cycle rate, i. e. δ = d/v. To represent the possibility that TD cells may be a source of both negative and positive feedback on P, we replace P with the expression p (11+γ χ1 (τ) ) (φ χ1 (τ) 1+φ χ1 (τ) ), where γ and ϕ quantify the strength of the negative and positive feedback, respectively, while p, a constant between 0 and 1, represents the maximum possible renewal probability. Expressing feedback as a product of Hill functions is convenient insofar as it guarantees a non-negative value of P, and exhibits saturation at large values of χ1, as would be expected for any biological process. However, the precise form of the feedback function has only minor effects on the qualitative results described below. Likewise, the interaction between negative and positive feedback terms can be represented in a variety of other ways without altering the qualitative behaviors that are observed (see Tables A and B in section 5 of S1 Text, which summarize outcomes of different forms of the feedback functions). For non-zero values of ϕ, system 1, with P is defined as above, has three equilibrium states. In these states χ1 is equal to either 0, −γ−φ+2pφ−−4γφ+ (γ+φ−2pφ) 22γφ, or−γ−φ+2pφ+−4γφ+ (γ+φ−2pφ) 22γφ, and χ0 is simply δ times the value of χ1 in each state, which is a consequence of the non-zero χ0 equilibrium solutions necessitating that P = 0. 5. For trivial reasons the solution in which both χ0 and χ1 are zero (i. e. both cell types become extinct) is necessarily stable. For the non-zero solutions, a saddle-node bifurcation occurs so that the second solution is unstable, while the third solution locally stabilizes as shown in Fig 1B–1E. The criterion for bistability—in which either of the two stable states is available for a single set of parameters—is the condition that p>0. 5 (stem cells must be capable of renewing more than 50% of the time) and: φγ>1+22p+2p1−4p+4p2. Using the maximum feasible value of 1 for p, we find the minimum of this ratio to be ≈5. 8. In other words, bistability requires that positive feedback on progenitor self-renewal be at least 5. 8 times stronger than negative feedback. Fig 1F plots the trajectories of solutions with different ϕ/γ ratios, starting from an initial condition of zero TD cells. In every case the system starts to generate TD cells, but when ϕ/γ is not large enough, TD production collapses and both cell types vanish. The switch between extinction and maintaining a non-zero steady state occurs, as predicted, at ϕ/γ ≈5. 8, and this holds regardless of the individual values of ϕ and γ (see S1 Fig). The fact that extinction occurs when positive feedback is too weak may seem surprising, given that spontaneous extinction was not observed in prior modeling studies that considered only negative feedback (e. g. [25]), but this is simply because such studies had to proceed from the assumption that SC are intrinsically specified to renew with high probability when there is no feedback. Here, we abandon that assumption, allowing the renewal and differentiation behaviors of SC to be largely determined by extrinsic factors, as is likely to be the case in vivo. Although the existence of a non-zero steady state depends strictly on ϕ/γ, whether that state is reached is a function of initial conditions, as shown in Fig 1G. This is to be expected for a bistable system. If SC numbers are initially too low, only extinction is possible. If initial SC numbers exceed a threshold, the non-zero steady state can be reached only if the initial number of TD cells is neither too high nor too low. This behavior can be rationalized: if TD numbers are too low, there is insufficient positive feedback to keep SC numbers from contracting to the point where SCs are incapable of sustaining themselves; if TD numbers are too high, the effect of positive feedback approaches saturation and negative feedback dominates, driving SC self-renewal below the threshold for sustainability. To explore the effects of stochastic fluctuations, we allow p, the maximum possible self-renewal fraction to be a random variable, while maintaining the same form of P given above. The reason for focusing on the self-renewal parameter is that the greatest source of stochasticity is likely to come from randomness in cell-division outcomes, whereby cells choose probabilistically between whether to generate zero, one or two differentiated progeny. In Fig 1H–1J, we take p = 0. 9+ζ, where ζ a random variable sampled every CP cell cycle from a normal distribution with mean zero and variance 0. 2 (p is truncated to ensure it lies between 0 and 1). The results, in terms of the distributions of the TD cell numbers (χ1) after 100 CP cell cycles, measured over a large number of independent trials, are presented in Fig 1H–1J for different values of ϕ/γ. In these cases, the initial conditions were also similarly distributed random variables centered around the values from Fig 1F. The results show that, when ϕ/γ is small, the most likely state is the zero steady-state (Fig 1H). When ϕ/γ is sufficiently large, the non-zero steady state is the most likely, although the distribution of populations is fairly broad (Fig 1J). For intermediate values, the distribution is bimodal with the zero and non-zero steady states being nearly equally likely (Fig 1I) due to noise-driven switching between the steady states. Because of the possibility of stochastic extinction of small populations, the values of ϕ/γ needed to achieve non-zero steady states in the stochastic case are larger than those required for the deterministic model. Nevertheless, the behaviors of the stochastic and deterministic models are qualitatively similar. The system of ODE in (1) models lineages in which TD cells constantly turn over, as occurs in hematopoietic, epidermal and many epithelial lineages. In such a system, steady state solutions can be calculated simply by setting the time derivatives in (1) to zero. Not all tissues reach steady states, however. In some, TD cells last for the lifetime of the organism (e. g. in the nervous system), with SC either disappearing, or becoming quiescent. Such tissues may be modeled by setting the turnover rate δ in (1) to zero, as illustrated in Fig 2A (see also [25]). In this case, the governing equations become dχ0 (τ) dτ= (2P (τ) −1) χ0 (τ) dχ1 (τ) dτ=2 (1−P (τ) ) χ0 (τ) (2) again allowing feedback to be captured by a dependence of P on the levels of cell types (e. g. , we again take P = p (11+γ χ1 (τ) ) (φ χ1 (τ) 1+φ χ1 (τ) ) ). The lack of TD cell turnover in (2) means that TD cells are now free to accumulate to a level that, through negative feedback, drives SC to extinction. A stationary state is thereby reached, but it is not “steady” in the usual sense of the word, i. e. there is no dynamic equilibrium; we therefore refer to it as a “final state”. Note that TD values at the final state cannot be calculated simply by setting the time derivatives in (2) to zero, because the time derivatives and SC value (χ0) both go to zero together, leaving P undetermined. Strictly speaking, final-state systems require TD cells that never turn over. However, if TD cells turn over slowly, the final state formulation will provide a better description of lineage dynamics on time scales that are short relative to the TD cell half-life. As tissue morphogenesis often occurs rapidly, relative to TD cell lifespans, final state models may thus be inherently better representations of lineage dynamics during morphogenesis, even in self-renewing tissues. We therefore explored the effects of mixed feedback in final-state systems. As with steady-state systems, we observed behaviors of final-state systems that depended strongly on the value of the ratio ϕ/γ, and exhibited abrupt switching between distinct states in response to small changes in the initial conditions for either χ0 or χ1 (see plots in Fig 2). As described in Box 1, the dynamics of TD production can be understood in relation to two critical TD population sizes: χ1critlow and χ1crithigh. When χ1 lies between these two values, χ0 grows; otherwise, χ0 shrinks to zero. Understanding the dynamics that drive TD cell numbers across these critical values allows us to understand why initial conditions, seen in Fig 2D, cause leaps in growth to coincide with neither critical values nor equilibrium solutions. There are essentially three regimes of behavior: Mechanistically, the abrupt transitions between regimes, as a function of initial conditions, are analogous to the abrupt transitions between steady states in Fig 1, although the precise relationships between initial conditions and final state are not identical. Final-state systems may not properly be called bistable, since they are not attracted to fixed, dynamic stable states. We therefore refer to them as bimodal, in the sense that they may be switched between slow and fast growth modes as function of the sizes of cell populations. Like bistability, bimodality creates an opportunity for ultrasensitive switching. Not only can small differences in initial conditions cause abrupt changes in the growth trajectories of feedback-controlled tissues, as shown in Fig 2D, so can small changes in the strength of either positive or negative feedback (see S3 Fig). Similarities and differences between the possible behaviors of steady state and final state systems, analyzed over all possible parameter values, and in the face of different kinds of TD-derived feedback, are summarized in Table 1, and Tables A and B in section 5 of S1 Text (see also S4 Fig in which positive feedback arises from SCs). These results show that the presence of positive feedback is required for bimodality, the presence of negative feedback is required for stability, and to achieve either bistability, or bimodality with a stable final state, both positive and negative feedback are required. Next, we consider a three-stage, steady-state lineage, in which we also alter the source and target of positive feedback. In this particular scenario (Fig 3A), an intermediate, committed progenitor (CP) is the target of all feedback, and a constant influx of progenitor cells from a stem cell source guarantees non-extinction of the CP pool (we assume that the SCs are in equilibrium). Negative feedback from TD cells is now mixed with positive feedback from CP cells themselves. The ODEs for this system are dχ0 (τ) dτ=0dχ1 (τ) dτ=ζ χ0 (τ) + (2P (τ) −1) χ1 (τ) dχ2 (τ) dτ=2 (1−P (τ) ) χ1 (τ) −δχ2 (τ) (3) in which ζ is the rate of proliferation of SCs, relative to that of CPs, and P = p (11+γ χ2 (τ) ) (ϕ χ1 (τ) 1+ϕ χ1 (τ) ). Such an arrangement models, for example, the apparent roles of GDF11 and FGFs in the mouse olfactory epithelium [25]. As observed in the analysis of the two-stage lineage, we see here that the ratio ϕ/γ, which captures the relative strength of positive versus negative feedback, plays a determining role in the emergence of bistability and ultrasensitivity (Fig 3B), although now both stable states are non-zero (Fig 3B and 3C). Further, when p and the initial data are taken to be random variables of the type used in Fig 1H–1J, the distributions of the TD cell population at a given time are similarly dependent on ϕ/γ (Fig 3D–3F). When ϕ/γ is small, the system is subject to stochastic extinction events even though the stable steady states are non-zero (Fig 3D). When ϕ/γ is sufficiently large, the most likely state of the system is the largest stable population, although the distribution of the populations is again fairly broad (Fig 3F). When ϕ/γ is an intermediate value, the distribution tends to be bimodal with peaks near the two steady states; the probability of stochastic extinctions decreases with increasing ϕ/γ. Although the three-stage lineage model is not as analytically tractable as the two-stage model, numerical solutions clearly demonstrate that the system exhibits bistability with respect to feedback parameters and initial conditions (see S5A–S5F Fig). Similarly, growth bimodality is also observed if cell turnover is set to zero, or made very slow (see S5G and S5H Fig). Real tissues are not “well-stirred”, but grow outward into physical space. This potentially influences the feedback regulation of cell lineages in three ways. First, diffusible molecules that mediate feedback will naturally form spatial gradients, so that the concentrations of feedback factors—and thus the effects of feedback—will differ at different locations. Second, growth may displace (advect) both cells and diffusible molecules to potentially different degrees at different locations. Third, boundary conditions (i. e. the disposition of molecules and cells when they reach natural tissue boundaries) may influence both cell behaviors and molecular gradients. To account for spatial variation, we re-interpret the variables in Eq (2) as representing a local measure of cell concentration, i. e. the volume fraction occupied by each cell type, and add additional equations for the spatial distributions of negative and positive feedback factors, which we refer to as G and F, respectively. Because the time-evolution of the distributions of these molecules is expected to be much faster than time scales of growth, we assume a quasi-steady state distribution for them. We allow G to be produced by TD cells; F to be produced by CP cells, TD cells or both; and in some cases also consider the potential effects of external sources of G or F (e. g. , stromal cells). The distribution of F is governed by Eq (7) in Methods; the equation for G is analogous. The forms of the external sources and sinks are given in Eqs (13) – (17) in Methods. The effects of F and G on CP self-renewal are modeled exactly as in Figs 1–3, with parameters γ and ϕ representing the strengths of negative and positive feedback, respectively (see Eq (6) in Methods). For simplicity, growth and diffusion are simulated here in two-dimensional domains, (extension to three dimensions is straightforward). As an initial condition we typically begin with tissues represented as thin, rectangular domains, as might be encountered in growing epithelia; see the schematic in Fig 3G. We denote one long side as apical [AP], and impose there a reflective (no-flux) boundary condition on molecular diffusion (representative of the tight junctions that exist at the apical surface of most epithelia). On the other side, where one would normally find a basement membrane [BM] that is not itself a barrier to diffusion, we impose no boundary condition, allowing molecules to move unimpeded into the adjacent space, which may be taken to represent the stroma upon which most epithelia sit. In some cases, we also model the presence of molecules in the stroma that act as absorbers, or sinks, for diffusible factors (e. g. follistatin in the stroma beneath GDF11- and activin-expressing epithelia [25,38]), see Eqs (15) and (16) in Methods. Diffusion, proliferation, differentiation and cell turnover are then modeled together with soft-tissue mechanics (as described in [39,40]). The dynamics of the epithelial cells are governed by Eqs (4), (5) and (8) – (12) in Methods. In addition, the SC population is restricted to the basement membrane (BM) and set to self-renew with probability 1/2. Temporal and spatial trajectories of tissue growth in such models were explored over a wide range of parameter values to identify the sorts of qualitative behaviors that emerge. A complete list of dimensionless parameters can be found in S1 Table in Supplementary Materials. The nondimensionalization is given in section 6 of S1 Text. See Methods for a discussion of the numerical implementation of the governing Eqs (4) – (12) and section 6 of S1 Text for additional model and numerical details. As in non-spatial models, spatial systems may be initiated from small size and grow to include large numbers of CP or TD cells. Interestingly, in spatial models the presence of negative feedback control no longer guarantees that a stable steady state, or final-state, will be reached for all parameter values. If the strength of negative feedback is not sufficient, tissue size will increase indefinitely (see Fig 3H, blue curve). This reflects the fact that diffusible signals have limited ranges of action, so that the incremental increase in feedback, averaged over a tissue, that accompanies any incremental increase in tissue size, gets progressively smaller the larger a tissue gets. This effect notwithstanding, we clearly observe that the mixture of positive and negative feedback can produce growth that is bimodal, switching ultrasensitively between fast and slow (or zero) growth rates as a function of the ϕ/γ ratio. For example, Fig 3H–3J compare the epithelia dynamics using positive feedback strengths ϕ = 3. 0 and ϕ = 4. 0 (holding other parameters fixed; see S1 Table). Note that in this example, initial growth is similar for both feedback strengths, but in the ϕ = 3. 0 case (green curves and insets), the lack of sufficient positive feedback causes increased differentiation of CP cells into TD cells, elevating the concentration of G, which causes CP numbers to decline even further (Fig 3I). The decline in CPs eventually causes TD production to fall to a low level that is insufficient even to balance TD turnover (Fig 3J). In contrast, when ϕ = 4. 0, the exact opposite sequence of events occurs, and CP levels rise to a value that sustains continual growth (see blue curves and insets in Fig 3H–3J). The distributions of the feedback regulators, the self-renewal fractions and the cells can be found in S6 Fig for ϕ = 3. 0 and S7 Fig for ϕ = 4. 0 in Supplementary Materials. Similar results are found when CP cells also produce G (S8 Fig), TD cells also produce F (S9 Fig), the ratio DF/DG is decreased (S10 Fig), and the sink of G in the stroma is removed (S11 Fig). Although such systems exhibit a sharp threshold in behavior at a critical value of ϕ/γ, the value itself is not easily calculated, as it depends not only on initial conditions, but also on the spatial gradients on G and F, and thereby on the tissue geometry, as well as on the relative locations of cell types (see S12A Fig for a plot of the feedback ratios for the simulations in Fig 3H–3J). We also explored the effects of stochasticity on the results in this spatial model. In particular, we considered two sources of random fluctuations—the self-renewal parameter p and the relative proliferation rate ζ of SCs. In S12B Fig, both are taken to be normally-distributed random variables in space and time with means given by their deterministic values. The results indicate that, when the variance is relatively small, simulations evolve similarly to what is seen in the deterministic case. When the variance is large, self-sustaining growth may occur spontaneously even when it would not, under deterministic conditions. When ϕ is large enough (e. g. ϕ = 4. 0), we found that growth may self-sustain even for small variances. One consequence of the dependence of growth trajectories on the spatial distributions of feedback molecules is that localized, exogenous sources of the same molecules can be used to create localized, persistent differences in growth rate. In Fig 4 we consider an epithelial tissue in which the ϕ/γ ratio would not normally be sufficient to sustain continued growth, but an exogenous point source of F is placed near the basement membrane (such a source might represent the production of a growth factor by a group of stromal cells, for example). In response, a local outgrowth appears and continues to grow, producing a stratified arrangement of cells, with CP cell densities highest near the F-source and TD cells highest farther away. Notice that the outgrowth becomes sharply demarcated, i. e. there is a distinct notch, or “neck” that separates it from the adjacent epithelium (see arrows in Fig 4B), despite the fact that the exogenous gradient of F declines smoothly and gradually across this zone. Because of this demarcation, we refer to the structure as a “bud”. Interestingly, once such a bud is induced, it may sustain its growth independently after the F-source is removed. This is shown in Fig 5A, in which the exogenous source is shut off after the 12th cell cycle (see S13 Fig for the corresponding evolution of the cell distributions). In this case, bud growth continues unabated, with an outcome almost identical to that in Fig 4A (see S21 Fig for a comparison). The reason is that the growing bud itself becomes an effective source of F, and it is this endogenous feedback that drives most of the growth (see S14 Fig, which shows the distributions of P). If exogenous feedback is removed too soon, however, before the bud is large enough, endogenous feedback is insufficient to sustain growth (see S15 and S16 Figs). In particular, CPs differentiate into TD cells, and growth stops. Just as a transient exogenous positive signal can initiate a self-sustaining process of local epithelial growth, an exogenous source of G can also shut off growth. The simulation in Fig 5B begins in the same way as Fig 5A, with a transient source of F applied for 12 cell cycles, but following its removal a transient source of G is applied from the 15th to 18th cell cycle. Growth of the bud gradually arrests, as CPs differentiate into TDs, and the level of endogenous feedback needed to keep P above ½ is no longer sustained. See S17–S20 Figs for plots of the cell distributions, the distributions of P, and an analysis of strength of the exogenous source of G on the growth arrest. We also tested the effects of stochasticity on p and ζ as above, and found it more difficult to trigger growth arrest using a transient exogenous source of negative feedback in the stochastic case, i. e. larger amounts of G were required to stop growth. These sorts of results (which are summarized in S21 Fig) are reminiscent of the observation in non-spatial systems that a single system can be moved in or out of the “growth-leap” region simply through an appropriate change in initial conditions (Fig 2D). In Fig 5A, the exogenous source of F drives growth of a localized region of the tissue to a point where it has sufficient endogenous F to enter a growth-leap regime. Indeed, as shown in Fig 6, short pulses of exogenous sources of F or G can also be used in non-spatial simulations to toggle, in an analogous manner, between fast- and slow-growing states. In both cases, it is the inherent bistability/bimodality of the systems that ensures an ultrasensitive dependence of such toggling on the levels of F and G applied. In spatial systems (Figs 4 and 5), that ultrasensitivity plays out as a sharp demarcation in the region of the tissue that responds to an exogenous cue, even if that cue forms a shallow spatial gradient over a larger region. Spatially demarcated behavior persists even after removal of an exogenous source because endogenous factors also form gradients, and the responses to them are similarly ultrasensitive. The sorts of buds that appear in the simulations in Figs 4 and 5 are simple, expanding balloon-like structures, with relatively little spatial organization. One reason for this is that the domains being simulated are relatively small compared with the ranges over which diffusible factors spread. To explore whether more intricate structures can be formed, it is necessary to “zoom out”, either by increasing the size of the simulated domain or—equivalently—reducing the spread of the diffusing molecules. In Figs 7 and 8 we do the latter, decreasing the intrinsic decay lengths of F and G, and adding uptake of the negative feedback factor G in the stroma (equivalent to including stromal expression of an irreversible inhibitor, such as follistatin [25,38]). In Fig 7, we start with initial conditions from which growth does not spontaneously occur, and initiate growth through a point source of F (Fig 7A and 7E). In panels A-C the ϕ/γ ratio is set at a level that does not allow for self-sustaining growth, and therefore the exogenous source is left on throughout the simulation (as in Fig 4). In panels E-G, the ϕ/γ ratio is slightly higher, so that growth, once initiated, will self-sustain, and thus the exogenous source is removed after 15 cell cycles (for further details, see S22 and S23 Figs, which show the distributions of the cells and F and G). In both cases, buds initiate, but they continue to grow outward as fingers of relatively constant width (Fig 7D and 7H). Interestingly, this width is at least 10 times the characteristic decay lengths of either F or G, implying that it is determined by characteristics of the system other than the spatial ranges of action of diffusible signals. The ability of a self-sustaining finger to grow in one direction, while maintaining a constant width in the other, suggests some sort of spontaneous self-organizing ability. One way to test for spatial self-organization is to see whether noise (random small inputs) gets transformed into quasi-regular patterns that are independent of the details of the noise. For example, reaction-diffusion systems that generate Turing patterns, when initiated from any of a variety of noisy initial conditions, evolve to produce patterns of stripes or spots with characteristics that directly reflect the underlying parameters of reaction and diffusion. In our simulations, a simple way to mimic noisy initial conditions is to make many small perturbations to the initial shape of the modeled tissue. In Fig 8, small, sinusoidal perturbations of a variety of arbitrary spatial frequencies were introduced along the basal surface. Values of the ϕ/γ ratio were chosen so that self-sustaining growth would not normally occur. However, the shape irregularities introduced by the initial conditions create tiny locales in which endogenous F can rise above or fall below the average. Some of these reach the threshold for locally self-sustaining growth, with the result that a series of buds begin to emerge. Eventually, buds fuse to produce wide fingers with a relatively regular spacing. Fingers may continue to fuse, or branch to create new fingers, depending upon parameter values. Fig 8A–8C, 8H–8J and 8O–8Q show the epithelia morphologies together with the distributions of P for different ratios of ϕ/γ and initial conditions. The corresponding distributions of F (from Fig 8A–8C and 8H–8J) are shown in S24A and S24B Fig. As in Fig 7, fingers are many times wider than the characteristic decay lengths of factors F and G (see scale bar, upper right). We find that finger widths are most readily influenced by changes to the diffusion coefficients DF and DG, with the widths increasing as a function of DF and decreasing as a function of DG (see S24C Fig). If DF gets too large, however, growth fails to self-sustain, because peak levels of F become too small. To show that these patterns are truly self-organized, and do not merely represent the expansion of pattern already present in the initial conditions, we used Fourier transforms of contours extracted from the images in Fig 8D, 8K and 8R to follow spatial frequency components over time. As shown in the kymographs in Fig 8E, 8L and 8S, as well as in individual periodograms excerpted from them (Fig 8F, 8G, 8M, 8N, 8T and 8U), simulations initially display power at frequencies corresponding to those used in the (sinusoidal) initial conditions, but quickly redistribute power to frequencies in the range of 2–4 (i. e. periods of 1/2 to 1/4 of the width of the modeled domain, which in these cases is ~12. 6 times the characteristic decay length of diffusible factors F and G). This occurs regardless of whether the dominant frequencies in the initial conditions are high or low (Fig 8V and 8W; see also S24D Fig), and are independent of both the exact value of ϕ/γ and the details of the boundary conditions. What these data indicate is that the role of initial conditions in these simulations is merely to get growth started; afterward, spatial interactions organize the behavior of the system so that fingers of similar size are produced. Such self-organization, through spatial interaction, is reminiscent of the evolution of Turing patterns. [21,22,41,42]. As with Turing patterns, this behavior depends upon diffusible, antagonistic factors (“activators” and “inhibitors”). Here, however, the patterns emerge as a spatial manifestation of bistability and bimodality generated by the diffusible factors and cell displacements. Thus, unlike Turing patterns, there is no need for the diffusivities of the factors to be different—indeed the diffusivities of the positive and negative factors are equal in Figs 7 and 8, and the characteristic scales of the patterns are much larger than the characteristic decay lengths of these factors.
  Here we show that feedback regulation of stem or progenitor cell self-renewal can lead to growth bistability or bimodality, with ultrasensitive switching between high and low (or zero) growth states. A necessary condition for such behavior is the simultaneous presence of saturable positive and negative feedback. Exactly which lineage stage serves as the source of feedback signals appears less important than the fact that the strength of positive feedback must significantly exceed that of negative. When this type of arrangement is implemented in space, the resulting ultrasensitivity can lead to sharply demarcated growth differences that may be self-sustaining, producing bulges, buds or fingers that could serve as building blocks for morphogenesis. Such behavior typifies self-organization, in that the forms and sizes of buds and fingers that emerge reflect the parameters and structure of the feedback system, not the details of the stimuli that initiated them. Because we only explored here a small number of very simple lineage architectures, boundary conditions, and sources of noise; and we considered only a limited number of mechanical constraints on cell behavior, it is possible that the diversity of forms that could be generated through lineage feedback is much greater than shown in this study. Even for the relatively simple buds and fingers that we observe here, it will be difficult to determine precisely what system characteristics are the primary determinants of shape and size. This is because the decision by cells whether or not to enter a “growth leap” regime at any point in time and space will be influenced not only by the parameters of the system, but also by the gradients of diffusible factors that are set up by dividing and differentiated cells nearby. So far, the idea that lineage negative feedback (i. e. negative feedback regulation of self-renewal) plays a major physiological role has come mainly from studies of size control in self-renewing tissues, although the molecular nature of the relevant feedback factors is known in only a subset of such cases [25,27,28,30,43]. In many of these tissues, positive feedback factors—secreted molecules that promote cell renewal—have been clearly identified, yet there has been little speculation about the purpose of such feedback. Here we raise the possibility that positive feedback is useful because it creates the potential for growth-driven, self-organizing morphogenesis. In support of this view, we note several in vivo reports of phenomena analogous to those discussed here, including growth bistability/bimodality, feedback ultrasensitivity, and feedback-regulated morphogenesis. For example, in the mouse olfactory epithelium (OE), a reduction in the strength of FGF signaling (due to loss of Fgf8 [44]), or an increase in the strength of GDF11/activin signaling (due to loss of Foxg1, which regulates responsiveness to GDF11 and activin [38]), can lead not just to a smaller OE, but to a complete absence of tissue (agenesis). Agenesis due to loss of Foxg1 can be rescued by inactivation of Gdf11, consistent with idea that a critical threshold ratio of positive-to-negative feedback determines whether a zero or non-zero growth state is achieved [38]. Interestingly, inactivation of even a single allele of Gdf11 in Foxg1 mutant animals restores the OE to full thickness, supporting the idea that growth state is an ultrasensitive function of feedback strength, as in Fig 1F. In mouse muscle, evidence for growth bistability can be found in transplantation experiments in which myofiber-associated satellite cells are allowed to re-populate injured muscle [45]. In such experiments the transplanted cells, which are FGF2-treated prior to transplantation, trigger an abnormally high rate of myoblast proliferation and differentiation, which becomes sustained, without further intervention, for years. As in the OE, muscle cell lineages are known to be controlled by mixed negative and positive feedback through secreted factors myostatin [28] and FGFs [46]. Triggered high rates of proliferation and differentiation that self-sustain for long periods are also observed in pathological skin conditions, such as psoriasis, acanthosis nigricans, and keloid formation [47–49]. Psoriasis is a particularly interesting example, because the excessive keratinocyte stem cell proliferation that characterizes psoriatic lesions typically occurs in domains that are sharply juxtaposed with normal skin, for reasons that are not yet understood. Ultrasensitivity that arises from a mix of intrinsic positive and negative lineage feedback could easily account for such sharply demarcated growth differences, as illustrated in Figs 4 and 5. Interestingly, therapeutic resolution of psoriatic lesions typically occurs only when stem cells are ablated below some threshold level [50], exactly what the lineage models presented here predict for a population that becomes “stuck” in a high growth state due to endogenous positive feedback. In acanthosis nigricans, which often occurs as a consequence of type II diabetes [48,51], the epidermis thickens due to the presence of deeply folded stem cell layers, producing structures much like those observed in Fig 8. It is thought that high insulin levels interact with insulin-like growth factor-1 receptors to inhibit epidermal differentiation, i. e. insulin acts as an exogenous signal that boosts P [52]. In Figs 4–8, just such a signal can clearly act as a trigger of budding and fingering. Of course, fingered growth also occurs in normal tissues, as exemplified by the palisades of Vogt of the corneal limbus; lingual papillae, hair follicles, and even intestinal crypts, which initially develop as downgrowths from a smooth epithelial lining [53–56]. The roles played by diffusible positive or negative factors in the morphogenesis of these structures are incompletely understood, although in several cases the existence of both positive and negative feedback has been documented [55,57]. Although these examples support the view that interesting, spatially dynamic phenomena do arise in vivo as a result of mixed positive and negative lineage feedback, the extent to which the morphogenesis that occurs during normal development relies upon such phenomena remains to be investigated. Some developmental phenomena do bear an intriguing superficial resemblance to observations made here. For example, vertebrate tailbud elongation proceeds through outgrowth, at a relatively constant rate, of a constant-width, finger-like projection, with self-renewing cells concentrated at the tip [58,59], much like what is observed in Fig 7. Moreover, tailbud outgrowth is known to be driven by high concentrations, at the growing tip, of renewal-promoting FGF [60,61], and terminated by differentiation-promoting retinoic acid, a diffusible factor that is both exogenously and endogenously provided, and acts in an ultrasensitive fashion to extinguish growth [61]. Obviously, feedback regulation of lineage progression is only one of many morphogenetic mechanisms at the disposal of developing tissues. Processes such as regulated cell cycle progression, directed cell migration, cell sorting, and active mechanical deformation clearly play important roles as well [e. g. 7,17–20,62–64]. What the present study elucidates are some of the advantages of lineage feedback-driven morphogenesis. These include the ability to specify spatial structure using the same molecules that regulate tissue homeostasis; the ability to coordinate the behaviors of cells over reasonably large spatial ranges; and the ability to create structure through self-organization, i. e. to form borders and boundaries spontaneously without the need for templates or pre-patterns.
  Analyses of the spatially homogeneous, ordinary differential equation models were performed using Wolfram Mathematica. Simulations of these models were performed using NDSolve. Stochastic simulations were carried out using MATLAB. The volume fractions of CP (χ1) and TD (χ2) satisfy the mass conservation equation ∂χi∂t=−∇⋅Jχi+Srcχi−∇⋅ (uχi), (4) where i = 1,2. The flux Jχi accounts for interactions among cell species, and u is the mass-averaged cell velocity. Here we have assumed that the cell densities are matched, that the cells are tightly packed and that the cells move with the same velocity following [39]. The source term Srcχi quantifies the net effects of the proliferation, differentiation and death of the lineage components: Srcχ1=v0χ0+v1 (2p1−1) χ1, Srcχ2=2v1 (1−p1) χ1−dχ2 (5) The self-renewal fraction of the CPs, p1, is regulated by a soluble negative feedback factor G (e. g. , GDF) and a positive feedback factor F (e. g. , FGF) instead of directly by the CP and TD populations in the spatially homogeneous model. Defining [G] and [F] as the concentrations of G and F respectively, we take p1=p1, min+ (p1, max−p1, min) φ[F]1+φ[F]11+γ[G], (6) which is analogous to the feedback in the ODE model. Here P1, min, P1, max are the minimum and maximum levels of self-renewal, and ϕ, γ are the positive and negative feedback gains respectively. The concentration [F] satisfies the quasi-steady reaction-diffusion equation 0=∇⋅ (DF∇[F]) +sF1χ1+sF2χ2− (dF+uFχ1) [F]+ExtSrcF (7) The equation for [G] is identical, but with different coefficients. We use a quasi-steady equation here because the diffusion time scales of G and F (e. g. minutes) are much faster than cell proliferation time scale (e. g. hours). Since it is not known which cell types produceF, the production of F may come from CPs or TDs or both. In particular, sF1 and sF2 are the production rates from CPs (χ1) and TDs (χ2) respectively, dF is the natural decay rate, uF is the uptake rate by CPs and DF is the diffusion coefficient. The last term represents the effects of an exogenous source. In the equation for [G] EXTSrcG also includes a sink that models the binding of G by Follistatin (FST) [65]. A complete list of dimensionless parameters and their values may be found in S1 Table. Following [39,66,67], we use a simplified model of soft tissue mechanics to obtain the cell velocity and flux. In particular, we use a generalized Darcy law to describe the passive motion of cells due to proliferation-induced pressure. This assumes that the tissue can be treated as a viscous, inertia-less fluid and also models flow through a porous media. Other constitutive laws and model formulations may be found in [39,66–72] and are expected to yield similar results. We take the cell velocity to be u=−κ (∇p−λ εμ∇χT), (8) where κ is a non-negative cell motility, p is the pressure generated by cell proliferation, μ is the chemical potential defined as μ= (dfdχ (χT) −ε2∇2χT) (9) where ε measures the thickness of the basement membrane and apical surfaces (both assumed to have the same thickness for simplicity). The parameter λ models the surface tension (stiffness) of the epithelia interfaces. The surface tension may be different for the BM and AP surfaces. The function χT denotes the epithelium volume fraction: χT = χ1 + χ2, and f (χ) =14χ2 (1−χ) 2 is a double-well potential. Note that the term λμε∇χT≈λBMHBMδΣBM+λAPHAPδΣAP where λBM, λAP are the surface tensions of the BM and AP surfaces, HBM, HAP are the corresponding mean curvatures of the surfaces and δΣBM, δΣAP are surface delta functions localized on the basement membrane and apical surfaces, e. g. , see [39]. Both terms in Eq (8) model passive cell motion with the first due to differential proliferation and the second due to cell-cell adhesion. Finally, the flux Jχ is modeled as a generalized Fick’s law [39] Jχ=−Mχ∇μ, (10) where M represents the magnitude of the flux. Together, these constitutive relations for the cell velocity and fluxes guarantee that the model is thermodynamically consistent. That is, the adhesion energy Eadhesion=1ε∫λ (f (χT) +ε22|∇χT|2) dx≈λBMLBM+λAPLAP, where LBM, LAP are the lengths of the BM and AP surfaces, is non-increasing in the absence of net sources, where the integration is taken over a domain containing the tissue, see [39]. To solve for the pressure, we assume that the mass conservation equation also holds for the volume fractions of the stroma (χs) and the fluid domain (χL) exterior to the AP surface. Defining χNE = χs+χL, we assume 1 = χT+χNE and correspondingly that JχNT=JχS+JχL=− (Jχ0+Jχ1+Jχ2). We also assume that there are no net sources in the stroma or fluid domains, i. e. Srcχs = SrcχL = 0. Then, summing the volume fraction equations for the lineage components, the stroma and the fluid, we obtain ∇⋅u=Srcχ0+Srcχ1+Srcχ2=SrcχT, (11) from which it follows that −∇⋅ (κ∇p) =SrcχT−∇⋅ (κλεμ∇χT) (12) We solve the system of equations in a rectangular computational domain Ω = { (x, y) |−20 ≤ x 20, −10 ≤ y ≤ 10} that contains the epithelial tissue. We use homogeneous Neumann boundary condition for cell types χ1, χ2 and feedback factors [F], [G]. In order to allow the tissue to leave the computational domain smoothly, we apply homogeneous Dirichlet boundary conditions for the pressure p and chemical potential μ, following [39]. The no-flux boundary condition at the AP is implemented by the diffuse domain approach (see section 6 of S1 Text). The stem cell population χ0≈A (t) δΣBM, where A (t) is chosen to ensure that numbers of SCs are maintained, even if the BM changes shape. The details about this approach are also found in section 6 of S1 Text. To solve the governing equations efficiently, we follow previous work in [66] and use an adaptive finite-difference method with an implicit time discretization. The equations at the implicit time step are solved using the nonlinear multigrid method, where the equations are reformulated as a system of second order equations. Additional details may be found in section 6 of S1 Text. The spatial distributions of cells and feedback factors are visualized using Matlab. Let T represent the current time in the simulation. Suppose that the exogenous source for F is introduced from time TF1 and TF2. Let ExtSrcF=αF (AF−[F]) χS (13) for TF1≤T≤TF2. Here χS is the volume fraction of the stroma, AF is a prescribed Gaussian source centered at the lower left corner of the computational region. F is produced in the stroma at a rate proportional to its difference from a target concentration AF with a proportionality constant αF. Given the computational domain { (x, y) |−20 ≤ x ≤ 20, −10 ≤ y ≤ 10}, we choose AF=MFexp{− (x+20) 210− (y−yC−1) 210}, (14) where MF measures the magnitude of the exogenous concentration, yC is the leftmost position of the BM. The feedback factor G has an uptake term in the stroma. Denoting uG, s as the uptake rate, we define ExtSrcG=−uG, S[G]χS (15) The definition of exogenous source of G is analogous. Suppose that the exogenous source for G is introduced from T=TG1 to T=TG2, we define ExtSrcG=−uG, S[G]χS+αG (AG−[G]) χS (16) for TG1≤T≤TG2. Similar to F, G is produced in the stroma at a rate proportional to its difference from a target concentration AG=MGexp{− (x+20) 250− (y−yC−1) 250} (17) with proportionality constant αG Here MG measures the magnitude of the exogenous concentration of G, yc is the leftmost position of the BM. The values for these dimensionless parameters are listed in S1 Table.
summary: A large variety of mechanisms have been proposed to explain how the cells of solid tissues work together to produce distinctive structures such as buds, branches, sacs and tubes. Here we identify a novel mechanism that emerges directly out of the processes that cells use to control proliferation and maintain overall tissue sizes. In particular, we focus on the diffusible molecules that regulate lineage progression—the decision by dividing cells either to renew themselves or produce cells of a more differentiated type. When such molecules are produced by cells at specific lineage stages, they create feedback circuits, and positive and negative lineage feedback circuits are often observed together in the same tissue. Here we show that the combination of positive and negative lineage feedback can lead to the emergence of patterned growth from which sharply-demarcated buds and fingers readily emerge, either spontaneously or in response to transient external signals. We show that the types of behaviors that emerge, and the conditions under which they do so, depend strongly on the relative strengths of positive and negative feedback.
keywords: cell death
  medicine and health sciences
  mechanisms of signal transduction
  cell cycle and cell division
  cell processes
  cell differentiation
  developmental biology
  stem cells
  morphogenesis
  feedback regulation
  animal cells
  muscle differentiation
  biological tissue
  signal transduction
  anatomy
  cell biology
  epithelium
  biology and life sciences
  cellular types
title: Feedback, Lineages and Self-Organizing Morphogenesis