Skip to main content
  • Research Article
  • Open access
  • Published:

Feedback control in planarian stem cell systems

Abstract

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.

Background

Stem cell systems operate by demand control [13] 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 planarian flatworms (Tricladdida), particularly Dugesia and Phagocata species, have been key models in the study of regeneration and wound healing for more than 100 years (see [46] for some classic studies; [711] for more recent ones). Their simplicity and the ease with which regeneration experiments can be performed make them an attractive system for understanding the fundamental mechanisms of regeneration. Recent advances in molecular techniques have allowed deeper understanding of these apparently simple organisms; it is now possible to study the stem cell system and its descendants from the 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 [1618]. There are however populations of transient post-mitotic 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, non-mitotic 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 modelsIn 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. 1 a 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 NN, N; NN, P; and NP, P. The transitions for progeny cells are PQ and PD 1,PD 2, or PD 3.

Fig. 1
figure 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 (NN, N), symmetric differentiation (NN, P), and asymmetric differentiation (NP, 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 NP, 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 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 the transitions NN,N,NN, P, and NP, 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

These transitions have a maximum rate that is modified by feedback control. In Fig. 1 b, we show the positive feedback control on these transitions. Neoblasts exert positive feedback f N (N) on the NP, 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 NN, N,NN, P and NP, P respectively. In Fig. 1 c, 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. 2, we show examples of the feedback functions and the resource dependent rate of apoptosis. The feedback function f D (D) (Fig. 2 a) from differentiated cells on the transitions of neoblasts (Fig. 1 c) falls from 1.0 as the number of differentiated cells increases. This function has a parameter α that controls the rate of decline, in the sense that larger values of α mean smaller values of f D (D) for the same level of differentiated cells. The feedback control from neoblasts to asymmetric differentiation (Fig. 1 a, Fig. 2 b) 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. 2 c, 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. 2 d we show the resource dependent rate of natural mortality, again in cross-section. Full details are in the Section “Methods”.

Fig. 2
figure 2

a The feedback function f D (D) from differentiated cells on the transitions of neoblasts (Fig. 1 c) falls from 1.0 as the number of differentiated cells (here represented as the fraction of the steady state value) increases. b The feedback control from neoblasts to asymmetric differentiation (Fig. 1 b) increases towards 1.0 as the number of neoblasts (here represented as the fraction of the steady state value) increases. c The feedback functions a i (N,D,Q) depend, for a model with three kinds of differentiated cells, on five variables. Here we show a cross-section of those functions by holding the number of neoblasts and differentiated cells constant and only varying the resource level. d Similarly, the resource dependent rate of mortality of differentiated cells depends on the number of neoblasts and differentiated cells through their metabolic requirements; here we show a cross-section holding those cell numbers constant and varying resource levels. The x-axis intentionally has no units since these images are intended to be schematic rather than accurate representations. 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”.

Results

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).

Fig. 3
figure 3

The steady state planarian size for the full model as a function of food in the environment (based on Eqs. 2126)

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.

Fig. 4
figure 4

The strength of feedback control (α i ) from differentiated cells to neoblasts emerges as food-dependent from Eq. 2. The feedback functions α 2 and α 3 sit on top of each other in the upper curve and α 1 is the lower curve

In Eqs. 1321, 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 71, 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.

Fig. 5
figure 5

The temporal pattern of size, total mortality of differentiated cells, total mortality of differentiated cells, and fraction of neoblasts under growth, degrowth, and regrowth for a situation in which the planarian starts in an environment with Y e =450 at 20 % of its steady state size. For simplicity of presentation, we scale time and divide all times in the model by 10,000 for presentation in figures. At scaled time t= 7 food is dropped to Y e =400 and then at scaled time t=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 the activity of neoblasts

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”).

Fig. 6
figure 6

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

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. 7 a. 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. 7 b). During this entire process, however, the fraction of neoblasts is nearly constant (Fig. 7 c). 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. 7 d and e respectively) respond to the food pattern in complex ways, in part because α in the feedback control \(f_{D}(D)=\frac {1}{1+\alpha D}\) is a function of external food level (Fig. 4) and because of additional mortality when resources are insufficient.

Fig. 7
figure 7

Pattern of growth and shrinkage using the simplified model. a The temporal pattern of food. b Predicted size of the planarian as a function time. c The fraction of neoblasts as a function of time. d The feedback control from differentiated cells to neoblast activity. e The feedback control from neoblasts to themselves on the NP, P transitions

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. 8 b), an increase in the feedback control function f D (D) (Fig. 8 c) and a slight drop in f N (N) (Fig. 8 d, 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”).

Fig. 8
figure 8

The in silico excision experiment. a After the planarian has reached its steady state size at scaled time t=7 differentiated cells are removed by excision, leading to a drop in size. b This leads to an increase in the fraction of neoblasts c, an increase in the feedback control function f D (D)d and a slight drop in f N (N) (lower panel), note the vertical scale

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.

Fig. 9
figure 9

The in silico x-ray experiment. After the planarian has reached its steady state size at scaled time t=7 neoblasts are removed. This leads to a decrease in the fraction of neoblasts (upper panel), an increase in the feedback control function f D (D) (middle panel), and a significant drop but then recovery of f N (N) (lower panel)

Discussion

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 semi-quantitative; 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 NN,P or NP,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, 3235]; most recently this has included study of single cells [13] and comparisons between wildtype animals and RNAi loss of function phenotypes [3639]. 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. 1.

    a curvilinear relationship between external food and planarian steady state size;

  2. 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. 3.

    the fraction of neoblasts in the steady state is constant regardless of planarian size;

  4. 4.

    planarians adjust size when food shifts first due to apoptosis and then through a reduction in neoblast activity;

  5. 5.

    a burst of mortality during regeneration as the number of differentiated cells are adjusted towards their homeostatic level;

  6. 6.

    following wounding or excision of differentiated cells, different time scales characterize the recovery of size and the two feedback functions;

  7. 7.

    the temporal pattern of feedback controls differs noticeably during recovery from a removal of neoblasts or a removal of differentiated cells; and

  8. 8.

    the signalling for apoptosis of differentiated cells depends upon both the absolute and relative deviations of differentiated cells from their homeostatic level.

  9. 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 t, D i (t) the number of differentiated cells of type i at time t, for i=1,2,…I (in computations we set I=3), and D(t)=(D 1(t),D 2(t),…D I (t)) the entire collection of differentiated cells.

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. NN, N; 2) symmetric renewal and progeny production, i.e. NN, P ; and 3) asymmetric progeny production, i.e. NP, 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 NN, N,NN, P, and NP, P transitions if they are sufficient. Thus, we assume that the fraction of neoblasts undergoing the NN, 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+\Sigma _{i=1}^{I} 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 NN, N divisions by a Hill-type function

$$ a_{1}(n,\textbf{d},q)=\frac{(q-q_{1})^{2}}{(q-q_{1})^{2} +(q_{12}-q_{1})^{2}} $$
((1))

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_{12}^{2}\) 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 NN+P and NP,P involve resource-dependent and cell-number dependent feedback control. We assume that the fraction of neoblasts undergoing NN+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

$$ f_{D}(\textbf{d})=\max_{\substack{i}}\bigg[ \frac{1}{1+ \alpha_{i} d_{i}} \bigg] $$
((2))

where α i sets the strength of control from differentiated cells of type i to neoblasts. The “ maxi” 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 NN+P transitions. We let q 23 denote the midpoint of q 2 and q 3 and set

$$ a_{2}(n,\textbf{d},q)=\frac{(q-q_{2})^{2}}{(q-q_{2})^{2} +(q_{23}-q_{2})^{2}} $$
((3))

as long as q>q 2; otherwise we set a 2(n,d,q)=0.

We assume that NP, P involves an additional feedback control from neoblasts; so that when neoblast numbers are low this transition is suppressed. We set

$$ f_{N}(n)=\frac{exp\bigg(\frac{n-N_{c}}{\sigma_{N}} \bigg)}{1+exp\bigg(\frac{n-N_{c}}{\sigma_{N}} \bigg)} $$
((4))

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 NP, 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

$$ a_{3}(n,\textbf{d},q)=\frac{(q-q_{3})^{2}}{(q-q_{3})^{2} +(q_{34}-q_{3})^{2}} $$
((5))

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/10000th 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

$$ \begin{aligned} N(t+1)&=N(t)[1+ p_{1} a_{1}(N(t),\textbf{N}(t),Q(t))\\ &\quad-p_{3} f_{D}(\mathbf{D}(t)) f_{N}(N(t)) a_{3}(N(t),\textbf{N}(t),Q(t))] \end{aligned} $$
((6))

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

$$ \begin{aligned} \frac{dN}{dt}&=r_{1} a_{1}(N(t),\textbf{N}(t),Q(t))N(t)\\ &\quad-r_{3} f_{D}(\mathbf{D}(t)) f_{N}(N(t)) a_{3}(N(t),\textbf{N}(t),Q(t))N(t) \end{aligned} $$
((7))

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

$$ \begin{aligned} \mathcal P(t)&= N(t) f_{D}(\mathbf{D}(t))\left[ p_{2} a_{2}(N(t),\mathbf{D}(t),Q(t))\right.\\& \quad+\left. 2 p_{3} f_{N}(N(t))a_{3}(N(t),\mathbf{D}(t),Q(t))\right] \end{aligned} $$
((8))

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−e x p(−β 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 \(\overline {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 \(\overline {\rho }_{i}=\frac {\overline {D}_{i}}{\Sigma _{k=1}^{I} \overline {D}_{k}}\) is the fraction of differentiated cells of type i in the steady state. Similarly, we let \(\rho _{i}(t)=\frac {D_{i}(t)}{\Sigma _{k=1}^{I} D_{k}(t)}\) represent the fraction of differentiated cells of type i at time t. We assume that the rate of mortality of differentiated cells depends upon i) how far D i (t) is from \(\overline {D}_{i}(t)\), ii) how far ρ i (t) is from \(\overline {\rho }_{i}(t)\), and iii) whether there are sufficient resources to maintain existing neoblasts and differentiated cells. In particular, we assume that the rate of mortality of differentiated cells of type i, given N(t)=n,D(t)=d,Q(t)=q is

$${} \begin{aligned} M_{i}(n,\textbf{d},q)&=\mu_{Q}(n,\textbf{d},q)\\ & \quad+\mu_{i}\left[\! \frac{exp\left(\frac{d_{i}-\overline{D}_{i}}{\sigma_{D_{i}}}\right)}{1+exp\left(\frac{d_{i}-\overline{D}_{i}}{\sigma_{D_{i}}}\right)} + \frac{exp\left(\frac{\rho_{i}-\overline{\rho}_{i}}{\sigma_{\rho_{i}}}\right)}{1+exp\left(\frac{\rho_{i}-\overline{\rho}_{i}}{\sigma_{\rho_{i}}}\right)} \!\right] \end{aligned} $$
((9))

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

$$ \mu_{Q}(n,\textbf{d},q)=\mu_{Q,max}\cdot \frac{(q-q_{\mu})^{2}}{(q-q_{\mu})^{2}+(q_{1}-q_{\mu})^{2}} $$
((10))

In the second term of Eq. 9 \(\phantom {\dot {i}\!}\mu _{i}, \sigma _{D_{i}}\) and \(\sigma _{\rho _{i}}\) are fixed parameters. Note that when the planarian is in homeostasis with sufficient resources (μ Q (n,d,q)=0), so that \(d_{i}=\overline {D}_{i}\) and \(\rho _{i}=\overline {\rho }_{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 \( \frac {\overline {D}_{i}}{D_{i}(t)}\) so that the relative need for the i th kind of cell is \(\frac {\frac {\overline {D}_{i}}{D_{i}(t)}}{\Sigma _{k=1}^{I} \frac {\overline {D}_{k}}{D_{k}(t)}} \). With these assumptions, the dynamics of differentiated cells are

$$ \begin{aligned} D_{i}(t+1)&=D_{i}(t)exp(-M_{i}(N(t),\textbf{D(t)},Q(t)))\\ \quad &\quad +f_{Q}(Q(t) \mathcal P(t) \frac{\frac{\overline{D}_{i}}{D_{i}(t)}}{\Sigma_{k=1}^{I} \frac{\overline{D}_{k}(t)}{D_{k}}} \end{aligned} $$
((11))

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

$$ \begin{aligned} Q(t+1) &= Q(t) + Y_{e} D_{1}(t)^{\delta} + \gamma_{p} (1-f_{Q}(Q(t))\mathcal P(t)\\ &\quad+\gamma \Sigma_{i=1}^{I} D_{i}(t)(1\,-\,exp(-M_{i}(N(t),\textbf{D(t)},Q(t)))\\ &\quad-m_{r}\left[N(t) +\Sigma_{i=1}^{I} D_{i}(t)\right]\\ &\quad-m_{d} N(t)\left[p_{1} a_{1}(N(t),\mathbf{D}(t),Q(t))\right.\\ &\quad+\left. f_{D}(\mathbf{D}(t))\left[p_{2} a_{2}(N(t),\mathbf{D}(t),Q(t))\right.\right.\\ &\quad +\left.\left. f_{N}(N(t))p_{3} a_{3}(N(t),\mathbf{D}(t),Q(t))\right]\right] \end{aligned} $$
((12))

The steady state under sufficient resources

Under sufficient resources, with overline denoting the steady state value of a dynamical variable, we assume \(f_{N}(\overline {N})=1,a_{i}(\overline {N},\overline {\mathbf {D}},\overline {Q})=1\) for all i, \(\mu _{Q}(\overline {N},\overline {\mathbf {D}},\overline {Q})=0\), and \(f_{Q}(\overline {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

$$ \overline{N}=\overline{N}\left[1+p_{1}-p_{3} f_{D}(\overline{\mathbf{D}})\right] $$
((13))

We assume that at the steady state all the feedback functions in Eq. 2 take the same value, so that for each i

$$ \bigg[ \frac{1}{1+ \alpha_{i} \overline{D}_{i}} \bigg]=\overline{\phi} $$
((14))

with \(\overline {\phi }\) to be determined. Using Eq. 14 in Eq. 13 we have

$$ \overline{N}=\overline{N}\left[1+p_{1}-p_{3}\overline{\phi}\right] $$
((15))

Note that this equation only makes sense if \(p_{1}=p_{3}\overline {\phi }\). However, in light of Eq. 14, \(\overline {\phi }<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,

$$ \overline{\mathcal P}=\overline{N} \cdot \overline{\phi}\left[p_{2} +2p_{3}\right] $$
((16))

and Eq. 11 becomes

$$ \overline{D}_{i}=\overline{D}_{i} exp[-\mu_{i}] +\frac{\overline{\mathcal P}}{I} $$
((17))

Using Eq. 16 in Eq. 17 we obtain

$$ \overline{D}_{i}=\overline{D}_{i} exp\left[-\mu_{i}\right] +\frac{\overline{N} \cdot \overline{\phi}\left[p_{2} +2p_{3}\right]}{I} $$
((18))

so that

$$ \overline{D}_{i}=\frac{\overline{N} \cdot \overline{\phi}\left[p_{2} +2p_{3}\right]}{I(1-exp\left[-\mu_{i}\right])} $$
((19))

We now define

$$ w_{i}=\frac{\overline{\phi}\left[p_{2} +2p_{3}\right]}{I(1-exp\left[-\mu_{i}\right])} $$
((20))

which allows us to write \(\overline {D}_{i}=w_{i}\overline {N}\). Consequently the fraction of neoblasts in the steady state is

$$ \overline{N}=\frac{N}{N+\Sigma_{i=1}^{I} \overline{D}_{i}}=\frac{1}{1+\Sigma_{i=1}^{I} w_{i}} $$
((21))

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

$$ \begin{aligned} \overline{Q}&=\overline{Q} + Y_{e} \overline{D}_{1}^{\delta}+\gamma \Sigma_{i=1}^{I} \overline{D}_{i} (1- exp[-\mu_{i}])\\&\quad-m_{r}\left[\overline{N} +\Sigma_{i=1}^{I} \overline{D}_{i}\right]-m_{d}\overline{N}\left(p_{1}+\overline{\phi}\left[p_{2}+p_{3}\right]\right) \end{aligned} $$
((22))

Since \(\overline {D}_{i}=w_{i}\overline {N}\) this becomes an equation for the number of neoblasts

$$ \begin{aligned} 0&=Y_{e} w_{1}^{\delta}\overline{N}^{\delta} +\gamma \overline{N} \Sigma_{i=1}^{I} w_{i} (1- exp\left[-\mu_{i}\right])\\&\quad-m_{r}\overline{N}\left[1 +\Sigma_{i=1}^{I} w_{i}\right]-m_{d}\overline{N}\left(p_{1}+\overline{\phi}\left[p_{2}+p_{3}\right]\right) \end{aligned} $$
((23))

which we solve to obtain

$${} {\fontsize{7.4}{8}{\begin{aligned} {}\overline{N}^{1-\delta}&\,=\,\frac{Y_{e} w_{1}^{\delta}}{m_{r}\left[1+\Sigma_{i=1}^{I} w_{i}\right]+m_{d}\left(p_{1}+\overline{\phi}\left[p_{2}+p_{3}\right]-\gamma \Sigma_{i=1}^{I} w_{i}(1-exp[\!-\mu_{i}]\right)} \end{aligned}}} $$
((24))

Once this equation is solved, we compute the number of differentiated cells from \(\overline {D}_{i}=w_{i} \overline {N}\) and then determine the α i by solving

$$ \overline{\phi}=\frac{1}{1+\alpha_{i} \overline{D}_{i}} $$
((25))

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

$$ log(S(t))=-3.8373+0.46582\cdot log\left(\Sigma_{i=1}^{I} D_{i}(t)\right) $$
((26))

(R 2 =0.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

$${} N(t+1)=N(t)\left[1+ p_{1} -p_{3} f_{D}(\mathbf{D}(t)) f_{N}(N(t))\right]\\ $$
$${} \begin{aligned} {} D_{i}(t+1)&=D_{i}(t)exp(-M_{i}(N(t),\textbf{D(t)},Q(t)))\\&\quad+ N(t) f_{D}(\mathbf{D}(t))\left[p_{2} + 2 p_{3} f_{N}(N(t))\right]\! \frac{\frac{\overline{D}_{i}}{D_{i}(t)}}{\Sigma_{k=1}^{I} \frac{\overline{D}_{k}}{D_{k}(t)}} \end{aligned} $$
((27))

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−e x p(−M i (N(t),D i (t),Q(t)) so that the per capita death rate of these differentiated cells is 1−e x p(−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

$${} \begin{aligned} N(t+1)&=N(t)\left[1+p_{1}a_{1}(N(t),D(t),Q(t))\right.\\ &\quad-\left. p_{3}f_{D}(D(t))f_{N}(N(t))a_{3}(N(t),D(t),Q(t))\right] \end{aligned} $$
((28))

and

$$\begin{array}{*{20}l} D(t+1) &= D(t)e^{-\mu-\mu_{Q}(N(t),D(t),Q(t))}\\ &\quad +N(t)f_{D}(D(t))\left[p_{2}a_{2}(N(t),D(t),Q(t))\right.\\ &\quad +\left.2p_{3}f_{N}(N(t))a_{3}(N(t),D(t),Q(t))\right] \end{array} $$
((29))

and the dynamics of the resource pool are

$$ \begin{aligned} Q(t+1)&=Q(t)+Y_{e} D(t)^{\delta}\\ &\quad + \gamma D(t)\left(1-e^{-\mu-\mu_{Q}(N(t),D(t),Q(t))}\right)\\ &\quad-m_{r}\left[N(t)+D(t)\right]\\ &\quad-m_{d} N(t)\left[p_{1}a_{1}(N(t),D(t),Q(t))\right.\\ &\quad+f_{D}(D(t))\left[p_{2} a_{2}(N(t),D(t),Q(t))\right.\\ &\quad+\left.\left.f_{N}(N(t))p_{3}a_{3}(N(t),D(t),Q(t))\right]\right] \end{aligned} $$
((30))

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 \(\frac {p_{1}}{p_{3}}=f_{D}(\overline {D})\) from which we conclude \(\alpha =\frac {1}{\overline {D}}\left (\frac {p_{3}}{p_{1}}-1\right)\) and

$$ \overline{N}=\overline{D}\frac{\left(1-e^{-\mu}\right)p_{3}}{p_{1}\left[p_{2}+2p_{3}\right]} $$
((31))

so that if we set \(\eta =\frac {(1-e^{-\mu })p_{3}}{p_{1}[p_{2}+2p_{3}]}\) the fraction of neoblasts in this steady state is \(\frac {\eta }{\eta +1}\).

In this simplified model, the steady state number of differentiated cells is determined from

$${} {\fontsize{8.9}{8}{\begin{aligned} \overline{D}^{1-\delta}=\frac{Y_{e}}{m_{r}\left[\eta+1\right]+m_{d}\eta\left[p_{1}+\frac{p_{1}}{p_{3}}\left[p_{2}+p_{3}\right]\right]-\gamma(1-e^{-\mu})} \end{aligned}}} $$
((32))

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 \(\frac {1}{\mu _{i}}\), which is another way of setting the link between one unit of time in the model and chronological time. We set \(\sigma _{D_{i}}\) and \(\sigma _{\rho _{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.15\overline {N}(Y_{e})\) and the parameter characterizing the spread of this function \(\sigma _{N}=0.15 \overline {N}(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 %.

Availability of data and materials

Code (written in C for computation and R for visualization) is included as a Additional file 1.

Endnote

1 As described below, for purposes of presentation we use a scaled time on the x-axis.

References

  1. Metcalf D. The Molecular Control of Blood Cells. Cambridge: Harvard University Press; 1988.

    Google Scholar 

  2. de Graaf CA, Kauppi M, Baldwin T, Hyland CD, Metcalf D, Willson TA, et al.Regulation of hematopoietic stem cells by their mature progeny. Proc Nat Acad Sci US. 2010; 107:21689–94.

    Article  CAS  Google Scholar 

  3. Mangel M, Bonsall MB. Stem cell biology is population biology: differentiation of hematopoietic multipotent progenitors to common lymphoid and myeloid progenitors. Theor Biol Med Modell. 2013; 10(5). [doi:10.1186/1742-4682-10-5].

  4. Morgan TH. Regeneration. London: MacMillan Company; 1901 Lemon CC: notes on the physiology of regeneration of parts in Planaria Maculata. Biol Bull. 1900; 1:193–204.

    Article  Google Scholar 

  5. Child CM. Senescence and Rejuvenescence. Chicago: University of Chicago Press; 1915.

    Book  Google Scholar 

  6. Brondsted HV. Planarian Regeneration. Oxford: Pergamon Press; 1969.

    Google Scholar 

  7. Aboobaker AA. Planarian stem cells: a simple paradigm for regeneration. Trends Cell Biol. 2011; 21:304–11.

    Article  CAS  PubMed  Google Scholar 

  8. Newmark PA, Sanchez Alvarado A. Not your father’s planarian: a classic model enters the era of functional genomics. Nat Rev Gen. 2002; 3:210–9.

    Article  CAS  Google Scholar 

  9. Adler CE, Sánchez Alvarado A. Types or states? Cellular dynamics and regenerative potential. Trends Cell Biol. 2015. doi:10.1016/j.tcb.2015.07.008;2015. Accessed 19 Dec 2015.

  10. Salo E, Agata K. Planarian regeneration: a classic topic claiming new attention. Int J Dev Biol. 2012; 56:1–4.

    Article  Google Scholar 

  11. De Mulder K, Pfister D, Kuales G, Egger B, Salvenmoser W, Willems M, et al. Stem cells are differentially regulated during development, regeneration and homeostasis in flatworms. Dev Biol. 2009; 334:198–212.

    Article  CAS  PubMed  Google Scholar 

  12. Wagner DE, Wang IE. Reddien PW.Clonogenic neoblasts are pluripotent adult stem cells that underlie planarian regeneration. Science. 2011; 332:811–6.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  13. van Wolfswinkel JC, Wagner DE, Reddien PW. Single-cell analysis reveals functionally distinct classes within the planarian stem cell compartment. Cell Stem Cell. 2014; 15:326–39.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  14. Aboukhatwa E, Aboobaker A. An Introduction to Planarians and Their Stem Cells. eLS; 2015. http://onlinelibrary.wiley.com/doi/10.1002/9780470015902.a0001097.pub2/abstract?userIsAuthenticated=false&deniedAccessCustomisedMessage=. Accessed 8 Feb 2015.

  15. Baguna J, Romero R. Quantitative analysis of cell types during growth, regrowth and regeneration in the planarians Dugesia meidterranea and Dugesia tigrina. Hydrobiologia. 1981; 84:181–94.

    Article  Google Scholar 

  16. Rink JC. Stem cell systems and regeneration in planaria. Dev Genes Evol. 2012. [doi:10.1007/s00427-012-0426-4].

  17. Almuedo-Castillo M, Crespo X, Seebeck F, Bartscherer K, Salo‘ E, Adell T. JNK controls the onset of mitosis in planarian stem cells and triggers apoptotic cell death required for regeneration and remodeling. PLoS Genet. 2014; 10(6):e1004400.

    Article  PubMed Central  PubMed  Google Scholar 

  18. Salo E, Abril JF, Adell T, Cebria F, Eckelt K, Fernandez-Taboada E, et al.Planarian regeneration: achievements and future directions after 20 years of research. Int J Dev Biol. 2009; 53:1317–27.

    Article  PubMed  Google Scholar 

  19. Oviedo NJ, Newmark PA, Sanchez Alvarado A. Allometric scaling and proportion regulation in the freshwater planarian Schmidtea mediterranea. Dev Dynam. 2003; 226:326–33.

    Article  CAS  Google Scholar 

  20. Pellettieri J, Sanchez Alvarado A. Cell turnover and adult tissue homeostasis: from humans to planarians. Annu Rev Genet. 2007; 41:83–105.

    Article  CAS  PubMed  Google Scholar 

  21. Pellettieri J, Fitzgerald P, Watanabe S, JMancuso J, Green DR, Sánchez Alvarado A. Cell death and tissue remodeling in planarian regeneration. Dev Biol. 2010; 338:76–85.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  22. Salo E. The power of regeneration and the stem-cell kingdom: freshwater planarians (Platyhelminthes). BioEssays. 2006; 28:546–59.

    Article  CAS  PubMed  Google Scholar 

  23. Reddien PW, Sánchez Alvarado A. Fundamentals of planarian regeneration. Ann Rev Cell Dev Biol. 2004; 20:725–57.

    Article  CAS  Google Scholar 

  24. Forsthoefel DJ, James NP, Escobar DJ, Stary JM, Vieira AP, Waters FA, et al.An RNAi screen reveals intestinal regulators of branching morphogenesis, differentiation, and stem cell proliferation in planarians. Dev Cell. 2012; 23(4):691–704.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  25. Lei J, Levin SA, Nie Q. Mathematical model of adult stem cell regeneration with cross-talk between genetic and epigenetic regulation. Proc Nat Acad Sci US. http://www.pnas.org/content/111/10/E880. Accessed 8 Feb 2016.

  26. Rørvik Høyem M, Måløy F, Jakobsen P, Brandsdal BO. Stem cell regulation: Implications when differentiated cells regulate symmetric stem cell division. J Theor Biol. 2015. [doi:10.1016/j.jtbi].

  27. Bardeen CR, Baetjer FH. The inhibitive action of the Roentgen rays on regeneration in planarians. J Exp Zool. 1904; 1:191–5.

    Article  Google Scholar 

  28. Takeda H, Nishimura K, Agata K. Planarians maintain a constant ratio of different cell types during changes in body size by using the stem cell system. Zool Sci. 2009; 26:805–13.

    Article  PubMed  Google Scholar 

  29. González-Estévez C, Felix DA, Rodrguez-Esteban G, Aboobaker AA. Decreased neoblast progeny and increased cell death during starvation-induced planarian degrowth. Int J Dev Biol. 2012; 56:83–91.

    Article  PubMed  Google Scholar 

  30. Wenemoser D, Redden PW. Planarian regeneration involves distinct stem cell responses to wounds and tissue absence. Dev Biol. 2010; 344:979–91.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  31. Umesono Y, Agata K. Evolution and regeneration of the planarian central nervous system. Develop Growth Differ. 2009; 51:185–95.

    Article  CAS  Google Scholar 

  32. Solana J, Kao D, Mihaylova Y, Jaber-Hijazi F, Malla S, Wilson R, et al.Defining the molecular profile of planarian pluripotent stem cells using a combinatorial RNAseq, RNA interference and irradiation approach. Genome Biol. 2012; 13(3):R19. [doi:10.1186/gb-2012-13-3-r19].

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  33. Blythe MJ, Kao D, Malla S, Rowsell J, Wilson R, Evans D, et al.A dual platform approach to transcript discovery for the planarian Schmidtea mediterranea to establish RNAseq for stem cell and regeneration biology. PLoS One. 2010; 5(12):e15617. [doi:10.1371/journal.pone.0.015617].

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  34. Sasidharan V, Lu YC, Bansal D, Dasari P, Poduval D, Seshasayee A, et al.Identification of neoblast- and regeneration-specific miRNAs in the planarian Schmidtea mediterranea. RNA. 2013; 19(10):1394–404. [doi:10.1261/rna.038653.113].

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  35. Labbé RM, Irimia M, Currie KW, Lin A, Zhu SJ, Brown DD, et al.A comparative transcriptomic analysis reveals conserved features of stem cell pluripotency in planarians and mammals. StemCells. 2012; 30(8):1734–45. [doi:10.1002/stem.1144].

    Google Scholar 

  36. Kao D, Felix D, Aboobaker A. The planarian regeneration transcriptome reveals a shared but temporally shifted regulatory program between opposing head and tail scenarios. BMC Genomics. 2013; 14:797. [doi:10.1186/1471-2164-14-797].

    Article  PubMed Central  PubMed  Google Scholar 

  37. Tu KC, Cheng LC, Tk Vu H, Lange JJ, McKinney SA, Seidel CW, et al.Egr-5 is a post-mitotic regulator of planarian epidermal differentiation. Elife. 2015; 4. pii: e10501. doi:10.7554/eLife.10501.

  38. Reuter H, März M, Vogg MC, Eccles D, Grífol-Boldú L, Wehner D, et al. β-catenin-dependent control of positional information along the AP body axis in planarians involves a teashirt family member. Cell Rep. 2015; 10(2):253–65. [doi:10.1016/j.celrep.2014.12.018].

    Article  CAS  PubMed  Google Scholar 

  39. Cowles MW, Omuro KC, Stanley BN, Quintanilla CG, Zayas RM. COE loss-of-function analysis reveals a genetic program underlying maintenance and regeneration of the nervous system in planarians. PLoS Genet. 2014; 10(10):e1004746. [doi:10.1371/journal.pgen.1004746].

    Article  PubMed Central  PubMed  Google Scholar 

  40. Salo E, Baunga J. Regeneration and pattern formation in planarians. Development. 1989; 107:69–76.

    Google Scholar 

  41. Wilkinson DJ. Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev. 2009; 10:122–33.

    Article  CAS  Google Scholar 

  42. Sun Z, Komarova NL. Stochastic control of proliferation and differentiation in stem cell dynamics. J Math Biol. 2014. [doi:10.1007/s00285-014-0835-2].

  43. Mangel M, Bonsall MB. Phenotypic evolutionary models in stem cell biology: replacement, quiescence, and variability. PLoS ONE. 2008; 3(2):e1591. [doi:10.1371/journal.pone.0001591].

    Article  PubMed Central  PubMed  Google Scholar 

  44. Gillespie DT, Mangel M.Conditioned averages in chemical kinetics. J Chem Phys. 1981; 75:704–9.

    Article  CAS  Google Scholar 

  45. Buzi G, Lander AD, Khammash M. Cell lineage branching as a strategy for proliferative control. BMC Biol. 2015; 13:13–28. [doi:10.1186/s12915-015-0122-8].

    Article  PubMed Central  PubMed  Google Scholar 

  46. MacLean AL, Kirk PDW, Stumpf MPH. Cellular population dynamics control the robustness of the stem cell niche. Biol Open. 2015; 4:1420–6.

    Article  PubMed Central  PubMed  Google Scholar 

  47. Goyal S, Kim S, Chen ISY, Chou T. Mechanisms of blood homeostasis: lineage tracking and a neutral model of cell populations in rhesus macaques. BMC Biol. 2015; 12:85–99.

    Article  Google Scholar 

Download references

Acknowledgments

MM was supported by NSF grant EF-0924195; AA by BBSRC Research Grant BB/K007564/1 and MRC Research Grant MR/M000133/1.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Marc Mangel.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

MM developed initial forms of the models, which were then modified through discussions with AA and MB. MM did the computations. All authors contributed to the manuscript, with MM leading the writing. All authors have read and approved the final version of the manuscript.

Additional file

Additional file 1

The supplementary materials contains the codes for both models. Four codes are included in one PDF document (so need to be identified and separated if used): C script for the full model; R script for visualizing output of the full model; C script for the simplified model; and R script for visualizing output of the simplified model. Computations were done on a 2012 MacBook Pro using the gcc compiler with home directory ’marco’. (PDF 323 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Mangel, M., Bonsall, M.B. & Aboobaker, A. Feedback control in planarian stem cell systems. BMC Syst Biol 10, 17 (2016). https://doi.org/10.1186/s12918-016-0261-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12918-016-0261-8

Keywords