Feedback control in planarian stem cell systems

Background In planarian flatworms, the mechanisms underlying the activity of collectively pluripotent adult stem cells (neoblasts) and their descendants can now be studied from the level of the individual gene to the entire animal. Flatworms maintain startling developmental plasticity and regenerative capacity in response to variable nutrient conditions or injury. We develop a model for cell dynamics in such animals, assuming that fully differentiated cells exert feedback control on neoblast activity. Results Our model predicts a number of whole organism level and general cell biological and behaviours, some of which have been empirically observed or inferred in planarians and others that have not. As previously observed empirically we find: 1) a curvilinear relationship between external food and planarian steady state size; 2) the fraction of neoblasts in the steady state is constant regardless of planarian size; 3) a burst of controlled apoptosis during regeneration after amputation as the number of differentiated cells are adjusted towards their homeostatic/steady state level. In addition our model describes the following properties that can inform and be tested by future experiments: 4) the strength of feedback control from differentiated cells to neoblasts (i.e. the activity of the signalling system) and from neoblasts on themselves in relation to absolute number depends upon the level of food in the environment; 5) planarians adjust size when food level reduces initially through increased apoptosis and then through a reduction in neoblast self-renewal activity; 6) following wounding or excision of differentiated cells, different time scales characterize both recovery of size and the two feedback functions; 7) the temporal pattern of feedback controls differs noticeably during recovery from a removal or neoblasts or a removal of differentiated cells; 8) the signaling strength for apoptosis of differentiated cells depends upon both the absolute and relative deviations of the number of differentiated cells from their homeostatic level; and 9) planaria prioritize resource use for cell divisions. Conclusions We offer the first analytical framework for organizing experiments on planarian flatworm stem cell dynamics in a form that allows models to be compared with quantitative cell data based on underlying molecular mechanisms and thus facilitate the interplay between empirical studies and modeling. This framework is the foundation for studying cell migration during wound repair, the determination of homeostatic levels of differentiated cells by natural selection, and stochastic effects. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0261-8) contains supplementary material, which is available to authorized users.


Background
Stem cell systems operate by demand control [1][2][3] in which the needs of the organism determine in large part the behavior of the stem cells. Indeed, both cancer and ageing may be understood as failures of this feedback control, albeit in different ways. The highly regenerative level of the single gene to the entire organism. The planarian life history provides the unique opportunity to take a systems approach to understanding stem cell dynamics in a whole organism.
In planaria, stem cells are called neoblasts and are defined collectively as the only dividing cells in the animal. Among these cells it has long been assumed that at least some cells are bona-fide pluripotent stem cells (see [9] for the most up to date review), capable of indefinite self-renewal and of producing all differentiated cell types in the adult animal; this was recently experimentally verified in the model species Schmidtea mediterranea [12]. A growing body of co-expression data shows that sub-populations of cycling neoblasts express lineage specifc mRNA markers [13]. Some of these co-expressed markers are functionally required for production of both the neoblast sub-population and the differentiated cell lineage in question; reviewed in [14]. This provides evidence for the existence of committed proliferating cells amongst the neoblast population but still awaits definitive experimental proof.
Fully differentiated cells in planarians have been divided into about 15 different classes, or 3 to 5 super-classes (e.g. cells associated with metabolism, muscle, nerve, and the epidermis), with the actual number of functional cell types likely to be much higher [8,15]. Unlike other stem cell systems such as the bone marrow stem cell system, in planaria there is still no conclusive evidence for mitotically active progenitor cells with strictly limited potency [16][17][18]. There are however populations of transient postmitotic stem cell progeny, and these cells either differentiate to a target lineage or potentially may apoptose rather than complete differentiation. We assume that the proportion of the various types of differentiated cells is regulated towards a homeostatic target [19,20] but in this paper do not model how that target is set (see [3] and "Conclusions" here).
The requirement for a given mix of differentiated cells and the simplicity of the system make planarians an ideal system for studying homeostasis and regeneration, including scaling, reproductive fission, and responses to wounding and amputation [21,22]. A minimum remaining tissue size is needed for such recovery after artificial amputation [23], but from that minimum size the entire organism can be regenerated through proliferation of stem cells to produce new tissue de-novo and remodeling of the remaining tissue [24].
We take a dynamical view of planarian stem cell system [25,26] and thus formulate our models as nonlinear dynamical systems in which the nonlinearity arises through feedback control. Although intraspecific variation and even intra-organism variation [9] exist, we ignore it for now.
We develop two models based on current knowledge of planarian stem cell biology. The first has neoblasts, nonmitotic progenitor cells, and three kinds of differentiated cells (three chosen for simplicity of setting parameters; the methods scale readily for an arbitrary number of differentiated cells). In the second model, we assume only one kind of differentiated cell and that the progenitor to differentiated cell transition is essentially instantaneous. This allows simplification that makes presenting some results clearer without losing any general principles.
We next give verbal and pictorial descriptions of the models, which are fully formulated in the "Methods" section. We then turn to Section "Results" beginning with steady state prediction of planarian size in relation to food in the environment and show that both the strength of feedback control in response to food in the environment and the constancy of the fraction of neoblasts are emergent properties of the model. We use the full model to explore growth, shrinkage, and regrowth under sufficient resources to maintain metabolism and to explore bursts of cell activity during remodeling following a fission or wounding. We use the simplified model to repeat the study of growth, shrinkage, and regrowth and then consider two in silico experiments. In the first experiment, we 'wound' the planarian by removing a large number of differentiated cells, while in the second experiment we simulate death of neoblasts as happens after γ irradiation [27]. Following this, in the "Discussion" section we focus on the comparison of our model to extant data and the potential for future experimental data. In the Section "Conclusions" we summarize the predictions of the model and look forward to future developments.

Verbal and Pictoral description of the models
In the Section "Methods" we describe the models as a discrete time dynamical system (obtained by writing a system of differential equations in a form suitable for numerical solution). Here we give both verbal and pictorial ( Fig. 1) descriptions of the model. In Fig. 1a we show the cell types and their transitions, along with the resource pool Q. We assume that neoblasts, N, are immortal and do not directly transition to differentiated cells, but rather go through a progeny P cell that is no longer mitotically active but which is not fully differentiated. Such a progeny cell may complete differentiation to one of the kinds of fully differentiated cells (in the figure, we show three kinds of such cells, indexed as D i ) or may apoptose and thus return to the resource pool. Thus, the transitions for neoblasts are N → N, N; N → N, P; and N → P, P. The transitions for progeny cells are P → Q and P → D 1 , P → D 2 , or P → D 3 .
These transitions have a maximum rate that is modified by feedback control. In Fig. 1b, we show the positive feedback control on these transitions. Neoblasts exert a b c Fig. 1 a Pictorial description of the model involving neoblasts (N), non-mitotic progenitor cells (P), three classes of differentiated cells (D i ), and a resource pool. The maximum rates of asymmetric renewal (N → N, N), symmetric differentiation (N → N, P), and asymmetric differentiation (N → P, P) are p 1 , p 2 , and p 3 respectively. As explained in the Section "Methods", a steady state exists only if p 1 < p 3 . Progenitor cells either complete differentiation or return to the resource pool. In the absence of external injury or death fully differentiated cells apoptose and return to the resource pool. These transitions have a maximum rate that is modified by feedback control. b Positive feedback controls superimposed upon the transitions. Neoblasts exert positive feedback f N (N) on the N → P, P transition in the sense that as the number of neoblasts increases, f N (N) increases, bounded by 0 and 1. The resource pool Q also has positive feedback on all transitions, in the sense that with larger resource pools the rate of transition is higher. However, we assume that the resource pool operates differentially on the different transitions, so use to indicate the feedback control of the resource pool on the transitions N → N, N, N → N, P, and N → P, P respectively. c Differentiated cells exert negative feedback control on the transitions, in the sense that as the number of differentiated cells increases, the rates of transitions decline, sharing a common feedback control f D (D).
Here we assume the that the differentiated cell in shortest supply sets the feedback control. Absent an external source of mortality, the only transition for differentiated cells is D i → Q through cell death, which occurs for cell type i at rate M i (N, D, Q). In addition, progenitor cells may either fully differentiate or return to the resource pool through apoptosis. We assume that the rate of the former is determined by a function f Q (Q) that increases as the size of the resource pool increases positive feedback f N (N) on the N → P, P transition in the sense that as the number of neoblasts increases f N (N) increases; it ranges between 0 and 1. The resource pool Q also has positive feedback on all transitions, in the sense that with larger resource pools the rate of transition is higher. However, we assume that the resource pool operates differentially on the different transitions, so use a 1 (N, D, Q), a 2 (N, D, Q), a 3 (N, D, Q) where D is the vector (D 1 , D 2 , D 3 ) to indicate the feedback control of the resource pool on N → N, N, N → N, P and N → P, P respectively. In Fig. 1c, we show the negative feedback control that differentiated cells exert on the transitions, in the sense that as the number of differentiated cells increases, the rates of transitions decline, sharing a common feedback control f D (D). We assume the that the differentiated cell in shortest supply sets the level of feedback control (details given in Section "Methods").
In the absence of an external source of mortality, the only transition for differentiated cells is D i → Q through cell death, which occurs for cell type i at rate M i (N, D, Q). In addition, progenitor cells may either fully differentiate or return to the resource pool through apoptosis. We assume that the rate of the former is determined by a function f Q (Q) that increases as the size of the resource pool increases.
In  (Fig. 1a, Fig. 2b) is a sigmoidal function that increases towards 1.0 as the number of neoblasts (measured as the fraction of the steady state value) increases. It is characterized by two parameters: the value of neoblasts at which f N (N) = 0.5 and the spread around that value. For a model with three kinds of differentiated cells, the feedback functions a i (N, D, Q) and the resource dependent rate of mortality of differentiated cells depend on five variables. In Fig. 2c, we show a cross-section of those functions by holding the number of neoblasts and differentiated cells constant and only varying the resource level. The values shown here are illustrative of the shape and relationship of the three feedback functions. The x-axis intentionally has no units since these images are intended to be schematic rather than accurate representations. In Fig. 2d we show the resource dependent rate of natural mortality, again in cross-section. Full details are in the Section "Methods".
To capture the transition from progenitor to fully differentiated cell, we make the relative need assumption that progenitors transition to differentiated cells according to how far they deviate from the homeostatic level, which is determined by both the number of differentiated cells and their relative proportion. The rate of apoptosis of differentiated cells is also determined by the number of differentiated cells and how far they deviate from the homeostatic proportion, and the resource pool. In particular, the rate of mortality increases for numbers of, and/or proportions of, differentiated cells above their homeostatic level and also increases as the resource pool declines. In the simplified model, we compress all of the differentiated cells into a single type and assume that the progenitor transitions are so rapid that they can be ignored. This allows clearer analytical and pictorial representation of the feedback controls. Details are given in Section "Methods".

Overview
We begin with steady state results, using the full model, showing the relationship between food in the environment and planarian size, and how the strength of feedback control from differentiated cells to neoblast activity emerges in response to food in the environment. We also demonstrate that in the steady state the fraction of neoblasts is independent of size (cf [28], and the "Discussion"). We then use the full model to study growth, shrinkage, and regrowth with sufficient resources. When resources are ample, feedback control is independent of resources and only depends on the number of differentiated cells (i.e., the a i (N, D, Q) = 1 and f Q (Q) = 1). We then turn to remodeling of a planarian following fission. To do, this we assume that the initial cell numbers are a small fraction of their steady state values and that neoblasts and differentiated cells differ from their homeostatic levels. Thus we anticipate that as the planarian regenerates it will need to change the number of neoblasts and mixture of differentiated cells (decreasing some that remain in excess and increasing those in shorty supply) so that there will be a burst of mortality following fission, increasing the resource pool which is then used to regrow towards the steady state. We use the simplified model to study cell activities during growth, shrinkage, regrowth, and regeneration, and in the in-silico experiments in which a fraction of the differentiated cells are removed, as would happen with wounding or amputation, or a fraction of the neoblasts are removed, as would happen with a relatively precise x-ray treatment at the center of the planarian.

Steady states
In Fig. 3 we show the steady state size of the planarian for the full model; this is determined by the steady state number of cells of all types and an allometric relationship between size and total cell number (cf. [29] and the "Discussion"). The curvilinear nature of this relationship is due to not all cells being able to acquire resources from the external environment (Eq. 12 below).
As described in detail in the methods, the strength feedback control emerges from the cell dynamics, depending upon the level of food in the environment and the maximum rates of asymmetric renewal and differentiation (Fig. 4). The strength of feedback control from differentiated cells to neoblast activity declines as food in the environment increases because the increasing food supports more cells overall. Thus, we predict that the strength of the signaling system between differentiated cells and neoblast activity will respond to the level of external food.
In Eqs. 13-21, we show that the fraction of neoblasts in the steady state is constant, independent of the size of the planarian. This proportion is determined by the transition rates of neoblasts and rates of death of differentiated cells (see the Section "Discussion" for comparison with data).

Dynamics with sufficient food resources (full model)
In Fig. 5, we show the temporal pattern of size, total mortality of differentiated cells, total mortality of differentiated cells, and fraction of neoblasts under growth, shrinkage, and regrowth for a situation in which the planarian starts at 20 % of its steady state size in an environment with food availability Y e = 450. Shrinkage can occur with sufficient resources due to apoptosis of differentiated cells that return to the resource pool. At scaled time 7 1 , food is dropped to Y e = 400 and then at scaled time 14 increased to Y e = 480. Even in the absence of a resource constraint, the cell population dynamics show a dependence upon the level of food in the environment. This is due to the feedback control on neoblast activity from differentiated cells and feedback control of neoblasts on themselves in asymmetric division.

Remodeling following a fission (full model)
Most organ systems in planarians are broadly distributed throughout the body [9] but not uniformly so. Thus, after a natural fission, the neoblast proportion is probably slightly higher than the normal steady state, but many differentiated cell types will be much lower than steady state and others will be higher. For example, in tail fragments neurons will be lower than the homeostatic level and gut cells higher. In Fig. 6, we show cell dynamics during remodeling following a fission. We assume that the initial cell numbers are 50 % of the steady state number of neoblasts, and 10 %, 30 %, and 5 % of the steady state values of the three kinds of differentiated cells respectively (rather than 25 % neoblasts and 40 %, 30 %, 30 % relative distribution of differentiated cells in the steady state). Starting from relatively small size, the planarian grows to a steady state (since food is constant), has a burst of mortality as the differentiated cells apoptose in order to achieve their target values (cf [30] and the "Discussion").

In-silico experiments (simplified model)
We first repeat growth, shrinkage, and regrowth using the simplified model (Fig. 7) to verify the same qualitative patterns. We started the planarian at 20 % of the steady state size associated with food level Y e = 240 and then varied food over time as shown in Fig. 7a. The planarian grows towards its new steady state until food is decreased at scaled time t = 7 and then again at scaled time t = 15 and in response the planarian size decreases (Fig. 7b). During this entire process, however, the fraction of neoblasts is nearly constant (Fig. 7c). The feedback controls f D (D) [note that there is only one kind of differentiated cell in the simplified model] and f N (N) (Fig. 7d and e respectively) respond to the food pattern in complex ways, in part because α in the feedback control f D (D) = 1 1+αD is a function of external food level (Fig. 4) and because of additional mortality when resources are insufficient.
In Fig. 8, we show the results of the in silico excision experiment in which we allow the planarian to grow to its steady state size for the given food environment and then at scaled time t = 7 remove differentiated cells via excision or wounding. This leads to an increase in the fraction of neoblasts (Fig. 8b), an increase in the feedback control function f D (D) (Fig. 8c) and a slight drop in f N (N) (Fig. 8d, note the vertical scale). Also note the different time scales in the recovery of size and the two feedback functions. Size is recovered by about scaled time t = 8.5 and although the feedback function f D returns to its previous value at approximately the same time, the feedback function f N takes much longer to return to its steady state. (cf [31] and "Discussion").
In Fig. 9 we show the results of a similar x-ray experiment: After the planarian has reached its steady state size at scaled time t = 7 neoblasts are removed. Since we assume that size is determined by differentiated cells only there is no change in size (see Section "Methods"), but an increase in the feedback control function f D (D) (middle panel), and a significant drop but then recovery of f N (N). Notice the very different patterns of the feedback control functions in Figs. 8 and 9, suggesting that different patterns of cell-signalling emerge from different in silico experiments and that these can be predicted.  Cell dynamics during remodeling following a division. We assume that the initial cell numbers are 50 % of the steady state number of neoblasts, and 10 %, 30 %, and 5 % of the steady state values of the three kinds of differentiated cells respectively (rather than 25 % neoblasts and 40 %, 30 %, 30 % relative distribution of differentiated cells in the steady state). Starting from relatively small size, the planarian grows to a steady state, since food is constant) (first panel), experiences a burst of mortality initially (second panel) as the differentiated cells apoptose in order to achieve their target values (fourth panel; colors are red: ρ 2 , black: ρ 1 , and blue: ρ 3 ), and a complicated trajectory for the fraction of neoblasts as the animal readjusts its composition

Comparison of the model to extant data and and the potential for future experimental data
Our eventual goal is to use experimental data to both test and expand the current model. Indeed one significant advantage of using the planarian experiments is that it is now quantitatively possible to assess stem cell self-renewal and differentiation events in vivo in different scenarios. To date, most experimental studies are qualitative or semiquantitative; we can use them to see if the current model can satisfactorily explain them or requires further innovation. While the current model has been developed with many biological observations in mind, such as the steady degrowth observed during starvation of planarian worms [15,29], the existence of large populations of post-mitotic progeny and that cell divisions N → N, P or N → P, P appear to increase during starvation [29], we ignored some other extant observations. Some of these, such as the consistency of the proportion of neoblasts irrespective of body size [15,29], are nonetheless predictions of the current model. Others are only partially represented or are currently absent and highlight areas for future refinement. For example, 4 h after amputation there is a large increase in apoptosis (cell death) close to the site of amputation [20]. A smaller and spatially broader pulse of cell death is then observed at 72 h after wounding/amputation  [20]. Both these events are thought to remove differentiated cells that are now inappropriately placed with respect to the missing tissue that needs to regenerate and remodel. Our model's behaviour allows for a similar adjustment of differentiated cells after fission (equivalent to amputation) and, in the absence of a spatial dimension in our model, reflects the same requirement to remove differentiated cells that are now over-represented. Furthermore, it is well established that neoblasts increase their proliferation rate high above basal levels in response to wounding forming  two temporally distinct proliferative peaks [30]. The first is spread spatially through the whole body between 6-12 h after wounding and the second is localised at the wound site at 48 h. This is an example of an observation that is currently not captured in the model, but can be incorporated into future iterations that expand on the details of the precise relative timing of stem cell behaviours and incorporate quantitative cell level data collected with this purpose in mind.
Genome wide studies of gene expression in planarian cells are now de rigueur [13,[32][33][34][35]; most recently this has included study of single cells [13] and comparisons between wildtype animals and RNAi loss of function phenotypes [36][37][38][39]. Many loss of function phenotypes provoke key changes in stem cell behaviour that can be monitored and subsequently correlated with gene expression changes. Our cell level model can be expanded to include the results of these experiments, particularly those that appear to impact transition states and feedback functions of in our model, or suggest the need to add additional dynamics. Thus our model will help experimental groups interpret the cell level phenotypes in the context of the stem cell behaviours. Eventually it should be possible to assign gene expression profiles to cell types and describe the genetic control of transitions and feedback controls in terms of quantitative gene expression profiles.
It is also possible to test the assumptions used in the model, which are described in detail below. These include i) determining the power of the Hill coefficient in the resource-dependent feedback control; ii) determining how the rate of apoptosis of fully differentiated cells depends upon their numbers and proportions; iii) determining how the rate of apoptosis of progeny that are not fully differentiated depends upon resources; iv) measuring the strength of feedback control from differentiated cells to neoblast activity; and v) determining the prioritization of resource use for cell divisions. The wound signal in planaria is currently unknown [9]; our model provides a conceptual framework (Eq. 2 below) relating the number of differentiated cells and the strength of that signal.
Thus, our current model predicts some extant experimental observations and more importantly provides a systems level framework within which to incorporate precise quantitative measurements spanning from the whole organism, through organs, tissues, cell types and eventually gene expression and function.

Conclusions
In summary, our models lead to testable predictions about the dynamics of size, the rates of mortality, and the signaling systems in planaria. These include: 1. a curvilinear relationship between external food and planarian steady state size; 2. the strength of feedback control from differentiated cells to neoblasts (i.e. the activity of the signaling system) and from neoblasts on themselves depends upon the level of food in the environment; 3. the fraction of neoblasts in the steady state is constant regardless of planarian size; 4. planarians adjust size when food shifts first due to apoptosis and then through a reduction in neoblast activity; 5. a burst of mortality during regeneration as the number of differentiated cells are adjusted towards their homeostatic level; 6. following wounding or excision of differentiated cells, different time scales characterize the recovery of size and the two feedback functions; 7. the temporal pattern of feedback controls differs noticeably during recovery from a removal of neoblasts or a removal of differentiated cells; and 8. the signalling for apoptosis of differentiated cells depends upon both the absolute and relative deviations of differentiated cells from their homeostatic level. 9. a whole-organism prioritization of resource use for cell maintenance, neoblast division, and progeny differentiation.
Much remains to be done, including comparing models to data explicitly, evolutionary origins of the steady state fraction [40], spatial dynamics of cells [41], and extension to stochastic models [42]. For example, the equations that we have developed here can be used in the formulation of a probabilistic model of homeostasis and regeneration [9], in which the deterministic framework emerges as the conditioned average of the stochastic system (see [43] for stem cells and [44] for the more general situation). A model such as the one developed here is the first step towards a full conceptual framework for planarian cellular dynamics and will complement the outstanding questions raised in [9], such as what is the wound signal stimulating neoblast differentiation, how do neoblasts sense the state of the organism, and what cellular machinery signals to the neoblasts.

Methods
We first describe the states that characterise the planarian, feedback control, and the dynamics of states. After that we determine the steady state and then discuss dynamics, and how in silico experiments can be performed to match empirical studies.

States and transitions
We let N(t) denote the number of neoblasts at time We denote the number of non-mitotic but undifferentiated progeny and resource pool at time t by P(t) and Q(t) respectively. We use lower case to indicate specifc values of these dynamic variables. We consider three kinds of transitions: 1) asymmetric neoblast renewal, i.e. N → N, N; 2) symmetric renewal and progeny production, i.e. N → N, P ; and 3) asymmetric progeny production, i.e. N → P, P. A progeny cell either returns to the resource pool or continues to complete differentiation into one of the types of differentiated cells, as explained below.

Feedback control
Each of the transitions 1-3 above are subject to feedback control [3,26,42]. Suppose that N(t) = n, D(t) = d and Q(t) = q.
We assume that there is a priority of resource use by the organism in which top priority is given to maintenance of existing cells and then resources are allocated to N → N, N, N → N, P, and N → P, P transitions if they are sufficient. Thus, we assume that the fraction of neoblasts undergoing the N → N, N transition is p 1 a 1 (n, d, q) where p 1 is a fixed value (we explain below how it is set, in the section on parameters) and a 1 (n, d, q) is determined as follows. If m r is the metabolic rate of neoblasts and differentiated cells (which we assume, for simplicity, to be the same) and m d is the cost of division, then q 1 = m r (n + I i=1 d i ) is the metabolic cost of maintenance of the existing cells and q 2 = q 1 +m d p 1 n is the level of resources needed to maintain all existing cells and support all asymmetric neoblast renewals. We let q 12 denote the average of q 1 and q 2 and model the feedback control on N → N, N divisions by a Hill-type function as long as q > q 1 ; otherwise we set a 1 (n, d, q) = 0. Whether the exponent of the Hill function is 2 is currently unknown and a topic for possible future experimental work as described in the Section "Discussion". The second term in the denominator of the right hand side Eq. 1 is clearly q 2 12 but we write the expression as above (and below for the other resource-dependent feedback controls) to make clearer that a 1 (n, d, q) starts at 0 when q = q 1 , reaches 0.5 when q = q 12 and asymptotes at 1. The transitions N → N + P and N → P, P involve resource-dependent and cell-number dependent feedback control. We assume that the fraction of neoblasts undergoing N → N + P is p 2 a 2 (n, d, q)f D (d) where p 2 is a  fixed value, a 2 (n, d, q) is determined in a manner similar to above and where α i sets the strength of control from differentiated cells of type i to neoblasts. The "max i " means that the differentiated cells that are in most demand set the level of feedback control. Equation 2 is consistent with the long held hypothesis that differentiated cells produce specific factors inhibiting their own growth (see [9], pg 7). Given the control from differentiated cells, q 3 = q 2 + f D (d)m d n is the level of resources needed to support all N → N + P transitions. We let q 23 denote the midpoint of q 2 and q 3 and set as long as q > q 2 ; otherwise we set a 2 (n, d, q) = 0. We assume that N → P, P involves an additional feedback control from neoblasts; so that when neoblast numbers are low this transition is suppressed. We set Since the exponent will be 0 when n = N c , f N (N c ) = 0.5. The parameter σ N controls the sigmoidal or S-shape of f N (n) . As σ N declines, f N (n) becomes more and more knife-edged, close to 0 when n < N c , close to 1 when n > N c but still 0.5 when equality holds. In the limit that σ N is very large (i.e. many times greater than n could be) f N (n) is close to 0.5 regardless of the value of n. We assume that σ N and N c are proportional to the steady state number of neoblasts, and are thus also environmentally determined by the level of food (see below, Parameters).
The two feedback functions f D (D) and f N (N) can be viewed as the result of transcriptional processes associated with homeostasis and regeneration/remodeling respectively (cf. [9]).
Thus, the fraction of neoblasts undergoing N → P, P transitions is p 3 a 3 (n, d, q)f D (d)f N (n) where p 3 is a fixed  value, and a 3 (n, d, q) is determined in a manner similar to above. That is, we set q 4 = q 3 + f D (d)F N (n)m d p 3 n, q 34 to be the average of q 3 and q 4 , and if q > q 3 and 0 otherwise. The correspondence between one unit of time in the model and physical time is set by the activity of neoblasts. In particular, the maximum fraction of neoblasts active in one unit of time is p 1 + p 2 + p 3 , which happens with abundant resources and feedback from differentiated cells and neoblasts both equal to 1. Thus a measurement of that fraction over a short interval of time provides a link between the cell cycle and the physical meaning of the time unit in the model. For results we have scaled time so that the numbers on the x-axis are 1/10000 th of what they are in the model.

The dynamics of neoblasts
We write the dynamics as difference equations rather than differential equations for two reasons. First, even if written as differential equations, the dynamics have to be solved numerically, requiring conversion to difference equations. Second, the use of difference equations makes the balance for the cell dynamics more explicit. With the assumptions given above, the dynamics of neoblasts are (N(t), N(

t), Q(t))
− p 3 f D (D(t))f N (N(t))a 3 (N(t), N(t), Q(t))] For those who prefer differential equations, one can proceed as follows. First replace N(t + 1) on the left hand side by N(t + t) where t is a suitably small unit of time. Second, define r i through the relationship p i = r i t + o( t) where o( t) denotes terms that are higher order in t. Then we have Thus, our difference equations are an Euler-method for the solution of the differential equations.

The production of progeny
As described above, we envision an intermediate progenitor cell between neoblasts and fully differentiated cells. This progenitor is not mitotically active and may continue development to a fully differentiated cell or may return to the resource pool. This is an inefficient process, but important if food finding is stochastic and lineage commitment must be made before food is searched for. In light of the feedback functions described above, the total production of progenitor cells at time t is We assume that a fraction 1 − f Q (Q(t)) of these progeny are returned to the resource pool and that the remaining fraction complete differentiation. For computations we assume f Q (q) = 1−exp(−β Q q), so that f Q (q) ≈ β Q q when q is small and f Q (q) → 1 as q increases. For cases in which we assume sufficient resources, we assume f Q (q) = 1.

Dynamics of differentiated cells
To model the dynamics of differentiated cells, we must capture two processes: the mortality of differentiated cells and the distribution of progenitors across the diversity of differentiated cells.
We let D i denote the number of differentiated cells when the planarian is in homeostasis, which we assume is set by natural selection and thus exogenous to the dynamics of cells within the life of a planarian. Then ρ i = D i where μ Q (n, d, q) is the resource-dependent rate of mortality, which we determine as follows. We set q μ = 1.5q 4 and assume that if q > 10q μ then μ Q (n, d, q) = 0 and otherwise In the second term of Eq. 9 μ i , σ D i and σ ρ i are fixed parameters. Note that when the planarian is in homeostasis with sufficient resources (μ Q (n, d, q) = 0), so that d i = D i and ρ i = ρ i , the second term on the right hand side of Eq. 9 is μ i .
To determine the allocation of progenitors across the different kinds of differentiated cells, we follow the relative need assumption as described above: at time t for the need of the i th kind of differentiated cell is D i D i (t) so that the relative need for the i th kind of cell is . With these assumptions, the dynamics of differentiated cells are The dynamics of the resource pool The resource pool increases by acquisition of resources from the environment, from progenitors that are directly returned to the pool, and from differentiated cells that die. It decreases due to metabolism of neoblasts and differentiated cells and through cell divisions. We assume that resource gain from the external environment is Y e D 1 (t) δ where Y e is a metric of food availability in the external environment, D 1 (t) is the number of differentiated cells used for food gathering and δ < 1 is a parameter accounting for not all cells being able to accumulate resources from the external environment. Resources returned to the pool from a progenitor that dies are γ p and from a neoblast or differentiated cell that dies is γ . Since m r and m d denote the resource cost of metabolism and division, the dynamics of the resource pool are

The steady state under sufficient resources
Under sufficient resources, with overline denoting the steady state value of a dynamical variable, we assume f N (N) = 1, a i (N, D, Q) = 1 for all i, μ Q (N, D, Q) = 0, and f Q (Q) = 1. Thus none of the neoblast divisions are resource limited and all of the progenitors continue to full differentiation, rather than returning to the resource pool.
The neoblast dynamics (Eq. 6) become We assume that at the steady state all the feedback functions in Eq. 2 take the same value, so that for each i with φ to be determined. Using Eq. 14 in Eq. 13 we have Note that this equation only makes sense if p 1 = p 3 φ. However, in light of Eq. 14, φ < 1, so we conclude p 1 < p 3 as a condition for the steady state. In this steady state, the production of non-mitotic progeny is, from Eq. 8, P = N · φ p 2 + 2p 3 (16) and Eq. 11 becomes Using Eq. 16 in Eq. 17 we obtain so that We now define which allows us to write D i = w i N. Consequently the fraction of neoblasts in the steady state is Thus we predict the same proportion of neoblasts in a planarian at the steady state, regardless of the number of cells and that this proportion is determined by the transition rates and rate of death of differentiated cells. In addition, the μ i will determine the relative abundance of differentiated cells in the steady state (and dynamically changing animal as well). For computations, we set μ i = s i μ 0 , where μ 0 is the baseline rate of mortality for differentiated cells and s i is a modulator according to the kind of differentiated cell.
In the steady state, Eq. 12 becomes Since D i = w i N this becomes an equation for the number of neoblasts which we solve to obtain Once this equation is solved, we compute the number of differentiated cells from D i = w i N and then determine the α i by solving Since the steady state number of neoblasts, and thus of differentiated cells depends upon the level of food in the environment, the feedback control parameter α i implicitly depends upon food, but this is emergent from the model, not an explicit assumption.

Size-cell number relationship
We compute the size of the planarian using the data from Table 1 of [15], which reports cell numbers and size for Dugesia mediterranea of 4, 7, 11, and 16 mm. We assume that size is determined only by the number of differentiated cells. Using those data, size S(t) at time t when the number of differentiated cells of type i is D i (t) is .99993 for the log-log plot).

Growth without resource constraints
If we assume that there are sufficient resources to support all divisions, then the full dynamics simplify considerably. That is, if we set a i (N(t), D(t), Q(t)) = 1, μ Q (N(t), D(t), Q(t) = 0, f Q (Q(t)) = 1, and set the resource dependent mortality in Eq. 9 equal to 0, then In this case, the dynamics of the resource pool are irrelevant. During such growth, the number of differentiated cells of type i dying in an interval of time will be D i (1 − exp(−M i (N(t), D i (t), Q(t)) so that the per capita death rate of these differentiated cells is 1 − exp(−M i (N(t), D i (t), Q(t)). This per-capita mortality is partitioned between the two types of non-homeostatic distributions of cells according to the relative values of the terms in Eq. 9. Our formulation allows us to separate mortality due to unbalance from the steady state level and mortality due to unbalance from the steady state proportions.

A simplified version of the model
For some questions, a version of the model with only one kind of differentiated cell suffices; a similar approach using only one kind of differentiated cell to investigate control of proliferation is taken in [45]. We briefly describe that version now. We replace the mortality function in Eq. 9 by μ Q (n, d, q) + μ where μ Q (n, d, q) is interpreted as before and μ is a baseline rate of mortality. With this assumption the dynamics of neoblasts and the single type of differentiated cell type become and and the dynamics of the resource pool are To compute the steady state without resource constraints, we set all the a i (n, d, q) = 1 and f N (n) = 1 and μ Q (n, d, q) = 0 (i.e., there are sufficient resources, neoblast transitions happen at their maximum possible value modified only by the feedback control from differentiated cells, and there are sufficient resources that resource dependent cell death does not occur). Following the procedure as above for the full model, we find as before p 1 so that if we set η = (1−e −μ )p 3 p 1 [p 2 +2p 3 ] the fraction of neoblasts in this steady state is η η+1 . In this simplified model, the steady state number of differentiated cells is determined from

Parameters
Although every parameter that is used in these models can be measured, most of them have not at this time. Indeed, one role of a paper such as this is to motivate empiricists to measure the parameters. We now explain how we determined the parameters. In general, we focus on the full model. When the simpler one differs from the full model, we explain the difference. Fundamental transition rates As described above, these parameters are connected to the physical interpretation of one unit of time in the model and are otherwise unconstrained except that p 1 < p 3 to ensure that a steady state exists. For computations here, we set p 1 = 0.0001, p 2 = 0.0005 and p 3 = 0.00015.
Food gathering and metabolic rates We assume that the exponent δ in Eq. 12 is described by the classic relationship between a linear variable and surface area, i.e. δ = 2/3 (qualitatively similar results are obtained with other choices, such as δ = 0.75). We choose metabolic rates in units so that the metabolic rate of a neoblast or differentiated cell is m r = 1.0 and assume that the cost of division is 4 times that, i.e. m d = 4.0. We assume that when a cell apoptoses 80 % of its resources return to the resource pool so that γ = 0.8m d . In this framework, we understand food in the environment, Y e in Eq. 12, to be multiples of m r .
Rates of cell death We assume that the μ i in Eq. 9 are multipliers of a basic mortality rate μ 0 so that μ i = s i μ 0 . We choose μ 0 = .00015 and s 1 = 0.75, s 2 = 1.0 and s 3 = 1.0. In the absence of resource constraints or deviations from homeostasis, the expected cell lifetime predicted from Eq. 9 is 1 μ i , which is another way of setting the link between one unit of time in the model and chronological time. We set σ D i and σ ρ i in Eq. 9 equal to 15 % of the steady values of D i and ρ i . For the results reported in this paper, we assume sufficiently large resources so that μ Q (n, d, q) = 0.
Feedback controls As described above, the parameters α i (or α in the simplified model) in Eq. 2 emerge from the steady state analysis. The choice of the functional form is somewhat arbitrary: we require that f D (d) declines as d increases and approaches 1 as d approaches 0 and the simple nonlinear form of Eq. 2 captures this idea without the risk of becoming negative as a linear function would; Taylor expanding these functions when α i D i << 1 gives 1 − α i D i (cf [46]). Similarly, the choice of q ij and the exponent 2 used in the feedback control functions a i (n, d, q) in Eqs. 1, 3, and 5 are arbitrary but capture the properties that we expect of such feedback functions. Finally, the feedback control f N (n) in Eq. 4 involves two parameters. We set the number of neoblasts at which the feedback control is 0.5 to N c = 0.15N(Y e ) and the parameter characterizing the spread of this function σ N = 0.15N(Y e ). These functions and parameters await experimental measurement. The analysis reported above shows the importance of combinations of parameters, perhaps even more than their individual values (cf [47]).

In silico experiments
With the full model, we do the following. First, we compute the steady state size as a function of food level. From that we compute the strength of feedback control. We then compute the size, total mortality of differentiated cells, and fraction of neoblasts under a temporal pattern of food in which food is dropped and then subsequently increased but there are sufficient resources for all the a i (n, d, q) = 1 and μ Q (n, d, q) = 0. Fourth, we follow the dynamics of cells during remodeling following a division, again under the assumption of sufficient resources. To do, this we assume that the initial cell numbers are 50 % neoblasts, and 10 %, 30 %, and 5 % of the three kinds of differentiated cells respectively (rather than 25 % neoblasts and 40 %, 30 %, 30 % relative distribution of differentiated cells in the steady state).
Using the simplified model, we do the following. First, we follow the dynamics during growth and regrowth. Second, we consider an excision experiment: we grow a planarian and then at scaled time t = 7 the number of differentiated cells is reduced by 25 %. For the x-ray experiment, we assume that at scaled time t = 7 the number of neoblasts is reduced by 25 %.