A model of stem cell population dynamics: in silico analysis and in vivo validation

Many tissues and organs contain self-renewing stem cell populations that are crucial for their maintenance. Improper control of stem and progenitor cell populations may underlie tissue degeneration and cancer. Two general mechanisms for stem cell renewal by cell proliferation are: (1) strict asymmetric divisions, in which one cell retains stem cell fate while the other cell or its progeny differentiate; and (2) renewal at the population level, whereby the sum total of asymmetric and/or symmetric divisions (symmetric in the sense that a cell division yields two stem or two non-stem cells) maintains the population. Two additional, not necessarily mutually exclusive, mechanisms of renewal are dedifferentiation and renewal by facultative division of differentiated cells (see Fuller and Spradling, 2007; He et al., 2009; Morrison and Kimble, 2006; Puri and Hebrok, 2010; Yamashita et al., 2010). Often, stem cells are associated with an anatomically distinct niche. The niche communicates with stem cells via cell-cell signaling to direct stem cell fate, survival and proliferation potential (see Drummond-Barbosa, 2008; Voog and Jones, 2010).

Arduous cell-autonomy and lineage-tracing studies are required to rigorously define renewal mechanisms. Ideally, lineally labeled cells are followed with non-invasive long-term imaging in both wild-type and informative mutant cell contexts. However, even when lineage information is available, the influence of anatomical constraints is difficult to assess in conjunction with signaling and cell division. Although lineage tracing is not currently feasible for many tissue types, anatomical, molecular and cell cycle information from static samples is often available. We reasoned that dynamic computational models with the potential to systematically manipulate different influences might facilitate an understanding of experimental studies on self-renewing cell populations.

We chose the C. elegans germ line as a general test for the utility of a dynamic modeling approach. Although this system is genetically and anatomically tractable, definitive stem cell markers are not available and lineage analysis is not yet feasible. A population of proliferating stem and/or transit-amplifying cells (here referred to collectively as ‘progenitors’) maintains the germ line. These cells are located in the distal-most region (the proliferative zone) of a blind-ended tube (the gonad arm, two of which are present in the hermaphrodite; Fig. 1). Progenitors undergo cell divisions in larval and adult stages, reaching and maintaining a population of ∼200 progenitors in the adult gonad arm. Starting in the third larval stage (L3) and continuing through L4 and adult, germ cells differentiate and ultimately produce gametes. Maintenance of an undifferentiated proliferation-competent progenitor population depends on signaling from a single cell, termed the distal tip cell (DTC), that caps each end of the blind-ended tube. The DTC produces ligands for the Notch family receptor GLP-1, which is expressed in the germ line (Austin and Kimble, 1987; Crittenden et al., 1994; Henderson et al., 1994). Withdrawal of Notch pathway signaling causes all germ cells to enter meiosis and differentiate (Austin and Kimble, 1987; Lambie and Kimble, 1991), whereas hyperactive signaling causes uncontrolled proliferation (Berry et al., 1997; Pepper et al., 2003a; Pepper et al., 2003b) (Fig. 1; supplementary material Fig. S1). GLP-1-mediated signaling opposes the activities of redundant genetic pathways that lead to meiotic entry, two of which are defined by the GLD-1 and GLD-2 proteins (Hansen et al., 2004a). Additional, non-DTC signals influence the establishment of the adult proliferative zone, including insulin/IGF-like signaling (Michaelson et al., 2010), and signals from the gonadal sheath cells (Killian and Hubbard, 2005). The progenitor population has been divided into subzones based on cell behavior and the expression of a number of genes and proteins (Cinquin et al., 2010; Crittenden et al., 2002; Crittenden et al., 2006; Hansen et al., 2004a; Hansen et al., 2004b; Hubbard, 2007; Jaramillo-Lambert et al., 2007; Lamont et al., 2004; Maciejowski et al., 2006; Merritt and Seydoux, 2010).

Although this system is relatively simple, many central features remain unresolved or controversial. For example, because the field lacks a cell-specific readout for GLP-1/Notch receptor activity, the precise range of signaling is unknown. The DTC body itself directly contacts germ cells over a distance of ∼2-3 cell diameters (CD). However, this cell bears long processes that contact germ cells further from the tip, and the role of these processes in signaling, if any, is unclear. Some processes extend to the proliferation/differentiation border (∼20 CD) in the early adult (Fitzgerald and Greenwald, 1995; Henderson et al., 1994), but their average length is shorter [∼8-10 CD (Hall et al., 1999)] and does not always correlate with underlying germline differentiation (Crittenden et al., 2006). The field also currently lacks markers to trace lineage and to distinguish putative stem cells from transit-amplifying cells. Finally, long-term imaging in real time is technically challenging. Nevertheless, complex cell behaviors have been tracked by arduous large-scale timecourse analysis of fixed specimens (Cinquin et al., 2010; Crittenden et al., 2006; Hansen et al., 2004a; Jaramillo-Lambert et al., 2007; Maciejowski et al., 2006; Michaelson et al., 2010).

Using this system as inspiration, we developed a general dynamic computational model based on a subset of available data that represents an anatomically restricted cell population responding to a localized signal from a niche, but for which the cellular basis for renewal is unknown. We analyzed simulations from this model to determine the cell population level outcomes resulting from the combined effects of signaling, cell cycle and anatomical constraints on individual cells. The model produced simulations that recapitulate basic behaviors of the system, including proper responses to genetic manipulations. Our analyses of model simulations and laboratory results suggest that: (1) when the ligand interaction occurs over a short distance (that is, reaching only the distal-most germ cells), small differences in this distance destabilize the system and introduce unexpected variability; (2) inherent differences between progenitor cell types need not be invoked to explain previously observed complex differentiation dynamics upon reduction of receptor activity; (3) population dynamics and anatomical constraints influence niche residence; and (4) the germ cell proliferation rate during larval stages influences the differentiation pattern in the adult.

The model was built using Statecharts (Harel, 1987), a visual language suited for the specification of complex reactive systems. We used an object-oriented version of Statecharts, as implemented in Rhapsody (Harel and Gery, 1997; Harel and Kugler, 2004), that graphically represents dynamics of classes and instantiated objects using states, transitions, events and conditions. As in classical state machines, behavior is described using states and transitions between states, which are triggered by the occurrence of events. However, Statecharts extend classical state machines allowing concise representations by supporting hierarchal states that define detailed behavior within the parental state. Statecharts also supports orthogonal concurrent states, allowing the system to reside simultaneously in several different states (see supplementary material Appendix S1).

Using Rhapsody (IBM; Harel and Kugler, 2004) as a reactive engine, we used code generation and compilation (C++ for Microsoft Visual Studio 2008) to generate executable code that responds to events, changes object properties and creates new objects at run time. Using generic reactive animation (Efroni et al., 2005; Harel and Setty, 2008) we linked the Rhapsody model with a real-time Flash-based front-end (Flash Professional CS3, Adobe). Each object in the simulation is represented as an animated figure, which changes properties such as size and color to designate changes during the lifetime of the object. The front-end reflects changes in real time, providing an animated representation of the simulation over time and space. The general approach was similar to that employed by Setty et al. in a model of pancreas development (Setty et al., 2008).

N2 Bristol (Brenner, 1974), glp-1(e2141) (Priess et al., 1987), glp-1(oz112) (Berry et al., 1997), lag-2(q420) (Lambie and Kimble, 1991) and daf-2(e1370) (Riddle et al., 1981) strains were grown according to standard methods (Brenner, 1974). For glp-1(e2141) mutants, after hypochlorite treatment of gravid adults, eggs were washed twice in M9 buffer and incubated on a platform shaker at 15°C overnight. L1 larvae were washed and transferred to plates containing OP50 bacteria at 25°C, harvested 48 hours later (young adult stage) and imaged live. Adult glp-1(oz112) animals were imaged live. lag-2 mutants were handled as for glp-1(e2141) except that they were hatched and raised at 20°C and harvested after 55 hours. After fixation and DAPI staining (Pepper et al., 2003a), gonad arms were imaged (Michaelson et al., 2010) and scored as fertile (normal pattern) or sterile (sperm-only Glp-1-like phenotype). daf-2 mutant animals were raised at 15°C and synchronized by L1 hatch-off (Pepper et al., 2003a). For 20°C conditions, hatched animals were immediately transferred to 20°C and scored at mid-late L4, early adult [staged as in Michaelson et al. (Michaelson et al., 2010)] and older adult (24 hours post-mid-L4 at 20°C). For temperature shifts, synchronized animals were raised at 15°C until early L3 (to avoid dauer formation) or early adult stage, and then transferred to 25°C. Older adult after L3 shift was 18.5 hours after the mid-L4 stage at 25°C (Hirsh et al., 1976). After fixation and DAPI staining, gonad arms were imaged and scored for the number of nuclei in the proliferative zone, the distance to the transition zone, and mitotic index as described (Michaelson et al., 2010).

We used MATLAB (MathWorks) to simulate the theoretical scenarios in Fig. 4A and to calculate the average distance of cells from the distal tip (see supplementary material Appendix S1). We averaged the distance of 15 individual precursor cells from the distal end up to 25 CD over 70 time steps. For Fig. 4B, samples from the Statechart-based model were taken each second throughout 4-minute simulations and cell counts of all cells within 25 CD from the distal tip were calculated. Three independent runs gave similar results.

We used Statecharts to specify a cellular decision-making program that captures germ cell behavior as defined by a subset of cell components and their dynamic interactions. These cells propagated within spatial constraints to simulate germline development within the gonad. In this approach, the cell population behavior emerged as a consequence of decisions made by individual cells. Based on a subset of the known developmentally regulated aspects of morphology, GLP-1/Notch signaling and proliferation, we built a baseline design that captured the general developmental progression of germline development (supplementary material Appendix S1, Movie 1; see below). This design specifies the cell as a basic information-processing unit that senses the extracellular space and acts on it over time. This design is an abstraction of the behavior of cells as the key unit within tissues (Cohen and Harel, 2007; Edelman et al., 2010).

This model of the C. elegans germ line allowed us to independently manipulate parameters including ligand interaction, receptor degradation, and influences on cell division and movement. Moreover, changes in these parameters were made at specific time points in the simulation. During the course of a simulation run, different aspects of cell behavior (e.g. ligand-receptor interaction and cell cycle) occur simultaneously as they do in the natural state (see below and supplementary material Appendix S1).

The model’s Statechart designs are organized by elements that specify: (1) ‘effectors’, comprising GLP-1/Notch pathway effector components and their activity; (2) ‘membrane’, comprising ligand-receptor interaction, receptor activity, receptor degradation and cell movement control; and (3) ‘cell’, comprising controls on proliferation and differentiation (see supplementary material Appendix S1). Here, we present these features in terms the model’s representation of space, movement, GLP-1/Notch signaling, proliferation/differentiation and time, as well as a baseline design. We note that the model is not intended to be comprehensive and simulations are not entirely realistic. Rather, they represent a subset of key events to a reasonable level of abstraction.

Although the C. elegans hermaphrodite germ line develops within a three-dimensional (3D) tube-like structure, the behaviors we wished to capture could be represented along the distal-proximal axis (the long axis of the cylinder). Therefore, we modeled within the simpler 2D space, using cell numbers similar to those observed in a surface 2D slice of the 3D structure (Fig. 1). Growth of the tube and cell numbers were based on live worm measurements. As the simulation advances, the gonad grows, filling a predefined gonad space within the overall grid and replicating the characteristic U-shape within which germ cells move and proliferate. With the exception of sperm-oocyte interactions at fertilization, cells can only move into free pixels, and they continually update their location within the grid during a simulation. A cell that vacates a grid pixel area changes the label of the pixel to ‘free’. This feature is important for the design of cell movement.

To simulate germ cell dynamics within the gonad, cells continuously seek possible moves by checking the x and y coordinates and, if applicable, decide to move left, right, up or down in the 2D space. When more than one direction of movement or division orientation is possible, this is a probabilistic decision. In larval stages of the simulation, cells fill the space and remain associated with the distal-most boundary of the gonad. In the adult stage of the simulation, proximal movement is prevalent due to space constraints and a higher probability for proximal movement (see supplementary material Appendix S1).

The ‘membrane’ element includes the GLP-1/Notch receptor component that interacts with an extrinsic ligand (see supplementary material Appendix S1, Fig. S1). When a cell interacts with the ligand LAG-2 in the model, a transition is enabled and the receptor moves from the unbound to the bound state. The LAG-2 interaction is triggered when a cell is positioned within the specified ligand interaction distance (e.g. 3 CD in the baseline design) from the distal end. The receptor active state is connected to the receptor absent state by the transition ‘degradation’, which is triggered when the cell crosses a specified threshold (e.g. 11 CD in the larva and 16 CD in the adult in the baseline design).

The behavior of intracellular pathway components downstream of the receptor (‘effectors’) is described by two states per effector, inactive and active, connected by transitions. All effector components are orthogonal states and thus run in parallel, so a cell is simultaneously in either the inactive state or the active state for each of the effector components. The effectors include: LAG-1, GLD-1 and GLD-2. LAG-1 is initially in an inactive state and its transition to the active state occurs in response to GLP-1 activation. Similarly, GLP-1 degradation inhibits LAG-1 activity. GLD-1 and GLD-2 components are negatively regulated by LAG-1 activity; they are inhibited if LAG-1 becomes active and promoted if LAG-1 becomes inactive.

We also included a mechanism (‘memory’) to allow the consequences of ligand interaction to be sustained and passed from mother to daughter cell (the daughter is designated as the cell that moves into a free position and the mother is the cell that retains its location during the division; both cells are free to move after the division). With memory implemented, in all but rare cases, once a cell interacts with the ligand it retains the memory of the interaction throughout subsequent movements by transitioning the receptor from unbound to bound regardless of the specified ligand interaction distance. Also, the ligand interaction status of the mother is inherited by the daughter. Once a cell crosses the threshold for receptor degradation, it differentiates. The current implementation includes a rare possibility for cells to differentiate prior to updating their receptor activity (for a more detailed explanation of memory and cell division see supplementary material Appendix S1). Memory is active in all simulations presented here. In the absence of memory, each time cells move they reassess whether they are within the assigned spatial limit for ligand interaction. In simulations without memory, the proliferative zone was maintained with ligand interaction distance at 16 CD and receptor degradation at 16 CD, but not for any other combinations of distances tested (3, 8 and 16 CD for each; data not shown). Roughly speaking, memory acts as an amplifier or feedback mechanism, ensuring maintenance of the proliferative zone beyond the position in which cells are in close proximity to the niche.

The initial cell state, ‘precursor’, denotes an undifferentiated germ cell. With the exception of the first few germ cell divisions (see Austin and Kimble, 1987), the ability to proliferate depends on GLP-1/Notch activity. If either GLD-1 or GLD-2 is active, the cell enters the early meiosis state (i.e. leptotene and zygotene). Once a cell enters pachytene, the differentiation component enters the meiotic state. Because differentiation is tied to distance from the distal tip in the model (see below), in situations in which loss of an established proliferative zone occurs, unoccupied space opens in the distal zone in the simulations. In reality, germ cells can differentiate while still occupying distal space.

The proliferation component specifies the mitotic cell cycle and differentiation mutually exclusively. Proliferative cells fill space over time by executing two orthogonal instructions: move and proliferate. An interplay between the cell movement and division maintains the proliferative zone. To facilitate spatial aspects of the modeling, we restricted cell division into an available (free) space, rather than cell division itself creating its own space by pushing its neighbors. This somewhat unrealistic solution nevertheless captures the essential elements of the biology because cells evaluate their ability to move far more frequently (every ∼1-5 mseconds) than they evaluate their ability to divide (every 110 mseconds in larvae and 550 mseconds in adults). Cells move whenever possible, opening space for an adjacent cell to move into the free place or for the daughter of a division to occupy the space. If a cell reaches a proliferation evaluation point but no adjacent space is available, it will try again in the next round of evaluation. Since the frequency with which cells evaluate available neighboring space is far more frequent than that with which they evaluate whether to divide into the available space (division evaluation frequency), the division evaluation frequency serves as a comparative proxy for cell division frequency or cell cycle duration. The 5-fold difference between larva and adult (110 versus 550 mseconds), which was determined empirically by identifying values that produced the proper developmental pattern, reflects a slowing of cell division cycle between larval and adult stages (Maciejowski et al., 2006; Crittenden et al., 2006; Michaelson et al., 2010). Thus, this ratio represents a qualitative change in the cell cycle and might not reflect a predictive ratio.

There is no specific correlation between simulation time and developmental time. As benchmarks, the simulation L1/L2 molt is at ∼5 seconds, L2/L3 at ∼10 seconds, L3/L4 at ∼30 seconds, and the L4/adult at ∼1:00 minutes. On average, this corresponds to 2 seconds/hour of development. In the adult, there are fewer benchmarks because the gonad size is fixed, but fertilized embryos appear at ∼30 seconds after the L4/adult molt, corresponding to 2-3 hours. This averages to 10-15 seconds of simulation/hour of developmental time. We tested simulations up to 10 minutes and observed maintenance of the proliferative zone.

The model is amenable to expansion (for further details, see supplementary material Appendix S1) and is thus a starting point for future work. Here, we designated a baseline design with the following starting points: memory implemented, ligand interaction distance of 3 CD, receptor degradation threshold of 11 CD (larva) and 16 CD (adult), and a division evaluation frequency of 110 and 550 mseconds for larva and adult, respectively. Under these conditions, the simulation produces and maintains a normal developmental pattern (when tested in more than 20 simulations). Starting from this design, we further manipulated pathway component activities, ligand interaction and receptor degradation distances, and division evaluation frequencies.

To test whether the baseline design responds appropriately to genetic perturbations, we altered the design to reflect known mutations in pathway genes (Fig. 1; supplementary material Fig. S1, Appendix S1). In short, we elevated or reduced Notch receptor activity and interfered with both positively acting (e.g. LAG-2 and LAG-1) and negatively acting (e.g. GLD-1 and GLD-2) components (Fig. 1; supplementary material Fig. S1, Appendix S1). In all cases, the design changes produced simulations that agreed with observed mutant phenotypes.

To examine the relationship between signaling and maintenance of the adult progenitor cell population, we ran simulations that altered the distance from the distal tip at which cells experienced ligand interaction and at which the receptor was degraded. We examined all combinations of 3, 8 and 16 CD for ligand interaction and receptor degradation to simulate short-, medium- and long-distance effects of Notch signaling (supplementary material Table S1A).

In one set of simulations we maintained all parameters as in the baseline model and only changed receptor degradation at the adult stage. Simulations with an adult GLP-1/Notch receptor degradation at 16 CD maintained the proliferative zone (Fig. 2A). By contrast, when receptor degradation was set to 8 CD, the zone became unsustainable over time with significant variability from individual to individual. Eventually, all progenitors differentiated. Finally, setting receptor degradation at 3 CD caused rapid loss of the proliferative zone (Fig. 2A; supplementary material Table S1A, Movies 1-3).

In another set of simulations, we altered the ligand interaction distance throughout germline development while maintaining all other parameters as in the baseline model. Under these conditions, a ligand interaction distance of 8 CD or 3 CD developed into a sustainable adult proliferative zone (Fig. 2B). However, when the ligand interaction distance was reduced from 3 to 2 CD (supplementary material Movie 4), considerable instability in progenitor maintenance occurred over time. In these conditions, only 37% of arms maintained the proliferative zone, while in the other 63% all progenitors differentiated [by 3:30 (minutes:seconds) of the simulation]. Interestingly, arms that were stable remained so when taken to time points beyond 3:30. Further, in simulations in which all progenitors differentiated, the time at which the zone ‘crashed’ was highly variable (Fig. 2B). When the ligand interaction distance was further reduced to 1 CD (supplementary material Movie 5), the proliferative zone was never sustainable and all progenitors differentiated on average earlier than in the 2 CD condition. Again, the exact timing of loss of the zone was highly variable (Fig. 2B).

We found that our simulation results were comparable with the effects of a mutation that reduces the activity of the ligand LAG-2. The penetrance of the loss-of-progenitors phenotype at the 2:00 (adult) time point in the simulation resembled that observed in the mutant (15% of arms lost all progenitors in the simulation and 20% in the mutant; Fig. 2C and supplementary material Table S1B). While analyzing the simulations, we also observed individuals in which one gonad arm was stable while the other was unsustainable. Interestingly, this phenotype is also observed in the mutant (Fig. 2D). Our results suggest that specific population dynamics within an individual gonad arm – which is likely to be related to a critical minimum number of progenitors within an individual gonad arm at a specific time – might contribute to an apparent all-or-nothing mutant phenotype. Taken together, our findings suggest that the exact penetrance of mutant phenotypes (the percentage of individuals that show the defect among the total) might be highly dependent on the exact time point at which the phenotype is analyzed.

Cinquin et al. (Cinquin et al., 2010) shifted a temperature-sensitive glp-1 mutant to the restrictive temperature and followed the dynamics of differentiation at subsequent time points. One particularly interesting result was that differentiation occurred progressively in a proximal-to-distal wave from the starting 12 CD to ∼5 CD, after which it occurred more or less synchronously. This result suggested that the distal zone could be divided into two subpopulations of cells based on the response to the loss of glp-1 activity: proximal cells that exhibited a graded response and distal-most cells that respond synchronously. However, the dynamics of the decay of glp-1 activity after the shift are unknown (and are not currently measurable). It therefore remained possible that other factors, such as population dynamics and/or the nature of the loss of temperature-sensitive receptor function, might underlie the observed temporal-spatial pattern of differentiation.

We used our model – in which all the cells have the same internal program – to determine which conditions, if any, would mimic the observed pattern of progressive differentiation followed by nearly synchronous differentiation. We created a starting point and temperature shift similar to that of the Cinquin et al. experiment and observed both intuitive and non-intuitive results (Fig. 3). Intuitively consistent with our model design, a sudden loss of all glp-1 activity depleted the entire zone immediately (0 CD, Fig. 3). Surprisingly, simulations from a design with a reduced area of receptor activity after the shift (receptor degradation threshold of 8 CD) followed the same complex differentiation dynamics as the in vivo temperature-shift experiment (Fig. 3). Specifically, the zone size decreased slowly and then (as indicated by the similar curve shapes and inflection points, most easily seen in the median curve with red circles, Fig. 3) decreased more rapidly within the distal-most 5-6 CD. These results suggest that the complex temporal pattern observed by Cinquin et al. (Cinquin et al., 2010) might result from the specific dynamics of receptor activity decay from the temperature-sensitive allele and/or the effects of a critical population size. If so, these results further imply that a reduction in glp-1 activity in the distal zone could produce complex differentiation dynamics without requiring inherent differences between cells. A direct test of these hypotheses in vivo awaits the feasibility of quantitative cell-specific readouts for Notch activity in this system.

In many stem cell systems, the stem cells adhere to the niche, and this adhesion is important for stem cell maintenance (Voog and Jones, 2010). It is not known whether adhesion is an important aspect of the C. elegans niche-stem cell interaction. Nevertheless, we interrogated our model to determine whether a subset of germ cells might display an emergent bias toward remaining in the vicinity of the niche (that is, a behavior of the cell population as a whole that is ‘emergent’ in that it is not explicitly specified by the model design). The rules governing cell mobility simply allow a cell to move into a vacated space whenever possible (assessing the chance to move every ≤5 mseconds; see model details above).

We first considered several theoretical possibilities for the overall dynamics of cell movement over time as a function of distance from a niche (Fig. 4A). If cell divisions occur such that the mother cell moves away from the niche and the daughter cells retain niche contact, eventually the average distribution of cell position over time is a narrow area half-way between the niche and the differentiation boundary (Fig. 4A, linear). Alternatively, if daughters divide away from the niche and mothers remain close to the niche, the distribution of cells over time averages to two areas: a smaller set of cells that remain close to the niche and a second set that average in a peak half-way from the niche (Fig. 4A, niche). If cell divisions can occur anywhere in the zone, cells will eventually be found in a broad distribution (Fig. 4A, random). Finally, if cell birth is restricted to the niche-adjacent region but then cells move freely, the average distribution is also found half-way from the niche, but over a broader area than in the linear or niche case (Fig. 4A, niche + walk).

We ran multiple simulations of our baseline design and determined the average distribution of cells and found that they approximate a combination of the theoretical niche and random behaviors. Specifically, under conditions when the zone is sustainable (ligand interaction ≥3 CD throughout development), ∼6% of the cells remain niche associated (Fig. 4B), while the others distribute away from the niche under a broad peak (for details of cell tracking in the simulation see supplementary material Appendix S1). As predicted, when the zone is not sustainable (ligand interaction 1 CD), the distal group of cells no longer stays tightly niche associated.

Although the actual pattern of cell divisions will require in vivo verification by cell lineage studies, the model predicts that a small subset of cells (possibly corresponding to a stem cell compartment) normally remains niche associated as a product of the interplay of anatomical constraints and cell division dynamics, perhaps in conjunction with active adhesion.

The model affords an opportunity to probe the relationship between developmentally regulated cell cycle dynamics and zone maintenance over time. Previous work indicates that the average larval cell cycle is faster than that of the adult (Maciejowski et al., 2006; Michaelson et al., 2010), and the difference is at least partially attributable to insulin pathway signaling (Michaelson et al., 2010). Insulin pathway mutants display a reduced larval cell cycle time and, as a result, contain fewer cells in the adult proliferative zone. The adult cell cycle in insulin pathway mutants is approximately the same as in wild type, however, and the zone, although smaller in cell number, is not greatly reduced with respect to the distance between the distal end and the start of differentiation. By contrast, reducing Notch signaling results in a substantial decrease in the distance between the distal tip and the first differentiated cells, but the larval cell cycle is similar to that of wild type (Michaelson et al., 2010) (data not shown). Therefore, although larval cell cycle clearly influences proliferative cell number into adulthood, we did not predict that such cell cycle changes alone would affect the adult differentiation pattern. The model provided a simple means to reassess these conclusions.

The baseline model includes an empirically derived 5:1 ratio of larval:adult division evaluation frequency that gave the proper developmental behavior (Fig. 5A,B). Reducing the frequency to 3.3:1 (supplementary material Movie 6) caused a slight delay in initial meiotic entry in the larva. An equal larval:adult ratio at a higher frequency than in the normal adult (that is, 2.5:2.5 versus 1:1; supplementary material Movie 7) produced a delay in initial meiotic entry and slightly reduced the numbers of proliferative zone cells in the early adult, phenotypes that we had observed previously in insulin pathway mutants (Michaelson et al., 2010). Taken together, these results suggest that not only the larval: adult ratio but also the absolute average cell cycle time is important.

Surprisingly, with 2.5:1 (supplementary material Movie 8) or 1:1 (supplementary material Movie 9) ratios, the zone of the late adult was not sustainable, suggesting that we had possibly missed a secondary effect of the larval insulin pathway on the maintenance of the later adult proliferative zone. These observations prompted us to re-examine later time points in the insulin mutants (Fig. 5). To test the hypothesis that a slower cell cycle time in the larval stages could influence the later adult differentiation pattern, we examined the later adult pattern in two conditions that provided weaker and stronger mutant phenotypes, respectively: the insulin receptor (daf-2) mutant raised continuously at a semi-permissive temperature, and the same mutant shifted to a higher temperature. The model predicts that the weaker mutant will show a reduced proliferative zone (in CD from the distal end) in the later adult and that this defect should be stronger under conditions in which the larval zone is more depleted. We found that the later adult mutants do indeed display a reduced distance between the distal end and the distal-most differentiated cells, and this defect is more pronounced in the more severe condition (Fig. 5; supplementary material Table S2). The zone is not lost completely in the daf-2 mutant, in contrast to what is observed in the model simulation, but this might reflect an unrelated anti-aging effect of loss of daf-2 (Luo et al., 2010). Thus, the model provided a testable hypothesis regarding the dynamics of the interplay between proliferation and differentiation and the data support the hypothesis.

Our model provided a hypothesis-generating dynamic view of the C. elegans germline progenitor population. In particular, the ability to study dynamic variability and track both single-cell and cell population behavior proved useful. The model recapitulated incomplete penetrance of all-or-none mutant phenotypes in lag-2 mutants and provided alternate considerations for complex temporospatial patterns of differentiation. It also offered a framework for considering the possible role of population dynamics in niche residence. Lastly, analysis of model simulations predicted that defects in larval proliferation may lead to reduced adult zone maintenance, an unexpected behavior that we validated with laboratory studies.

Our model complements a variety of modeling efforts aimed at understanding and predicting stem cell population dynamics. Wang et al. (Wang et al., 2010) employed a general mathematical model based on differential equations to estimate the number of epithelial stem cells. Others have used mathematical models to infer patterns of stem cell divisions from clone size distributions in order to analyze systems in which cell labeling and tracking are available (Clayton et al., 2007; Lopez-Garcia et al., 2010). Another line of work has focussed on simulation of intestinal crypt stem cells (Gerike et al., 1998; Meineke et al., 2001; Buske et al., 2011). Our approach differs from the latter in that we use a modular visual language (Statecharts, see Materials and methods). In addition, although our approach could be combined with lineage information, it is applicable to systems in which cell labeling and tracking are not yet feasible. To share our model with the community and allow a more rapid development of stem cell models for other systems, we provide a description and the source code plus links to tools under development in supplementary material Appendix S1 (see also http://research.microsoft.com/celegans).

The value of the model derives from its ability to integrate intracellular and intercellular level information and to obtain a dynamic view of the cell population and phenotypes in a way that allows comparison with experimental results. Key features are the ability to conveniently specify a wide range of concurrent processes (signaling, movement, cell cycle, differentiation, division), to test different parameter values and different hypotheses when there is mechanistic uncertainty, and to systematically compare simulations with experimental results. Another useful aspect is the ability to perform in silico mutations and other perturbations. For example, mutations can be mimicked by modifying only the region of the Statechart specifying the corresponding rule. Double- or triple-mutant combinations can similarly be specified. In addition, changes in dynamic behavior (e.g. cell cycle) can be made alone or in combination with mutations. Finally, and perhaps most importantly for cell populations limited to fixed-tissue analysis, the cumulative effects of cell population dynamics over time – resulting from a combination of influences – can be studied using this approach.

Our analysis provides insights into possible general design principles of stem cell systems, including robustness and the influence of developmentally regulated cell cycle control. We observed that a small reduction in the ligand interaction distance throughout development, even in the presence of feedback, had a strong impact on the critical number of progenitors necessary to maintain the zone. This result (loss of the zone and variability in timing of the loss) emerges from the interplay of cell cycle, receptor activity and movements. It is likely that these factors together (as opposed to molecular threshold effects alone, for example) also contribute to incompletely penetrant all-or-none phenotypes. The model also suggests that developmentally regulated cell cycle changes might be an important design principle that impacts the maintenance of stem cell systems. It will be of interest to determine whether these principles are relevant to other stem cell systems.

In short, our approach provides a powerful tool for furthering our understanding of stem cell population dynamics. In silico ‘experiments’ can be performed rapidly, and multiple influences of anatomy, cell cycle and cell fate can be manipulated independently. Although modeling cannot replace traditional experimental approaches, in cell population systems that are influenced by multiple parameters over time it might prove particularly useful.

We thank Hubbard lab members for discussion.