Mechanistic regulation of planarian shape during growth and degrowth

Multicellular organisms regulate the growth of their whole body across different axes to maintain appropriate shape, size and proportions, a process that is not completely understood (Harmansa and Lecuit, 2021; Kicheva and Briscoe, 2023; Texada et al., 2020). Planarians are ideal model organisms for studying shape regulation owing to their extraordinary plasticity; they can grow or degrow (shrink) their adult bodies as they are fed or starved, respectively (Birkholz et al., 2019; Felix et al., 2019; Lobo et al., 2012). During growth and degrowth, the planarian body dynamically changes its length at a linear rate (Almuedo-Castillo et al., 2014; Oviedo et al., 2003; Schad and Petersen, 2020) and, correspondingly, its body area at an exponential rate (González-Estévez et al., 2012a; Thommen et al., 2019). Although several essential signaling pathways in planarians balance mitosis and apoptosis depending on nutrition intake (González-Estévez, et al., 2012a; Miller and Newmark, 2012; Pascual-Carreras et al., 2020; 2021; Peiris et al., 2012; Ziman et al., 2020), how their body proportions change during growth and degrowth and which mechanism coordinates this differential growth across the axes remains poorly understood.

Whole-body rescaling in planarians is mainly due to changes in cell number rather than in cell size (Baguñà et al., 1990) and is facilitated by a population of pluripotent stem cells called neoblasts that are distributed throughout the body of the worm (Scimone et al., 2014). During growth and degrowth, the differentiation of neoblasts results in a constant ratio of specific cell types (Hill and Petersen, 2015; Oviedo et al., 2003; Takeda et al., 2009). Crucially, neoblasts and their progeny are regulated by position control genes (PCGs) that form graded morphogen signals across the different body axes (Bonar et al., 2022; Reddien, 2021; Witchley et al., 2013). PCGs in the Wnt pathway are necessary to correctly define the anterior-posterior (AP) axis (Gurley et al., 2008; Scimone et al., 2016; Stückemann et al., 2017; Sureda-Gómez et al., 2016), slit and wnt5 are necessary to establish the medio-lateral (ML) axis (Adell et al., 2009; Cebrià et al., 2007; Gurley et al., 2010), and bmp4 to pattern and coordinate the dorsal-ventral (DV) axis (Clark and Petersen, 2023 preprint; Molina et al., 2007; Reddien et al., 2007). However, understanding at a mechanistic level how these signals coordinate differential growth across different axes remains a current challenge.

Mechanistic models – specifically, those based on mathematical descriptions of the underlying causes – have been proposed for explaining planarian body patterning during homeostasis and regeneration. These dynamic models can provide mechanistic hypotheses for the formation of the planarian poles and AP patterning by reaction-diffusion (Meinhardt, 1982; Schiffmann, 2011), including the rescaling of head and tail patterns (Werner et al., 2015) as well as the location of fission planes (Herath and Lobo, 2020). Conversely, inhibitory signals diffusing from the worm AP-ML border can produce the planarian midline gradient forming from the ML axis (Meinhardt, 2004). In addition, models explaining head-trunk-tail patterning by PCGs and other morphogens can explain correct and aberrant body patterns after genetic, pharmacological and surgical manipulations (Lobo et al., 2016; Lobo and Levin, 2015; 2017; Pietak et al., 2019), as described in formalized experimental datasets (Lobo et al., 2013). However, although these mechanistic models can predict the formation and re-establishment of key patterning signals controlling body morphology, even during changes in body size, no mechanistic model exists for explaining how these signals regulate and coordinate the dynamics of planarian whole-body shapes.

The regulation of biological shape is a complex problem owing to the dynamic feedback loop between diffusible regulatory signals and the tissue locations where they act (Sharpe, 2017). In planarians, PCG signals regulate growth and differentiation, which, in turn, determine the expression locations and diffusion domains of the PCG signals themselves. This feedback regulation involves dynamic elements at multiple scales, including cellular mechanical forces, tissue growth, and diffusion of signals, and hence systems-level approaches are required to gain a mechanistic understanding (Dillon et al., 2003; Germann et al., 2019; Marin-Riera et al., 2015; Mirams et al., 2013), a further challenge when considering whole organisms. Additionally, for such mechanistic models to recapitulate experimental phenotypes precisely, numerical values for all unknown parameters are needed. Thus, computational inference methods must be combined with the simulation of multiscale models in order to calibrate such models with relevant experimental data (Crocker et al., 2016; Francois and Siggia, 2010; Ko et al., 2022; Mousavi et al., 2021; Schnell et al., 2007; Uzkudun et al., 2015).

Here, we experimentally study the dynamics of whole-body proportions in growing and degrowing planarians and demonstrate how the observed shapes can result from the feedback loop between pole-border signaling and whole-body shape. We show how planarians scale isometrically at similar rates during both growth and degrowth yet maintain different proportions. We propose and study a mechanistic regulatory model of planarian AP-ML body proportions based on the diffusion of pole and border signals regulating tissue growth and the resultant body shape. Using a novel parameter optimization approach combining descriptive and mechanistic multi-level simulations, we successfully calibrated the proposed model to recapitulate whole-body shape dynamics during planarian growth. Then, we employed the model to study how this regulatory mechanism can also explain planarian degrowth dynamics. We show how an increase in apoptosis and pole signaling can explain the transition from the whole-body shape dynamics seen during growth to those observed during starvation-induced degrowth. More generally, this study demonstrates how the complex feedback between diffusible signals controlling growth and the subsequent emergent tissue spatial dynamics can regulate the coordination and scaling of multicellular shape.

Planarians increase their body size when fed (growth) and shrink when starved (degrowth). To understand the dynamics of their shape and body proportions, we imaged for 9 weeks two isolated planarian populations that were either fed (n=20) or starved (n≥36), respectively. Fig. 1 shows a set of representative images of growing (Fig. 1A) and degrowing (Fig. 1B) worms. Next, we analyzed the images with an in-house automatic algorithm to compute the width and length of each worm. The results show that the length and width dynamics each followed a similar linear trend during feeding or starvation, with the magnitude of the slope for both populations being 0.2 mm/week for length and 0.04 mm/week for width, the slopes being positive during growth and negative during degrowth (Fig. 1C,D). Analysis of the body proportion dynamics (Fig. 1E) revealed that the planarian body scaled isometrically, with allometric coefficients of 0.91 during growth and 0.98 during degrowth. Furthermore, the linear slope of the body proportions during feeding and starvation were 0.17 and 0.18 width/length, respectively, which are not significantly different (P=0.624, one-way ANCOVA). However, the body proportions during feeding and starvation were significantly different (P<0.0001, one-way ANCOVA), with fed worms being 0.2 mm wider than starved worms. To assess the transition in proportions from growth to degrowth, a third population of fed worms of a variety of sizes (lengths of 2.6-7.8 mm) were isolated from a regularly fed colony and then starved for 3 weeks (Fig. 1F). The results show that the initial relative width in this fed population (n=47) was not significantly different (P=0.056, unpaired, two-tailed Student's t-test) from that of the worms during the previous 9-week-period feeding experiment (n=200). However, when starved they rapidly transitioned to degrowth proportions. After 3 weeks of starvation (n=33) their relative widths were not significantly different (P=0.018, unpaired, two-tailed Student's t-test) from those of the worms during the previous 9-week-period starvation experiment (n=389). The difference in relative width between the fed and starved worms after 3 weeks was significantly different (P<0.001, unpaired, two-tailed Student's t-test).

The inherent variability of the planarian body shape across different worms, together with the discontinuous timing of the microscopy images (as taken once per week), hinders our ability to extract mechanistic knowledge automatically from the experimental data. To standardize the planarian body shape during growth and degrowth, we developed a continuous descriptive model of worm shape dynamics. The model abstracts the planarian body along the AP and ML axes as an obround geometrical shape (also called a stadium), which is composed of a rectangle with two semicircular end caps (Fig. 2A). This standardized shape is specified with two parameters: the width w and length l. Then, the corresponding obround is defined as the set of all points within some radius, ⁠, of a line segment of length a=l − 2r. Continuous obround functions were defined using the linear regression for width and length obtained from the experimental data of both growth and degrowth (Fig. 1C,D). In this way, the descriptive model based on these fitted continuous equations can recapitulate the dynamic planarian shapes observed during both growth (Fig. 2B) and degrowth (Fig. 2C).

The descriptive model fitted to the experimental data recapitulates in a continuous fashion the planarian shape dynamics during growth and degrowth, but does not explain the mechanism controlling such dynamics. Hence, we sought to develop a mechanistic model that can explain the regulatory interactions controlling the dynamics of planarian shape. We hypothesized that this regulation involves feedback between morphogen signals, controlling tissue growth, and the existing whole-body tissue, providing a domain through which such morphogens can diffuse. To test this hypothesis, we developed a mathematical model based on partial differential equations incorporating both planarian whole-body shape and diffusible morphogen signals (Fig. 2D). We modeled the cells in the planarian body as a continuous tissue with variable cell density held together by cell–cell adhesion and through which morphogens can diffuse (Ko and Lobo, 2019). The model includes morphogens expressed by planarian AP pole organizers, similar to how genes in the planarian Wnt pathway establish expression domains at the poles in intact worms and shortly after amputation (Gurley et al., 2010; Petersen and Reddien, 2009, 2011). The model also includes a planarian border organizer similar to the expression pattern of multiple planarian genes, including wnt5 (Cebrià et al., 2007; Gurley et al., 2010). In addition, all cells can express a morphogen that function as a growth signal, inducing tissue growth. Thus, we hypothesize that four morphogens regulate each other in a feedforward pathway: the anterior and posterior poles (A and P in the diagram) inhibit the border morphogen (B), which in turn inhibits the growth morphogen (G). This simple regulatory pathway can establish the planarian AP and ML positional information gradients and produce a growth morphogen expression pattern that is high along the planarian midline and low elsewhere, inducing the elongated whole-body shape characteristic to planarians. Crucially, this whole-body shape feeds back into the patterns formed by the morphogens themselves as both their diffusion domain and the location of the organizers directly depend on the shape and size of the worm. We propose that these regulatory interactions could explain the spatial dynamics of the planarian shape during growth and degrowth.

We simulated the model to test its capacity to produce AP and ML gradients as well as the elongated growing shapes similar to planarian worms. The initial state for cell density was generated based on the shape of the descriptive model of planarian growth. The parameters for the model were obtained from the literature or manually estimated by trial and error until the mechanistic model produced shape dynamics that qualitatively resembled elongated growth dynamics (Table S1). Fig. 2E,F shows snapshots of the model during a time-course simulation, demonstrating its capacity to produce growing elongated body shapes as well as AP and ML gradients (see also Movie 1). The expression pattern of the growth morphogen is higher at the midline and decreases towards the edges as a consequence of the inhibition by the border morphogen. The anterior and posterior morphogen gradients inhibit the expression of the border morphogen, resulting in higher tissue growth at the poles and hence an elongated shape. However, the resultant simulated whole-body shape deviates substantially from the descriptive model (Fig. 2B). Manually finding the parameters that can precisely recapitulate the observed planarian shape dynamics is extremely difficult owing to the multiple levels of regulation and feedback loops in the model. Instead, to optimize the parameters, we developed a parameter optimization approach that combines the descriptive and mechanistic models.

To automate the discovery of model parameters that can recapitulate the observed planarian whole-body shape dynamics, we developed a novel optimization approach based on comparison of descriptive and mechanistic models, as well as early termination of poor simulations to improve performance. The goal of the method is to find a parameter set for the mechanistic model that, when simulated, recapitulates the shape and size dynamics during 9 weeks of growth as encoded by the descriptive model. Thus, the output of this algorithm is a calibrated model that can explain the experimentally observed planarian whole-body dynamics. The method is based on evolutionary computation, an optimization algorithm that mimics biological evolution (Miikkulainen and Forrest, 2021). A random initial population of model parameter sets iteratively evolves via the application of mutation, crossover and selection operators (Fig. 3A; see Materials and Methods for details).

To quantify the ability of a parameter set to recapitulate planarian growth dynamics, we first define the shape error as the difference between cell density in a simulation of the mechanistic model and the shape defined by the descriptive model at the same time point (Fig. 3B). Each parameter set in the evolving population is scored (fitness) as the simulation time before exceeding a specified shape error threshold (Fig. 3C). In this way, simulations of parameter sets with excessive errors are terminated before reaching the total 9 weeks of simulation, increasing the performance of the method. The shape error threshold decreases dynamically during evolution when parameter sets are found that can complete the whole simulation period and hence stay below the current shape error threshold. Therefore, the average simulated time increases as better model parameter sets evolve, but decreases when a lower shape threshold is selected after a parameter set completes the 9 weeks of simulation (Fig. 3D). This dynamic shape threshold results in more stringent simulations as evolution progresses, balancing evolutionary pressure with an enhanced performance of the method.

To assess the ability of the proposed parameter optimization methodology to find model parameter sets that can recapitulate the planarian growth dynamics, we ran the algorithm ten independent times for 10 days of wall clock time each. Note that, although the mechanistic model is fully deterministic, the evolutionary algorithm includes stochasticity in the crossover and mutation operators. Each of the ten runs converged to a parameter set with maximum shape error below ∼0.1. The results demonstrate how the optimization methodology ensures that the shape error of the best model in the population decreases monotonically during evolution and that all runs converge to model parameters that closely recapitulate the shape dynamics of the descriptive model (Fig. 3E).

We simulated and analyzed the mechanistic model using the parameter set with the lowest error discovered by the parameter optimization methodology (Table S2). The time-course simulation of planarian growth shows how the anterior and posterior morphogen gradients (Fig. 4A) reduce the expression of the border morphogen near the poles (Fig. 4B, red). As a result, the border morphogen inhibits the expression of the growth morphogen all around the border of the worm except near the anterior and posterior poles (Fig. 4B, green). Comparing the simulated mechanistic shape (Fig. 4C) with the descriptive model shape (Fig. 4D) during growth reveals a close match during all 9 weeks of growth (Fig. 4E).

Next, we compared the planarian length, width and body proportions resulting from the mechanistic model during growth directly to the experimental measurements of planarian worms. These results show that the mechanistic model correctly recapitulated the width and length dynamics during the 9 weeks of growth experimental data (Fig. 4F,G). Furthermore, the planarian body proportion dynamics were not significantly different between the mechanistic model and experimental measurements (Fig. 4H), both in terms of rate and proportions (P=0.340 and P=0.287, respectively; one-way ANCOVA). Thus, the mechanistic model with the automatically discovered parameters is able to reproduce the planarian whole-body shape dynamics during growth.

We performed a sensitivity analysis to determine the robustness of the mechanistic model of planarian growth and the effect of each of its parameters on the resulting shape. We perturbed each parameter calibrated with the parameter optimization method by applying factors between 0.5 and 1.5. For each perturbation, the growth simulation was repeated from the original initial state and the resultant whole-body shape analyzed after three simulated weeks (Fig. S1). The resultant whole-body shapes show how each of the individual parameters can affect different morphological aspects of the growing worm. Several parameters, such as the diffusion rate of the anterior and posterior morphogens (D_AP) or the growth morphogen production rate constant (b_G), can control the convexity of the shape. In addition, parameters can also affect the overall length and width of the simulated worm, with apoptosis (λ), dispersion (k_p), and cell adhesion strength (k_a) having the strongest effect (Fig. S2). In this way, the mechanistic model is robust and can produce a variety of growing shapes depending on the values of the different biological parameters.

We sought to validate the calibrated model trained with the planarian experimental growth data by recapitulating the observed dynamics during starvation (Fig. 1). When planarians are starved, the rate of cell apoptosis increases (Pellettieri et al., 2010) and the rate of cell mitosis decreases (González-Estévez et al., 2012a), compared with fed worms. Indeed, an increase in the rate of cell apoptosis, λ_, in the model results in simulations with a decrease in body size, as shown by the sensitivity analysis (Figs S1 and S2). However, this perturbation alone is not sufficient to recapitulate the observed shape dynamics from the experimental data of starved worms, as it distorts the convexity of body shape. Hence, we performed simulations varying the rate of apoptosis, λ_, together with one other parameter in the model, testing all such combinations (Fig. S3). The results showed that only varying both the rate of apoptosis, λ_, together with the pole coefficient, α (which modulates the regulation of the border morphogen by the pole morphogens), produced degrowing obround whole-body shapes during the simulation. Hence, we hypothesized that increasing both the rate of apoptosis, λ_, and the pole coefficient, α_, in the calibrated mechanistic model could be sufficient to cause the observed degrowth shape dynamics during starvation. To test the hypothesis, we performed a parameter scan varying both the rate of apoptosis, λ (increasing up to 10%, at 2% increments), and the effect of the pole signaling, α (increasing from − 1 to 1, at 0.2 increments), for a total of 66 simulations, each starting with the body shape and size defined by the descriptive model for the planarian degrowth experiments (Fig. 2C).

The results of the parameter scan show how the calibrated model could produce a variety of degrowth behaviors with different body proportion dynamics (Fig. 5A). Increasing the pole coefficient, α_, results in higher production of the border signal B as regulated by the pole morphogens A and P (Fig. 5B), and hence a localized lower tissue growth at the poles that changes the resulting whole-body shape. The pole coefficient correlated more strongly with the resulting body length (Fig. 5C), whereas the rate of apoptosis influenced more strongly the body width (Fig. 5D). Furthermore, the simulations that resulted in correct body proportions also matched the body shape and size defined by the descriptive model derived from the experimental degrowth data (Fig. 5E).

Among all the simulations in the parameter scan, an 8% increase in the rate of apoptosis ,λ_, and an increase of the pole coefficient, α_, to 0.6 resulted in the closest match between the simulated (Fig. 5A, black line) and the experimental (Fig. 5A, red dashed line) whole-body proportion dynamics during degrowth. The simulations with this model also revealed that the whole-body shapes closely matched those in the descriptive model from the experimental degrowth dynamics (Fig. 6A-E). Furthermore, the predicted dynamics with the degrowth model also matched the experimental data of planarian degrowth in terms of length (Fig. 6F) and width (Fig. 6G). In addition, the proportions of the simulated dynamic shapes (Fig. 6H) were not significantly different from those obtained experimentally in terms of both their rates (P=0.220, one-way ANCOVA) and magnitudes (P=0.036, one-way ANCOVA). Thus, we have demonstrated that the proposed model can recapitulate both in vivo growth and degrowth dynamics (Fig. S4), although only experimental data of growth was used for the parameter optimization.

Here, we experimentally studied the dynamics of planarian body proportions during growth and degrowth and propose a mechanistic mathematical model of how whole-body shape can be regulated. First, we imaged fed or starved planarians in vivo for several weeks and quantified the microscopy images with an automated analysis pipeline. The results showed that their body proportions are significantly different – growing worms being wider than degrowing ones – yet they scale isometrically at the same rate. To quantify the transition from growth to degrowth whole-body shapes, a third population of fed worms of different sizes were subjected to starvation, resulting in a rapid reconfiguration from growth proportions to degrowth proportions. Next, we created a descriptive model to standardize the changes in planarian whole-body shape and size observed during growth and degrowth. Then, we developed a computational approach to understand the regulatory mechanisms that could control the differential tissue growth and degrowth dynamics required to produce the observed whole-body shape dynamics. Owing to the complex feedback between shape and the expression and diffusion of signaling morphogens, we employed a multi-level mathematical approach whereby cells are modeled as a dynamic continuum. In the proposed computational model, the tissue shape is subjected to cell adhesion forces and can exhibit growth in response to morphogen concentration and degrowth due to cell apoptosis. In turn, the locations where regulatory morphogens are produced and the diffusive domains of those morphogens depend on the dynamic shape of the tissue, completing a mechanochemical feedback loop. We showed how this mechanistic model including the regulation of growth signals by border and AP-pole morphogens could produce growing obround-like shapes.

Next, to test whether the mechanistic model could recapitulate the observed planarian shape dynamics, we developed an evolutionary computation approach to calibrate the model parameters. For this, we designed a novel fitness function combining descriptive and mechanistic models and optimized for computationally demanding simulations. Prior work has used the fast Fourier transform to extract radial coordinates to analyze yeast cell shapes (Chen et al., 2014), but this method is restricted to circular morphologies. Additionally, morphometrics studies of fungi, fish and flies have successfully used Procrustes analysis (Corral and Aguirre, 2019; Klingenberg and McIntyre, 1998; Perner, 2018), but this requires expert labeling of key landmarks and assumes that the organism is scale invariant. In contrast, our data suggest that planarians become proportionally thinner on the ML axis as they grow. Hence, as a generic solution for calibrating models of shape, we defined a descriptive model that abstracts the organism's whole-body shape data as a continuous geometric shape that changes over time. In particular, we used an obround to capture the shape dynamics of the planarian worm. Future work could extend this standardization approach to gene expression patterns, for further calibration of the model (Lobo et al., 2013; Roy et al., 2020). Taking this continuous descriptive model of shape as input, the proposed fitness function minimizes, through simulation, the difference between whole-body shapes produced by the descriptive and mechanistic models of planarian growth. Crucially, the fitness function makes use of a dynamic threshold approach that allows the computation of a model fitness without the need to complete the whole simulation. This computational methodology was successful, resulting in a calibrated mechanistic model that can recapitulate the planarian proportions initially observed during growth.

To understand how planarians shift between regulating growth when fed to degrowth when starved, we used the calibrated mechanistic model to perform simulations increasing cell apoptosis and varying one other model parameter at a time. This parameter scan showed that an increase in cell apoptosis alone cannot explain the transition to degrowth shape dynamics during starvation. Instead, the combination of increasing apoptosis and pole signaling was the only combination of two model parameters that resulted in a simulation that recapitulated the observed initial planarian proportions during starvation, suggesting that the positional control genes patterning the AP axis are actively regulated during degrowth. Furthermore, the model predicts spatially uniform rates for both tissue growth and apoptosis during growth and degrowth (Fig. S5), which are similar to the uniform distribution of mitosis (González-Estévez et al., 2012a,b) and apoptosis (Pellettieri et al., 2010) observed in vivo. The differences in these rates during planarian growth and degrowth result in either a positive or negative net growth, respectively, at the lateral border and AP poles (Fig. S5C,C′). These localized net growth differences are higher at the AP poles than the border, which results in the elongated planarian whole-body shape observed in the growth and degrowth simulations. Thus, even though planarian proportions during growth and degrowth are significantly different, the same proposed computational mechanistic model serves as a combined explanation of how the planarian body shape is regulated during both growth and degrowth.

The proposed pole and border morphogens in the mechanistic model could be implemented by several planarian genes that act as axis organizers (Sureda-Gomez and Adell, 2019). In particular, diffusive Wnt ligands are expressed at the planarian anterior and posterior ends (Gurley et al., 2010; Petersen and Reddien, 2009; 2011), similar to the pole organizers proposed in the mechanistic model. The non-canonical Wnt ligand wnt5 is expressed at the planarian AP-ML border (Cebrià et al., 2007; Gurley et al., 2010), like the border signals included in the model, and is antagonistic to medial slit. Finally, pathways such as EGFR and TOR control planarian cell proliferation (Fraguas et al., 2011; González-Estévez et al., 2012b; Peiris et al., 2012), and, as proposed in the model, may need to be coordinated with AP pole genes in the Wnt (and dorsoventral BMP) pathway controlling differential growth (Clark and Petersen, 2 preprint023).

Although the proposed mechanistic model can recapitulate the changes in planarian shapes and sizes during growth and degrowth, it includes simplifications to facilitate the mathematical simulations. First, the model is based on direct pathways from organizers to tissue growth. However, the equivalent planarian pathways in vivo likely involve multiple interacting genes and signals that also participate in other processes, such as cell differentiation and body patterning. In order to capture this level of regulation, the model could be extended with Turing systems (Turing, 1952), which has been proposed as a self-regulatory mechanism for establishing planarian positional information gradients (Herath and Lobo, 2020; Meinhardt, 1982; Stückemann et al., 2017) and for the patterning scaling with respect to whole-body size (Werner et al., 2015). The model also excludes cell migration for simplicity and instead abstracts tissue growth as the effect of both cell proliferation and migration. Indeed, the differentiation of neoblasts involves their migration towards the target tissues (Lapan and Reddien, 2011), effectively increasing cell density in those locations and hence resulting in tissue growth. Finally, the model is focused on the AP and ML axes, but future work could extend this approach to three dimensions to incorporate the important contribution of DV patterning in the regulation of planarian whole-body shapes (Gaviño and Reddien, 2011).

In conclusion, this study determined experimentally the different dynamics of planarian whole-body shape during growth and degrowth and developed a computational model suggesting that its regulation could be due to a feedback loop between AP-ML morphogens controlling differential growth and the emergent tissue shapes where they are expressed and diffuse. Our approach combining experimental assays, image analysis, multi-level mathematical modeling, and parameter optimization enabled the development of a computational model of the mechanochemical feedback regulating planarian shape. Crucially, by combining a descriptive and mechanistic model, we generated an intuitive continuous definition of model error grounded in experimental data, which served as the training basis for the inference methodology. The resultant mechanistic model is interpretable and hence could be used to shed light on plausible regulatory mechanisms controlling planarian shape during growth and to predict the process by which it can transition into degrowth. Further work will be essential to validate experimentally the predictions posed by the computational model, but this work could pave the way towards a general understanding of relationship between morphogens and shapes in animal development.

The asexual strain CIW4 of the planarian Schmidtea mediterranea was maintained in 1× Montjuic water in an incubator at 20°C in the dark as previously described (Merryman et al., 2018). Two planarian populations were kept separately, each initially composed of either small animals of ∼3 mm length fed beef liver paste twice per week (n=20), or large animals of ∼6 mm length that were starved (n=50). Planarians that fissioned in the starved population were removed (final n=36). No planarians fissioned in the growth population during the reported period. A third population of fed worms of different sizes (lengths between 2.6 and 7.8 mm) were isolated (n=47) and then starved for 3 weeks; worms that fissioned were removed (final n=33).

Each week, three images (technical replicates) were taken for each planarian when gliding in a straight line using a Zeiss SteREO Discovery.V20 zoom stereoscope with dark-field illumination. Planarian body proportions were automatically measured for each image with a custom pipeline in MATLAB (The MathWorks, Inc.) by binary thresholding and then computing the minimum bounding rectangle of the largest connected component, so small debris and visual artifacts were discarded.

Model parameters were derived from the literature, manually estimated, or calibrated with an automated optimization methodology. Parameters not derived from the literature were first manually estimated by trial and error until the mechanistic model produced elongated shape dynamics (Table S1). Then, those parameters were calibrated with an automated computational methodology to recapitulate the experimental data of planarian growth (Table S2).

Simulations were computed in a two-dimensional domain in a uniform square lattice of 130×130 elements corresponding to 6.5×6.5 mm using the explicit upwind finite volume method with flux limiting and a zero-flux boundary condition. The nonlocal integral term for adhesion is discretized into angular and radial components using bilinear interpolation, as described by Ko and Lobo (2019). The system was numerically solved using ROWMAP (Weiner et al., 1997) and implemented in MATLAB (The MathWorks, Inc.). Body proportions of simulations were computed as the bounding box of cell density, where the cell density in border elements signifies sub-spatial-resolution width. The colors in the simulation panels and movies represent each variable normalized separately relative to its maximum value throughout space and time.

We developed a novel parameter optimization approach to infer the model parameter values that can recapitulate the dynamics of experimental whole-body shapes. The method uses a continuous descriptive model as intermediary and is optimized for better performance by terminating low-fitness simulations early. The method is based on evolutionary computation, an optimization algorithm that mimics biological evolution. A population of 36 individuals each representing a different set of model parameters evolve by applying mutation, crossover, and selection operators towards a parameter set that can recapitulate the shape dynamics in the descriptive model, and hence the experimental data.

The initial population consists of 36 random individuals, where each individual consists of a set of parameters, and each parameter value is chosen from separate uniform random distributions (see ranges in Table S2). For each generation (except the first, where the initial population is directly simulated), 18 random pairs of parents are uniformly randomly chosen from the population, such that no individual breeds with itself, but a single individual may breed multiple times or may not breed at all. Each pair of parents produces a pair of children with the same parameter values, except that corresponding parameter values can be swapped between the two children (crossover) or replaced by a new random value (mutation) from the same distribution used to create the initial population (see ranges in Table S2). Both crossover and mutation events can occur with a probability of 10% (Table S3). The parameters for each child are used to simulate the mechanistic model and a fitness score is computed for every individual by comparing the shape error between the simulation and the descriptive model (see below). These 36 new children are then added to the population and out of these 72 individuals, the algorithm selects the 36 with the highest fitness (lowest error) as the population from which the next generation of parents is chosen. The evolutionary algorithm iterates in this way until a predetermined amount of time has elapsed, after which the best parameter set is returned.

To improve the performance of the algorithm, the fitness of an individual is defined as the normalized simulated time before the shape error is higher than some threshold. In this way, simulations with high errors can be terminated as soon as they cross the threshold, resulting in a fitness less than 1, and only those staying within shape errors below the threshold are simulated for all the experimental data, resulting in a fitness of 1. Crucially, the shape error threshold changes dynamically during evolution. It starts with a high value of 1 but decreases when a new individual completes the whole simulation period (fitness of 1). The threshold is then updated to the maximum shape error during the simulation of this individual. In this way, the shape error threshold dynamically decreases as better models evolve and can complete the simulation. This increases the stringency of the fitness computation and hence reduces the computational cost of simulating models while increasing the evolutionary pressure towards better models. Additionally, simulations that have numerical errors or hit the boundary of the domain are also terminated early. Notably, shape error values are only calculated once per individual, because after the threshold changes the new fitness values can be recalculated from the shape error values of the original simulation (as the threshold is strictly decreasing).

We thank Shane Miller and Alejandro Sánchez Alvarado for providing the planarian worms and husbandry protocols. We thank the members of the Lobo Lab and the planarian community for helpful discussions. Computations used the UMBC High Performance Computing Facility (HPCF) supported by the NSF MRI program grants OAC-1726023, CNS-0821258, and CNS-1228778, and the SCREMS program grant DMS-0821311.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.202353.reviewer-comments.pdf