Neural control of body-plan axis in regenerating planaria
Control of axial polarity during regeneration is a crucial open question. We developed a quantitative model of regenerating planaria, which elucidates self-assembly mechanisms of morphogen gradients required for robust body-plan control. The computational model has been developed to predict the fraction of heteromorphoses expected in a population of regenerating planaria fragments subjected to different treatments, and for fragments originating from different regions along the anterior-posterior and medio-lateral axis. This allows for a direct comparison between computational and experimental regeneration outcomes. Vector transport of morphogens was identified as a fundamental requirement to account for virtually scale-free self-assembly of the morphogen gradients observed in planarian homeostasis and regeneration. The model correctly describes altered body-plans following many known experimental manipulations, and accurately predicts outcomes of novel cutting scenarios, which we tested. We show that the vector transport field coincides with the alignment of nerve axons distributed throughout the planarian tissue, and demonstrate that the head-tail axis is controlled by the net polarity of neurons in a regenerating fragment. This model provides a comprehensive framework for mechanistically understanding fundamental aspects of body-plan regulation, and sheds new light on the role of the nervous system in directing growth and form.
Humanity’s inability to regenerate significant portions of anatomy (e.g. lost limbs, severed spinal cord, damaged organs) has prompted decades of intensive study into organisms that can. The planarian flatworm is one such model organism, as it can completely regrow from a piece as small as ~1/250th of the original worm [1, 2]. Crucially, this regenerating organism exhibits complex behaviors and properties functioning well beyond the level of a single cell. It can maintain body-plan homeostasis under normal circumstances, including body-wide allometric remodeling during growth and shrinkage depending on availability of food, while also detecting injury to precisely reproduce missing features, and stopping regeneration once the correct body structure (the target morphology) has been restored [3, 4]. Therefore, a deep and functional understanding of regeneration requires, not only an account of single cell activities such as gene expression [5], but also of the regenerating organism as an intricate system exhibiting complex dynamics and coordination over multiple levels of scale. How do patterning control systems scale to enable correct regeneration of fragments of vastly different sizes and shapes? What tissue-level properties underlie long-range coordination of anatomical outcomes? What is the role of the nervous system in regeneration? All of these open questions have major implications, not only for evolutionary and developmental biology, but also for implementing rational repair strategies in the regenerative medicine of complex organs.
Body-plan consists of precisely-organized collectives of different cell types, where cells receive cues to proliferate and differentiate by responding to morphogens, including molecular-genetic [6–8], mechanical [9–11], and bioelectrical [12–14] signals. Therefore, spatial patterns of morphogens can serve as instructive pre-patterns to induce changes to single cell identity and ultimately specify aspects of the final body-plan [6, 15–18] (Fig 1). Planaria maintain a population of pluripotent adult stem cells (neoblasts), which migrate to wounds to form a blastema, where subsequent proliferation and differentiation ultimately regenerates lost portions of head, trunk, and tail [19–22]. Neoblasts respond to morphogenetic gene products [19, 20, 22–24] and chemical messengers [25–27] to produce the needed cell types to transition the blastema at the former anterior portion of a wound into a new head, while creating tails at wound sites facing the original posterior.
For whole planaria in morphological homeostasis, factors such as β-Catenin (β-Cat), Wingless/Integrated (Wnt), extracellular-regulated receptor kinase (ERK), Notum, as well as other position control genes (PCGs), are found as concentration gradients with characteristic polarities with respect to the anterior-posterior axis of the worm [18, 24, 28–35]. Remarkably, during early regeneration (from 15 to 72 hours after injury), these same morphogens spontaneously reform their original polarity in each fragment of a worm cut into multiple pieces [24, 28–34]. Consistent with the concept of positional information regulating the body-plan, various heteromorphoses can be generated through experimentally-applied genetic and pharmacological interventions, which act by influencing levels of morphogens involved in anterior-posterior axis control [32, 36, 37] (Fig 2, Table 1).
A major knowledge gap exists in determining how wound-induced repolarization of morphogens occurs in each fragment of a cut planaria. Reaction-diffusion mechanisms (e.g. diffusing proteinaceous transcription products or other signaling factors that regulate each other’s expression and activity levels), have long been proposed to account for self-assembly of molecular pre-patterns in a variety of biological systems [6, 15, 17, 43, 44], including the gradients observed in planaria homeostasis and regeneration [45–47]. However, the dependence of reaction-diffusion mechanisms on the size of a biological system—which is an intrinsic feature of most reaction-diffusion models—may limit their ability to account for chemical gradients in planaria fragments of diverse size [48, 49]. Specifically, for a system of interacting molecules with fixed rates of diffusion and other physically-constrained factors, a gradient may spontaneously develop for a system at a certain size, but this pattern will drastically change (e.g. a categorical change in pattern type from smooth gradient to stripes or spots) when the system size increases [48, 49]. This scale-dependence of reaction-diffusion schemes can be mitigated by incorporating directional transport, here termed vector transport, of substances into the reaction-diffusion scheme, which has been shown to robustly account for self-assembly of concentration gradients in a manner essentially independent of size scale [48, 50, 51]. The complex nature of emergent outcomes driven by such processes underscore the essential need to simulate models of regeneration in a spatialized framework to fully determine their large-scale patterning properties, extract testable predictions, and ultimately develop strategies for rational modulation of outcomes in biomedical and synthetic biology settings.
Here we present, and quantitatively analyze, an inclusive model of the molecular regulation underlying planarian regeneration, synthesizing much published work on the results of physiological and molecular-genetic manipulations (Tables 1 and 2). Our model combines a molecular signaling network at the cellular level with vector transport of morphogens (e.g. the directional transport of morphogens by a vector field), and correctly predicts regenerative patterning outcomes in a variety of settings, including editing of regenerative pattern without use of molecular or genetic manipulations. Our model also explains several outstanding puzzles in planarian regeneration, such as how self-assembling, self-scaling morphogen gradients can form in pieces of different size and location across the worm; the mechanistic origin of various heteromorphoses (e.g. one-headed, two-headed, missing-head or missing-tail outcomes) under a range of pharmacological and RNAi interventions; and the instructive role of the nervous system in regeneration. We also tested a number of unique predictions of this model experimentally (Table 2). Hypothesizing that vector transport can take place on the aligned microtubules of nerves enabled the experimental discovery of the crucial role for nerve directionality in determining the AP polarity of fragments, including the re-setting of pre-existing polarity. We show how cooperation of biophysical and molecular-genetic processes can functionally link single cell physiology and long-range anatomical patterning.
Importantly, while this work explores mechanisms enabling body-plan regeneration and target morphology in the context of planarian regeneration (Fig 1), the underlying concepts reveal a powerful, multi-scale patterning system that is applicable to a wide range of biological systems and contexts, such as embryogenesis, regeneration, cancer, and the bioengineering of organoids [52, 53].
Our model is hierarchical and multifactorial in nature, connecting several major concepts. Details of our theory, model, and software are described in S1 Text.
The planaria models presented herein were formally implemented and explored using the Planarian Interface for Modeling Body Organization (PLIMBO), a 1D and 2D finite volume method simulator written in Python3 with open-source tools utilized from Scipy, Numpy, Matplotlib, Scikit-learn [59, 60], and BETSE [61, 62]. PLIMBO was developed specifically to study all aspects of the planaria regeneration modeling reported in this manuscript. PLIMBO allows quantitative testing of the behavior of the regulatory network model reported on herein, in both 1D and 2D contexts, and under a range of experimental perturbations, with vector (convective) transport of morphogens on imported nerve polarity vector fields, as well as the extraction of novel testable predictions. PLIMBO also has the capacity to run parameter searches to automatically iterate parameters to assist model parameterization, to perform sensitivity analyses, and to perform scaling analysis of the model (where a body-shape is progressively scaled); these tools all assist in the development and exploration of complex biological models. PLIMBO is freely available from: https://gitlab.com/betse/plimbo.
Our planarian regeneration model requires the biophysical mechanism of vector transport to create virtually scale-free self-assembly of morphogen gradients instructive for planarian body-plan homeostasis and regeneration, a premise which extends conventional reaction-diffusion mechanisms by adding a convective transport term, the fundamental properties of which have been explored elsewhere [48, 50, 51]. Here we further propose transport of morphogens on polarized microtubules of nerve axons distributed throughout the planarian tissue (i.e. not just the ventral nerve cords, VNC) as the specific vector transport mechanism (Fig 1B, 1C and 1D). Vector transport of Hh on the VNC has been previously hypothesized and qualitatively reported, but not subjected to a formal analysis to determine feasibility and predictions [30, 31, 39].
Experimental investigations into neuronal polarity and axoplasmic transport networks of planaria [63] are reportedly very difficult undertakings, and have produced results that are not useful for the scope required by this present work. Therefore, we estimated maps of nerve axon polarities (Fig 1B), which were derived by manually tracing synapsin stains to determine approximate axon angle at various locations in the planaria tissue cross-section (Fig 1A). The manual traces were assumed to be vector fields supporting directional transport of morphogens (Fig 1B). A summary of all nerve map models is shown in S1 Fig.
Two-dimensional planaria models were created by tracing synapsin stains (S1 Fig columns A and E) to obtain realistic estimates of body shapes, proportions, and axon angles throughout the tissue. Experimentally-derived synapsin stains were straightened along the central midline using an ImageJ feature (’straighten’). In PLIMBO, a clipped and optimized Voronoi grid (1500–2000 grid elements per model) was generated within the domain prescribed by the desired planaria body shape, thereby creating planaria models with a modeled length from 0.5 to 2.5 cm (all models shown in main figures measured 1.20 cm in simulated length, which is consistent with average worm size used in experiments). A finite volume technique was used to model transport of substances within the body (as described in detail elsewhere [61]), where the edge of the planaria body or wound was treated as a closed boundary. Three replicates each of 1H and 2H worm heteromorphoses were generated (S1 Fig). A production gradient, G(x, y), intended to represent the distribution of nerve bodies in the planaria, was used to modulate production of morphogens subjected to axoplasmic transport. The production field was estimated by taking the divergence of the vector transport field, u→(x,y), defining the neural traces (S1 Fig columns C and G). The production gradient G(x, y) was normalized to a maximum of 1.0 and used as a multiplier of chemical dynamics equations. The normalized vector transport field, u→(x,y), in combination with its associated production gradient G(x, y), and the model’s regulatory network (Fig 1E, with details in S1 Fig), produced self-assembling morphogen gradients consistent with observed regeneration outcomes (e.g. Figs 1F, 1G and 2).
A regulatory network at the core of the model included major molecular signals known to be involved in regulation of anterior-posterior axis regeneration in planaria (Fig 1E and Tables 1 and 2). Importantly, the nodes and relationships shown in Fig 1E are consistent with existing experiments described in the literature (see Tables 1 and 2), and represents the simplest model that can account for changes to planarian body-plan, which occur following major known RNAi and pharmacological intervention experiments. The morphogen Hh is assumed to be produced in neurons, and to be moved away from neural cell bodies by kinesin [64, 65], while a presently unidentified substance, Notum Regulating Factor (NRF), is moved back towards neural bodies by dynein [66]. Factors ERK, Wnt, adenomatous polyposis coli (APC), cyclic adenosine monophosphate (cAMP), and Notum are assumed to have the capacity to be expressed homogeneously in all cells. Substances of the regulatory network have production/decay conditions that are often modulated by other substances of the regulatory network, as shown schematically in Fig 1E, with details in S1 Text. The outer edges of the planaria shape (or wound) are considered a closed (i.e. zero-flux) boundary condition. The initial conditions for all molecular factors of the model were set at zero values, except for cAMP which was set at 1.0 for all cells of a simulation not involving a cAMP-based intervention.
Each substance of the regulatory network could contain terms of the following form (shown in Eq 2 for a hypothetical concentration Ci regulated by substances Cj and Ck), where specific equations and free parameters for each substance of the regulatory network are listed in S1 Text: ϕi︸Massflux=-Di∇Ci︸Diffusion+αiu→Ci︸Motorproteinconvection(1) dCidt︸RateofChange=rimax(11+(CjKj)nj)︸GrowthInhibitor((CkKk)nk1+(CkKk)nk)︸GrowthActivator−δiCi(11+(CkKk)nk)︸DecayInhibitor((CjKj)nj1+(CjKj)nj)︸DecayActivator−∇·ϕi︸Divergenceofflux(2)
Concentrations were updated in time using the explicit Euler method. Model stability with time-step (Δt) and minimum grid resolution (Δx) were rigorously tested and results were found to be stable and invariant for time-steps smaller than △t<△x23Dmax, where Dmax represents the largest diffusion constant in the network (Dmax = 1.5x10−11 m2/s). All simulations used time-steps half the size of the stability limit to ensure results were produced in a stable and invariant range.
Detailed information about model parameterization can be found in the S1 Text. Biologically realistic values were chosen for all parameters, based on thorough literature searches (see S1 Text for details). The planaria model was first parameterized in 1D, which features rapid simulation times that allow for high-throughput model testing. Notably, while using a simple 1D transport field, the 1D model shows similar results to that of the 2D model (see Supporting Figures and S1–S4 Videos). After achieving adequate results in 1D, the model was transitioned to the full 2D version using a body shape, production gradient G(x, y) and neural transport map u→(x,y) for one of the models shown in the S1 Fig. A local sensitivity analysis was performed on both 1D and 2D models to identify parameters with the highest impact on model outcome, and to screen out parameters with the lowest impact (see S1 Text for details). The sensitivity analysis also showed the model results to depend strongly on the majority of parameters of the model (see S1 Text for details). Note that units stated on modeled concentrations are model-specific/arbitrary and are intended to be relative to the biological reality, rather than an absolute measure. While concentration magnitudes are relative, they were largely sourced from reports of canonical Wnt/β-Cat signaling pathways [67, 68] and other studies [69, 70].
Cutting of planaria was dynamically simulated by directly removing material from models at steady-state and re-establishing closed boundary conditions at the wounded boundaries. Models were further perturbed by knocking out or altering levels of regulatory network nodes to simulate a variety of RNAi or pharmacological interventions (see Fig 2). To examine the effect of an intervention (e.g. RNAi experiment or pharmacological agent), models were first initialized for 4.5 days of simulated time with original parameter values to obtain steady-state morphogen patterns, at which time the model was cut into fragments, with the introduction of the intervention to the model. All intervention models were examined 4.5 simulated days after cutting, when morphogen patterns were re-established. RNAi simulations (RNAi to Notum, Wnt, Hh, Ptc, and APC) were modeled by setting the growth rate of the inhibited substances to zero (specifically by multiplying their growth constant by 0.0).
Note that due to the long decay time of β-Cat in the system, for a ‘full’ β-Cat RNAi, as shown in Fig 2Ai, inhibition of β-Cat signaling was applied to the whole worm after the initial 7 day model initialization, and a ‘re-initialization’ phase consisting of a further 7 days of simulation of the whole worm with inhibited β-Cat growth was conducted. The re-initialization reduces levels of β-Cat further by allowing them time to decay in the whole worm, and is consistent with experimental technique for RNAi to β-Cat. The model was then cut into fragments and simulated for 4.5 days after the re-initialization phase, with β-Cat expression continued to be inhibited. For ‘partial’ RNAi block of β-Cat, the intervention was introduced at the time of cutting, and was applied with a shorter re-initialization phase of only 48 hours. Pharmacological intervention simulations involved modeling changes to cAMP levels (representing application of pharmacological agents such as IBMX [71], MDL12330A [72], neurotransmitter signaling enhancers such as bromocriptine [42, 73], or direct application of serotonin and dopamine [25, 42]), and block of ERK signaling via U0126, were modeled by introducing the intervention to a 7-day initialized and cut model at the time of cutting, and applying the intervention for 4.5 days of simulated time after cutting. A decrease in cAMP (and application of dopamine, a factor known to reduce cAMP via D2R receptor activation [74]) was modeled as a 0.25x decrease in the growth rate of cAMP, whereas an increase in cAMP (and application of serotonin, a factor known to increase cAMP via 5HTR7 activation [25, 42]) was modeled as a 5x increase in the growth rate of cAMP. ERK was inhibited by multiplying the growth rate of ERK by 0.0. Morphogen patterns were reported 4.5 days after cutting with natural decay rates of all substances preserved.
The resulting model produces stable, self-assembling morphogen gradients as outputs, which examine the functional relationships between morphogen patterns, cell/tissue patterns, and body-plan in homeostasis (e.g. Fig 1F) and with perturbation of the model by cutting events (e.g. Fig 1G). Supporting S3 and S4 Videos show the dynamics of key morphogen gradients during model initialization and after a cutting event in 2D, highlighting the ability for the system to maintain a steady-state, and yet to adapt to a perturbation such as cutting to acquire a new steady-state capable of directing regeneration. Note that the PLIMBO software includes configuration files and information to run all planaria models and associated analyses.
The final component of our model works from the fundamental hypothesis that molecular morphogen gradients ultimately specify the anterior-posterior head/tail decision by changing blastema identity (via morphogen influence on neoblast differentiation). Conceptually, neoblasts located at the blastema of a wounded region were assumed to probabilistically differentiate towards anterior or posterior lineages, thereby shaping the identity of the blastema towards the regeneration of a head or tail at the wound (Fig 1G). Consistent with experimental reports, ERK concentrations were assumed to promote development of an anterior/head identity [24, 30, 31], while β-Cat was assumed to promote establishment of posterior/tail identity [36, 37]. A Markov Chain Model of Regeneration was developed to predict regeneration outcomes in a population of planaria under the influence of steady-state morphogen concentration levels at wound sites of the model. The premise of the Markov Model is summarized in Fig 3. The transition probabilities from blastema to head (αBH) and from blastema to tail (αBT) were estimated in terms of averaged concentrations of ERK and β-Cat (respectively) at wound edges (Fig 3A and 3B): αBH=11+exp(-(CERK-C1κ1))(3) αBT=11+exp(-(CβCat-C2κ2))(4) Where C1 = 0.75 nM, C2 = 300 nM, κ1 = 0.05 nM, and κ2 = 40.0 nM. Values for the reverse transitions were set as constants, where βHB = βTB = 5.0x10−3. These parameter values were chosen based on inspection of ERK and β-Cat morphogen gradients produced by model output under the normal regeneration condition (i.e. without an applied intervention) and choosing values that appropriately prescribe head and tail regions at anterior and posterior wound sites.
At the site of a wound, the probability of a regenerated head (pH), tail (pT), or failure to regenerate either head or tail (pB), were computed in terms of transition probabilities between states under the constraint that all probabilities must sum to 1.0, where the probabilities of a blastema proceeding to different outcomes were dynamically expressed as: dpHdt=(αBH-TαBH-H(αBH+βBH))Cwound(5) dpTdt=(αBT-HαBT-T(αBT+βBT))Cwound(6) pB=1-pH-pT(7) Cwound=1-11+exp(-g(t-tf))(8)
Equations describing anatomical outcome probabilities were solved dynamically (i.e. at each time-step of the iterated model) to yield regeneration probabilities as a function of time, which were sampled for cells at a particular cut line to determine the probability of head/tail at the site of a wound (Fig 3C).
Note that regeneration was assumed to be a time-dependent process induced by wounding and proceeding at the site of the wound. Fixing a time-dependence to regeneration was necessary to describe RNAi interventions such as partial RNAi to β-Cat, where morphogen levels continue to decrease to zero over time. Regeneration time-dependence was modeled by modulating probability for head and tail developments by a molecular wound signal factor (Cwound), which is proposed to be a molecular substance transiently produced by cells at the time of wounding, and to modulate the ability for head or tail probabilities to change.Cwound was described by a simple pulse function (Eq 8), with t representing time after wounding and g = 12 hrs and tf = 90 hrs. The biological identity of the Cwound substance may be a factor such as reactive oxygen species (ROS), which is strongly produced by wounding and has been previously found to be required for planaria regeneration and remodeling to proceed [75, 76]. These parameter values were chosen based on experimental observations that pharmacological interventions applied 72 hrs after amputation had no effect on regeneration.
The probability of morphological regenerates of a particular fragment were then computed in terms of the probabilities of head, tail or regenerative failure at each cut line of the fragment (Fig 3D). An existing, non-wounded head or tail on a fragment was assumed to be terminally differentiated tissue that would not change during the regeneration process, and these head or tail endpoints of the original body were assigned probabilities of 100% head and 100% tail outcomes, respectively (Fig 3C). The above-described Markov Model was used to estimate frequencies of heteromorphoses expected in a particular fragment appearing in a population of regenerates by multiplying the probabilities of each possible outcome at each cut line of a fragment (Fig 3D). The respective predicted fraction of heteromorphoses appearing in the population were inferred from the probabilities of combinations at both cut-lines of a fragment, where the possible heteromorphoses and their probabilities are shown in (Fig 3E). The Markov Model, including identification of fragments after cutting, identification of cells at the wound, calculation of probabilities at the wound, and calculation of heteromorphoses probabilities in each fragments, was performed in an automated fashion in PLIMBO.
Planaria of the species Dugesia japonica were used in all experiments. Planaria were maintained as described in [77]. Animals were kept at 13 °C in Poland Spring water with weekly feedings of calf liver paste and twice-weekly cleanings. Animals were starved for at least 7 days before experiments were performed.
Planaria were cut as described in Nogi et al [78], using a scalpel on moist and cooled layers of filter paper. 2H worms were generated via 3 day incubation of cut fragments in 127 μM octanol as described in [79, 80]. A 2H worm was characterized by at least one eye on each end of the worm [79].
A variety of cutting/amputation scenarios were tested, each of which were designed to cut the proposed neural transport map in key geometries to influence outcomes according to the main hypothesis.
In order to explore the role of proposed VNC cord polarity in directing regeneration outcomes, a nerve-deviation experiment inspired by the work of Kiortsis in Mediterranean Fanworm (Spirographis spallanzanii) [81], was designed. Trunk fragments of planaria were excised and immediately further cut using upwards or inverted L-shaped incisions. Specifically, fragments were cut perpendicular to the head-tail axis halfway from margin to midline of the worm and further cut along the midline towards either the head (upwards) or tail (inverted), causing tissue movement exposing the VNC in a forward- or reverse- polarity, with N = 43 and N = 46 replicates in upwards and inverted L-cut groups, respectively. The cuts were re-enforced every day for 7 days while worms regenerated at 20 °C, to prevent simple wound healing. All worms were monitored once a day for 14 days, with specimens reserved from each group for synapsin staining every other day for 11 days. Computational modeling using the same upwards and inverted L-shaped incisions on planaria models was performed to obtain predictions for these cutting scenarios. Control L-cuts were made similar to the L-cuts above, but not cutting so far from the margin to the midpoint as to including the VNC, thereby creating an exposed side area free of VNC but attached at the base to the main worm, with N = 30.
To explore the role of net nerve alignment in fragments, rectangular fragments with the ~1 mm long-edge corresponding to the original anterior-posterior worm axis (thereby allowing for identification of axis-polarity in the regenerates) were cut to include the ventral nerve cord (VNC-containing) or from the side margin of the worm, excluding the VNC (VNC-free), with N = 94 and N = 103 replicates in the VNC-free and VNC-containing groups, respectively. The fragments were cut from the margin of 1H worms, with the worm placed ventral side up, which allowed the VNC to be seen and to thereby be avoided in the VNC-free fragments. The fragments were checked 1 day after cutting to ensure that all 4 sides were cut. Fragments that did not fulfill this condition were discarded. All fragments were monitored once a day for 14 days, with additional specimens from each group reserved for fixation and synapsin staining at days 1 through 7, 10 and 14. The net orientation of the synapsin signal was determined using the Directionality plugin in ImageJ. Specifically confocal stacks of the whole fragments with a z-resolution of 2 μm were z-stacked using the average intensity projection in ImageJ and oriented so that the long edge of the fragment was vertical in the images. Fragments in which no long edge could be determined were not analysed. The Directionality plugin was run using the Fourier component method, 90 bins and a histogram start at 0. Resulting data was averaged for all fragments, binned into 10° sections and plotted using Excel. Brightfield images and respective synapsin stains taken over multiple days and showing further replicates, can be found in S1 and S2 Datasets.
To explore regeneration outcomes for different cutting scenarios in two-headed worms, 2H worms were generated as described above, and had completed regeneration at least 1 month prior to subsequent amputation. The 2H worms were cut using either decapitation, or cuts through the central symmetry line of the worm, cuts along the diagonal of the worm [80], or a combination of transverse cuts producing fragments. The number of replicates varied in each cutting scenario, with N = 401 used for short fragments containing a head, N = 90 for medium fragments with head and second head amputation cut not crossing the central symmetry line, N = 140 for fragments with one head and second amputation crossing the central symmetry line, N = 90 for dual head amputation containing the central symmetry line, and N = 41 for fragments with dual head amputation and not containing the central symmetry line. Diagonal cuts of 2H worms were performed through the midpoint of the worm (as estimated by the pharyngial opening) and the base of the two heads (N = 97). Cut fragments were allowed to regenerate for 1 week at 20°C and 1 week at 13°C (to prevent fissioning) before scoring regenerative outcomes.