The Role of Bioelectrical Patterns in Regulative Morphogenesis: An Evolutionary Simulation and Validation in Planarian Regeneration

1. Introduction

Morphogenesis, a central aspect in developmental biology, is the processes that governs the formation and development of an organism's large-scale structural features. These processes involve large numbers of cells that create and repair a species-specific target morphology [1,2]. Understanding this process is critical for advancing the study of evolutionary developmental biology and for pioneering new approaches in regenerative medicine. Specifically, knowledge of how groups of cells interpret biophysical signals to determine the anatomical structure they will cooperate toward is crucial for addressing issues like birth defects, recovery from traumatic injuries, degenerative diseases, cancer, and aging.

One complex and poorly understood aspect is the cessation of growth once the desired structure is fully formed. For instance, axolotls can regenerate perfect replicas of lost limbs, tails, eyes, jaws, and other body parts, with growth and remodeling halting precisely when the correct form is achieved [3]. Moreover, regulative morphogenesis can reach the target morphology whether the starting point is a typical anatomical scenario or a novel configuration. For instance, tadpoles that start with experimentally-induced disordered placements of craniofacial features still develop into normal frogs {Vandenberg, 2012, 22411736}{Pinet, 2019, 30974103;Pinet, 2019, 31253636}. This occurs because the eyes, nostrils, and mouth shift along unconventional paths from their initial, misplaced positions to align correctly with standard frog anatomy. This ability is often described as a form of anatomical homeostasis, where deviations from the intended morphology are gradually corrected to achieve the desired anatomical structure, or in other words progressively reducing the error between the current state and a specific anatomical setpoint (the target morphology) [4,5,6].

How do cellular collectives reliably navigate morphospace [7] to form correct anatomical structures despite perturbations? Cells coordinate their actions during morphogenesis through a mix of chemical, biomechanical, and bioelectrical signals [1]. Rather than a precise anatomical blueprint, the genome provides the cellular building blocks, such as proteins, and it's the role of developmental physiology to translate this information into morphogenetic actions. Here, we focus on one layer of that control information: bioelectrical prepatterns [1,8,9,10,11], e.g. spatiotemporal distributions of resting potential, established by ion channels and shaped by electrical synapses known as gap junctions.

Bioelectricity refers to the endogenous electric potentials that are generated within cells and tissues, which are essential for coordinating cellular activities and control the large-scale patterning of anatomical structures [12]. These bioelectric signals, characterized by resting transmembrane voltage potential gradients across all cells, act as critical regulators of cellular behaviors formation [13]. Bioelectric circuits, which are found in both neural and non-neural tissues [2,14] and are facilitated by electrochemical synapses such as gap junctions for non-neural tissues, play a crucial role in regulating large-scale cell growth and overall tissue morphology [2,14,15]. These circuits are essential for making collective cellular decisions about migration, differentiation, and gene expression that affect the entire organism and guide it towards a specific geometry. Additionally, bioelectric signals regulate species-specific anatomical features during critical processes such as embryogenesis and regeneration interacting with molecular-genetic pathways to ensure the proper coordination of tissue morphology and maintain tissue homeostasis [8].

Recent advancements in the study of developmental bioelectricity have enabled specific manipulation of cellular voltage states in vivo by manipulating the endogenous sources and conduits of bioelectrical state, affording much greater specificity than classical techniques dependent on electrodes [13,16,17,18,19]. Optogenetic and chemogenetic tools have been developed to control ion channel activity, and thus V_mem patterns, with high precision [20]. These manipulations can induce regeneration [21,22,23,24,25,26], de novo growth of entire organs such as eyes [27], repair of birth defects induced in the brain, heart, and face by chemical or genetic teratogens [28,29,30,31], and the normalization of cells expressing powerful human oncogenes [12,20,32,33,34].

To increase impact of bioelectric modulation in fields like cancer, regeneration, embryonic malformations and aging [8,35,36], as well as in bioengineering, it is essential to crack the bioelectric code [26,37]. This means, moving beyond existing information on how voltage levels manipulate second messengers and regulate gene expression [10,38,39,40,41,42,43,44,45,46,47] and stem cell behavior [48,49,50,51,52,53] to achieve a computational understanding of how specific bioelectric patterns map onto anatomical outcomes. Such knowledge could then be used to derive efficacious interventions for biomedical settings [18,54,55], as well as a better understanding of the contribution of bioelectric pattern control mechanisms in an evolutionary-developmental biology context [56,57].

The key to understanding any biological code lies in deciphering how it is read or decoded [58,59,60,61] — what information does it contain and how is this information interpreted by cells and tissues? The bioelectric code has been shown to operate through at least three main types of patterns: directly encoding, indirectly encoding, and discrete triggers {Lobo, 2014, 24402915}. Directly encoding pre-patterns map 1:1 with the specific target morphology that they encode. For example, the “Electric Face” is a pre-pattern that arises during embryogenesis of Xenopus laevis and instructively demarcates the location of future tadpole facial structures [62]. It is a direct pattern because the regions of bioelectric state indicate, in a clear way, the actual locations of the craniofacial organs.

In contrast, indirectly encoding pre-patterns do not map to the specific structures that they encode in any obvious manner. For example, a specific V_mem pattern induces the growth and patterning of a new tail [63]. Similarly, a V_mem gradient controls the size [64] and specific shape [65] of planarian head regeneration. In both of these cases, bioelectric pre-patterns specified growth of structures whose geometry could not have been guessed by inspection of the bioelectric pre-pattern – the voltage patterns don’t look like the anatomy they encode.

Finally, trigger pre-patterns have even less similarity to the final outcome: they do not seem to encode specific structures and are very low information-content stimuli. Instead, they can initiate comprehensive morphogenetic programs that result in the formation of complex structures (such as organs and appendages) without requiring micromanagement of the specific structure. For example, inducing sodium influx or proton pumping to tadpole limb or tail amputation sites has been shown to trigger the regeneration of the entire appendage—with correct vasculature, innervation, and muscle—even during periods when such regeneration is not typically possible [21,26]. And, the same stimulus (triggered by the ionophore monensin) induces tails in a tadpole and legs in a froglet – the specificity is not in the signal but in the interpretation machinery triggered by it [26]. These examples indicate that bioelectric patterns can serve as master regulators, kick-starting the growth of complex structures, and all the myriad required molecular processes, from a relatively simple upstream event that itself doesn’t have the informational bandwidth to specify the details.

While computational models exist of bioelectrical dynamics [66,67,68,69,70,71,72], there are currently no computational platforms that facilitate the study of the different modes of bioelectrical encoding in the context of morphogenesis, to generate testable predictions about interventions and these decodings’ evolutionary properties. To bridge these knowledge gaps, we built simulation software that enables exploration of how artificial organisms develop morphogenetic competencies over developmental and evolutionary timescales with the different types of bioelectrical patterns.

We employed evolutionary algorithms to simulate the evolution of a neural cellular automata (NCA) [73] model for regulative morphogenesis, and studied its various properties under the 3 different models of bioelectric encoding. We found that: (1) All three types of bioelectrical codes allow the reaching of target morphologies; (2) The resetting of the bioelectrical pattern and the change in duration of the binary trigger alter morphogenesis; (3) The direct pattern organisms show an emergent robustness to changes in initial anatomical configurations; (4) A The indirect pattern organisms show an emergent robustness to bioelectrical perturbation; (5) Direct and indirect patterns organisms show an emergent generalizability competency to new bioelectrical patterns obtained via rotation of the original ones; (6) Direct patterns organisms show an emergent repatterning competency in post-developmental phase. Because our simulation was fundamentally a homeostatic system seeking to achieve specific goals in anatomical state space, we sought to test the effects of anxiolytics that reduced concern over the distance to those goals. We developed a selective serotonin reuptake inhibitor simulation, which predicted a diminished ability of artificial embryos to correctly interpret bioelectric patterns due to an altered organismal reward machinery, leading to higher variance of developmental outcomes, global morphological degradation, and anatomical bistability. These computational findings were then validated by data collected from planaria regenerating with SSRI exposure.

2. Related Work on Artificial Embryogeny

The field of artificial embryogeny has been historically centered on the French flag model [74,75], with a focus on how cells interpret positional information and differentiate into specific cell types based on morphogen concentrations to form a specific anatomy. This model relies on the existence of thresholds in the morphogen gradient to guide cell fate determination and spatial patterning. The foundational work of Lewis Wolpert introduced the French Flag model, which explained how morphogen gradients can provide positional information to cells, resulting in distinct cellular fates based on concentration thresholds (Wolpert, 1969). Subsequent studies have developed mathematical and computational models to simulate the French Flag pattern formation. For example, Bowers [76] demonstrated the formation of modular structures in a computational evolutionary model of embryogeny, highlighting the phenotypic robustness in generating French flag-like patterns [76]. Recent studies have integrated the French Flag model with other theoretical frameworks in developmental biology, such as the Turing mechanism of pattern formation and gene regulatory networks. These integrative approaches have provided deeper insights into the robustness and flexibility of developmental patterning mechanisms [77]. Other work also included growth and repair simulation models, reaction-diffusion experimental models, and computational frameworks that aim to create robust and tunable axial patterns [77,78,79]. Chavoya and Duthen utilized a genetic algorithm to evolve cellular automata capable of generating various two-dimensional and three-dimensional shapes. They also evolved an artificial regulatory network (ARN) to generate cell patterns, effectively addressing the French flag problem [80,81]. The French flag problem has also been studied in the context of self-organizing paradigms, where the challenge lies in achieving robust and tunable axial patterns without the need for global signaling cues [82]. The French Flag model remains a pivotal framework in developmental biology and artificial embryogeny, facilitating our understanding of how cells interpret positional information to create intricate patterns during development.

Recently, Neural Cellular Automata (NCA) models have been introduced in the field of in silico embryogeny [83,84,85]. Pio-Lopez et al. conducted a study on the scaling of goals from cellular to anatomical homeostasis using evolutionary simulations to resolve the French-flag problem [85]. The results suggested that the collective problem-solving of cells during morphogenesis can evolve into traditional behavioral intelligence by scaling up during evolution homeostatic competencies via stress sharing in metabolic space. Grasso (2022) explored the concept of Empowered Neural Cellular Automata in the context of morphogenesis. The study found that incorporating empowerment as a secondary objective in the evolution of NCA for morphogenesis tasks leads to higher fitness compared to evolving for morphogenesis alone [83]. This highlights the importance of incorporating empowerment as a guiding principle to enhance the efficiency and effectiveness of NCA in performing morphogenetic tasks. Very recently, Hartl et al. developed an NCA implementing a multi-scale competency architecture to resolve a more complex version of the French flag model, the Czech flag and a smiley face as another target morphology [84]. Other works extended modelling approach of morphogenesis with NCAs to the encoding of a manifold of NCA, each of them capable of generating a distinct image [86], growing 3D shapes [87], or hierarchical NCAs for morphogenesis [88].

However, none of these models incorporated the use of the bioelectrical pattern for regulating morphological behavior. Manicka et al. studied how non-neural cells collectively make decisions in the context of morphogenesis, by constructing a minimal biophysical model of a bioelectric network that integrates fundamental components and processes of bioelectrical signaling [89,90]. Although computational models of bioelectrical dynamics exist [66,67,68,69,70,71,72], there is no platform to analyze modes of bioelectrical decoding by cell groups in morphogenesis, or to predict their evolutionary effects. To fill this gap, we created simulation software by extending previous work using evolutionary NCA for morphogenesis [84,85], that explores how artificial organisms develop morphogenetic abilities over developmental and evolutionary times, driven by different types of bioelectrical patterns: direct, indirect and the binary trigger.

3. Materials and Methods

3.1. In-Silico Evolutionary System for Bioelectrically-Regulating Morphological Behavior

We simulated the evolution and development (Figure 1) of 9x9 cells artificial organisms on a morphogenetic task driven by different types of bioelectrical patterns using a genetic algorithm. Using this simulation, we aimed to gain insight on how biological organisms utilize bioelectric patterns during pattern formation and morphogenesis during regulative development, as well as to explore the developmental properties that emerge from organisms evolved with three classes of bioelectric patterns. Due to the large space of possible anatomies that exist in nature, we simplified the model to explore the development of two organisms (the “Tadpole”, inspired by Vandenberg et al.’s tadpole craniofacial rearrangement experiment [91], and the “French Flag”, as described by Wolpert’s original French flag problem [74,75]), each with their own starting and target morphologies. The bioelectric patterns were designed by us, so that we could test our hypotheses on the emergent properties of organisms utilizing various patterns. As with all computational models of biology, we made choices about which features to simulate, and did not explicitly include all the mechanisms by means of which cells communicate, perceive and act on their environment, differentiate, and perform basic metabolic functions. While these details are important to biological systems, we did not require more advanced modeling to answer our specific questions.

The following sections detail the structure and behavior of our artificial organisms and cells, the developmental cycles that enable embryos to reach adult morphologies, the evolutionary cycles by which a population’s best morphogenetic strategies are discovered, and the implementation of bioelectricity as a goal-encoded layer.

3.1.1. Organisms

Each digital organism possesses a distinct artificial neural network [2], similarly to the gene regulatory networks controlling the decision-making pathways in real cells [92]. Organisms are represented as 2-dimensional arrays of cells, depicted as colored circles, with each cell containing a copy of the ANN (see Figure 1D). Organisms have two stages: embryo and adult. The adult phase starts once the target morphology has been reached. Embryos are initialized as a 9 by 9 matrix of cells and undergo a developmental cycle to become an adult. Adults undergo fitness evaluations for selection at 100 steps.

There are two species of artificial organisms in this work: “Tadpole” and “French Flag” (Figure 2A). Tadpoles’ embryonic craniofacial morphology begins as deformed facial features but its adult morphology is a rearranged, symmetrical face, modeling the remarkable observed craniofacial repair of “Picasso tadpoles” [91] and the “electric face” prepattern known to drive face development [62]. The embryonic morphology of French Flag organisms begins as an undifferentiated monolayer of blue cells, which develops into the likeness of the French Flag. This models Wolpert’s famous French Flag problem [75], in which 3 colors are used to represent three distinct cell types in a positional information axis.

3.1.2. Cells

Representation and Attributes

Cells are represented as objects containing energy (a floating point number ranging between 0 and infinity, initialized at 50), molecules (a floating point number greater than or equal to 0, initialized to an amount depending on the starting cell type (Figure 2A)), gap junctions (percentage between 0 (fully closed) to 100 (fully open), initialized at 0), a cell type that depends on the species (Figure 2A), a voltage (integer between 0 and -100, initialization depends), and an ANN [93] (Figure 1A).

An energy threshold greater than zero is required for cells to stay alive. At 0, the cells die. Molecules are exchangeable with neighboring cells and determine cell type. Gap junctions are channels between cells that regulate the flow of molecules according to their opening percentage. A cell’s type is its identity. The available cell types vary by species. Tadpoles possess 4 cell types (tan, red, black, and blue) while French Flags contain only three cell types (white, red, and blue) (Figure 2A). Voltage is another type of information that can be seen as a goal-encoded layer [26] in the form of different bioelectrical patterns (direct, indirect and the binary trigger, see Figure 2B).

Cellular Perception-Action Cycle

An activated cell behaves in a precise series of steps, similar to Algorithm 1 from Pio-Lopez et al. [85]. First, a bias of 0.5 and cellular attributes from times t and t_-1 are passed as inputs to the cell’s ANN, which then returns the outputs.

The inputs are energy and molecule levels, cell type, and voltage. In addition, French Flag cells have one extra input which is the current fitness of their stripe, similarly to [85].The outputs are the amount to open gap junctions and the number of molecules to distribute to neighboring cells (Figure 1A). These outputs are then processed by the cell’s communication module, similarly to Algorithm 2 from Pio-Lopez et al. [85]. The cell’s gap junctions are opened, and its molecules are transferred to the neighbors. The number of molecules distributed depends on both outputs. Molecules are distributed to neighboring cells in a randomized order. If the number of molecules a cell must distribute times the gap junction opening times the number of neighboring cells is greater than the cell’s molecule stock, then the cell will keep distributing molecules until it runs out of molecules. Then, the cell differentiates into a new cell type if it has met the new cell type’s threshold for molecules (Figure 2A). For details of various differentiation mechanisms explored in this work, look at (Supplementary Material, Figure S1). Every differentiation costs the cell 0.50 energy.

Cells have a metabolic cost for survival, between 0.65 and 0.95 energy per activation, depending on the difficulty of the learning task. This energy loss per activation is counteracted by an energy gain of between 0 and 1 energy per activation. This energy gain is in proportion to the organism’s fitness at the beginning of the current timestep. Energy acts as a reward incentive for cells to cooperate towards correct morphogenesis. In sum, cells rapidly lose energy and die unless the organism reaches close to the target morphology by adulthood. Thus, an organism-level anatomical homeostatic loop is implemented from the binding of cellular-level metabolic homeostatic loops (Figure 2B).

3.1.3. Development and Fitness

Each embryo underwent a 100 timestep developmental cycle into an adult. Organisms learned to differentiate their cells in a way that increases fitness. For every timestep in a developmental cycle, all cells are activated in a random order. At the end of the 100^th timestep, the adult’s fitness is evaluated out of 100 points. Fitness is quantified by how close the adult’s morphology is to its pre-defined target morphology (either the Tadpole or the French Flag). Each species has its own specific fitness evaluation. For tadpoles, they received 30 points if the eyes are in the correct quantity and position, 28 points for the mouth, 6 points for the nose, and 36 points for the skin. Skin was weighted less to enable other organs to form despite the prevalence of skin cells.

For the French Flag, organisms receive a fitness score according to how much they reached the appropriate morphology with 33.3 points for each stripe being in the correct quantity and position, each stripe matters equally.

3.1.4. Evolution and Genetic Algorithm

Each population of organisms underwent a 250-generation evolutionary cycle, with each generation composed of hundreds of developmental cycles (Figure 1C). A population contains 350 organisms. To create and evolve populations, we used the ES-HyperNEAT genetic algorithm [94]. HyperNEAT uses geometric regularities of the task domain to evolve large-scale neural networks that emulate neural connectivity patterns in the brain [95], but cells in non-neural tissue have similar connectivity. There are several stages of evolution. First, a population of 350 organisms was initialized, each with different ANN weights. For each generation over 250 generations, each organism receives a fitness equal to the average fitness from three development cycles. Averaging fitness controls for intra-organism fitness variability stemming from random activation. The top 20% most fit organisms reproduce to populate the next generation. These offspring are also vulnerable to random mutations. We set an overall mutation rate of 15%.

3.1.5. Bioelectric Patterns

Every cell of an organism is assigned a voltage integer. A bioelectric pattern is defined as the 2d-array of cell voltages that maps to an organism’s 2d-array of cells (Figure 1D). Organisms may decide to use the bioelectric pattern present to them, or not, to achieve their morphogenetic goals during evolution. We did not force or incentivize them to use their bioelectric patterns, in order to observe the evolutionary-developmental conditions necessary for the adoption of bioelectricity.

There are three types of bioelectric patterns in this work: direct, indirect, and the binary trigger. They are named based on their relationship with the target morphology. Direct bioelectric patterns have voltages that map one to one with cell types. This is not the case for indirect bioelectric patterns. Binary triggers have a 5 timestep rest period, followed by a 20 timestep “trigger” period, followed by a rest period for the remainder of the simulation. The binary trigger has a homogenous voltage.

Tadpole organisms only had access to a direct pattern during evolution. And there are two sets of French Flag organisms: (1) those who had access to four indirect patterns, or (2) those who had access to one binary trigger pattern during evolution. All bioelectric patterns used for the evolutions are in Figure 2B.

3.1.6. Parameters and Code Access

Genetic algorithm parameters used in this work are the same as in the previous work [85]. All code used in this work is freely available at https://github.com/sammyhansali/ScaleFreeCognition. Agent-based modeling framework, Mesa [96], was used to visualize organismal development along with dashboards for tracking fitness, gap junction connectivity, molecule exchange activity, energy, and more. The MultiNEAT package was used to implement ES-HyperNEAT [93].

3.2. Material and Method for the Planaria SSRI Experiment

A clonal strain of Dugesia japonica was kept in Poland Spring Water (Framingham, MA, USA), and maintained as in [97] and was starved for >7 days before cutting, a standard method used in this model system to reduce variability of data by controlling for the metabolic status of individuals [97]. To assess the impact of transient fluoxetine treatment on bioelectric states, middle third (trunk) fragments from D. japonica were allowed to regenerate in a 2uM Fluoxetine or Sertraline solution for 3 days, at which point the solution was washed out and the samples were placed in Poland Spring water for 1 week (at 13°C to prevent fissioning). The were then photographed and scored for the number of heads, and to identify any changes in head morphology (which were not observed). Scoring was performed under a Zeiss SV6 dissecting microscope.

4. Simulation Results for the 3 Types of Bioelectric Patterns to Reach a Target Morphology: Direct, Indirect, and the Binary Trigger

4.1. All Three Types of Bioelectrical Codes Allow the Reaching of Target Morphologies

To determine whether the artificial organisms learned to navigate the morphospace by following the bioelectric information available to them, we trained population of organisms per conditions (20 organisms for the direct bioelectrical pattern, 20 organisms by patterns for the indirect patterns organisms and 20 for the binary trigger organisms). Each species determines the embryonic morphology, target morphology, and cell types (Figure 2A) for the developmental task of that trial (tadpole for the bioelectrical pattern, French flag for the indirect and binary triggers). Each trial was repeated at least 20 times.

We found that at population level of the direct pattern organisms, the fitness is on average of 91 (σ=5%) demonstrating a good solution of the morphological task using direct pattern encoding in tadpole organisms (Table 1). In the representative direct pattern example shown in Figure 3A, the organism was able to reach the target morphology and achieved a fitness of 91.4 (σ=11.2%) over 20 developmental runs of the same organism. It displayed an initial burst of molecule exchange activity that set the stage for organogenesis. This was preceded by all of its gap junctions opening fully for one timestep, before closing off by 80% for the rest of the simulation. By time 50, the location for the future eyes, nose, and mouth were demarcated. By time 75 and 100, the resolution of the face increases. Activity plateaued over time as fitness increases and error decreases. Notably, the organisms reached a perfect fitness of 100, by time 200, demonstrating allostasis and stability capabilities seen in the previous similar models [98]. This is a common strategy for direct pattern organisms, especially with the gap junction opening fully for 1 timestep before closing off dramatically for the rest of the developmental period. Some direct organisms deviated from this strategy and instead learn to assign gap junction opening amounts specific to each cell type.

At the population level for the indirect bioelectrical pattern organisms, the fitness scores were: 83 (σ=3%) for “Half and Half” bioelectrical pattern, 77 (σ=3%) for “Cross”, 80 (σ=5%), for “Triangles, and 78 (σ=4%) for “LL” (Table 1 and Figure 2B for the indirect bioelectrical patterns). Overall, the morphological task is harder than with the direct bioelectrical pattern with an average fitness among all indirect patterns of 79.5 (Table 1). The representative indirect pattern example (Figure 3B), which utilized a “Triangles” pattern, was also able to complete the developmental task. It achieved a maximum fitness of 96% during development and a fitness of 80.9 (σ=5.2%) over 20 runs. From timestep one onwards, its gap junctions were permanently open and molecule exchange was frequent. By the time 25 the blue stripe is clearly demarcated. By 75, the French flag starts to become discernable as the red cells colonize the right stripe. Finally at step 100, following several spikes of activity, the target morphology is realized with an error of 7 cells in the red stripe.

Finally, at population level of the binary trigger organisms (Figure 4C), the fitness is on average of 81 (σ=5%) (Table 1). In the example displayed in Figure 3C, we found that it solved the French flag with a fitness of 80.4 (σ=4.2%) over 20 runs. It immediately focused on isolating the blue stripe and white stripe, usually achieving that by the end of the binary trigger, step 25. From then on, its strategy was more stochastic. The organism started to have a three-cell wide red stripe on the right side, but it was not successful at retaining the red stripe for a long period. But it was enough to reach 90 fitness in the developmental run shown in Figure 3C. The cells repeatedly differentiated from white to red and red to blue. This strategy essentially continued in the post-developmentary period, causing the decay of the tissue by step 200. Generally, there are two types of binary trigger organisms we observed. Those that began development before the binary trigger was presented at time 5, and those that began at time 5.

Thus, organisms with all three types of bioelectric patterns learned to navigate anatomical space to develop from early-stage embryo to the adult target morphology using their bioelectrical pattern. However, the task is much harder for the indirect patterns organisms and the binary trigger compared to the direct pattern organisms given their fitness difference, probably as learning a direct mapping is a linear task easier than learning an indirect one, with a non-linear relationship.

4.2. Regulative Morphogenesis Depends Fully on the Direct Pattern, Partially for the Indirect Pattern, and Depends on the Duration of the Pattern for the Binary Trigger

The digital organisms can reach their anatomical goals via high-level bioelectric patterning information. However, we did not specifically reward for utilization of bioelectricity data because we wanted to understand the developmental conditions which may have caused cellular collectives, over the course of evolution, to co-opt spatial gradients of membrane voltages for the purpose of morphogenesis. Without forcing the adoption of bioelectricity, would artificial organisms learn on their own to utilize bioelectric patterns to traverse anatomical space, or rely on the other ANN inputs (energy, molecules, cell type)?

In order to answer this question, we devised the “Resetting Bioelectricity” experiment (Supplementary Material, Figure S2). We reset the native bioelectric pattern of artificial embryos to a homogenously hyperpolarized bioelectric pattern (in every step, except for the binary trigger, where each cell has a membrane voltage of -100 mV, at timestep zero. This was done half-way during a simulation as a loss-of-function experiment to determine whether or not they were relying on bioelectrical mechanisms for their effective morphogenetic performance. This pattern was selected because it shares no resemblance to any of the developmental bioelectric patterns used for direct and indirect pattern organisms (Figure 2B). Therefore, if an organism relied on its native bioelectric pattern for correct anatomical development, this novel pattern should disrupt that. This “Resetting Bioelectricity” experiment was modified for binary trigger organisms, as it is already a homogeneous pattern, in order to understand their temporal dependence on the binary trigger duration, to know if they are really dependent on a presentation of the binary trigger to activate the downstream pathways to resolve the morphogenetic task. The native bioelectric trigger lasts for 20 timesteps during development. So, for each organism we reset their bioelectric pattern for 1 timestep, 2 timesteps, all the way up to 20 timesteps. Then we observed how changing the duration of the binary trigger affected fitness (Supplemental Material, Figure S3).

We found that resetting bioelectricity inhibited correct development, but the impact varies with the different types of pattern. Indeed, by comparing reset scores between direct and indirect organisms, the direct organisms have population scores of 0 (σ=0) while indirect organisms don’t fail completely on the morphogenetic task: they have higher scores but still low in terms of development and with a lot of variance (population fitness score of 31 (σ=24%) for “Half and Half” pattern, 51 (σ=16%) for “Cross”, 48 (σ=15%) for “Triangles and 52 (σ=17%) for “LL” (Table 1). The high variance in resetting scores for indirect organisms supports the fact that different indirect organisms have different levels of dependence on bioelectricity for morphogenesis. Therefore, organisms with the direct patterns (Figure 4D) are fully dependent on bioelectricity for survival, but indirect types have less dependence.

We also found that binary trigger organisms varied in their bioelectrical dependance (see Figure 4F and Supplemental Material, Figure S3B-C). Some required only 1 timestep of their binary trigger in order to reach a fitness comparable to if they had the full 20 timestep trigger. For others, 5, 10, or 15 timesteps is necessary. Some did not require the bioelectric pattern at all to reach the target morphology. This is shown in Figure S3A, where resetting the whole binary trigger does not affect 60% of the organisms, while the rest either die or develop incorrectly. At the population level, when removing the entire binary trigger (all 20 timesteps) the fitness is on average of 58 (σ=27%). This demonstrates how resetting the full binary trigger can severely inhibits the ability for binary trigger organisms to resolve the French Flag but the high variance shows too that some organisms in the population don’t rely fully on it to resolve the morphogenetic task (Table 1). The distribution of fitness scores for binary trigger organisms undergoing full resetting (Figure 4F) shows a bimodal distribution. The majority of the population didn’t survive the resetting of the pattern, while the others can reach up to fitness=80, showing a low dependency on the bioelectrical pattern.

Direct, indirect, and binary trigger pattern organism examples shown in Figure 5A, 5B, and 5C, respectively, achieved fitness of 0 (σ=0%), 61.4 (σ=3.84%), and 43.6 (σ=2.29%) over 20 runs after resetting bioelectricity. Given that these specific organisms have baseline fitness of 98, 83%, and 76% without any perturbation, respectively, these organisms rely on their native bioelectric patterns for development. Generally, what occurs for these types of organisms that have their native patterns reset is the following. First, they show rapid cellular differentiation. This is what is seen in practice – cell collectives typically undergo global, rapid, and stochastic cellular differentiation. Since cellular differentiation has a metabolic cost in our model, this accelerates the demise of the organism in a positive feedback loop.

The mapping between indirect and binary trigger patterns with their downstream morphologies is complex, as opposed to the linear mapping of direct organisms. As shown in the population statistics of indirect and binary trigger organisms, some organisms with non-linear goal-encoded pattern (indirect and binary trigger) learn to rely less on their membrane voltages for development and may rely on other things, like local information or the starting state of the tissue, or a memory of the goal encoded in the neural networks of the cells.

Therefore, we conclude that organisms with direct patterns are fully dependent on the bioelectric pattern (they perished immediately upon removal), while those with indirect patterns and binary triggers can evolve to different degrees of robustness to loss of bioelectrical information, to differing extents based on the specific pattern and the duration of the binary trigger.

4.3. An Emergent Robustness to Changes in Initial States Configurations for the Direct Pattern

Previously we showed that artificial organisms rely on bioelectricity to reach their morphogenetic setpoints (to various degrees depending on their type of bioelectric pattern) but how robust are these developmental processes to external perturbations? Indeed, did organisms learn to map the native bioelectric pattern to the target morphology, so that they may reach it from various initial states? Do they show developmental robustness [98,99]?

To address this question, we devised the “Morphological Robustness” experiment (Supplementary Material, Figure S4), in which embryos were tasked with developing into the correct target morphology despite starting off in different morphological states. We created 20 new starting morphologies for each organism (Supplementary Material, Figure S5) with direct and indirect bioelectrical patterns. Each organism was evaluated on these 20 developmental cycles for their morphological robustness. We didn’t apply this experiment to the organisms with a binary trigger as we want to understand the mapping between the goal-encoded bioelectrical pattern and the morphological initial condition, and the binary trigger is more a higher level bioelectrically encoded switch activating the memory of the cellular network to reach a specific morphology.

We discovered that at the population level, the direct pattern organisms were morphologically robust, by reaching 65 (σ=16%) of fitness score on average, with the 20 new starting morphologies (Table 1). We also observed that several organisms even reached more than 90 of fitness score after perturbation (Figure 4G). The representative example of the organism with direct pattern (see Figure 6A) scored a fitness of 97.7 (σ=3.00%) over 20 runs. The target morphology was almost complete by time 75, which is 25 timesteps earlier than the baseline (Figure 3A). The common process by which the organisms reach the anatomical setpoint from a new starting point is identical to the framework laid out in (Section 3.1).

Indirect organisms were not morphologically robust (Figure 4H), reaching at population level a fitness score of 20 (σ=12%) for “Half and Half” indirect pattern, 6 (σ=11%) for “Cross”’, 2 (σ=4%) for “Triangles” and 1 (σ=2%) for “LL” for a general average of 7.25 of fitness score for organisms with indirect pattern (Table 1). The representative example of the indirect pattern organisms, utilizing the “Half and Half” pattern and a starting morphology biased towards white and blue cells, achieved a fitness of 31.7 (σ=7.72%) over 20 runs (Figure 6B). By around timestep 25, the organism reached its peak fitness level and the morphology remained mostly unchanged until death. No coherent stripe was formed. Activity ceased and the organism died after timestep 110.

Thus, we conclude that direct pattern organisms show emergent morphological robustness to external morphological reconfiguration as a byproduct of their linear mappings between voltage patterning and target anatomy. Indirect pattern organisms displayed almost non-existent morphological robustness to changes in initial states configurations in our experiment.

4.4. Emergent Robustness for the Organisms with Indirect Patterns to Bioelectrical Perturbation

We investigated earlier how completely resetting the bioelectric pattern inhibits correct development of the organisms with direct patterns and with different extent with organisms with indirect pattern. Would it be possible that organisms rely only on specific part of the indirect bioelectrical pattern? What about smaller perturbation, such as resetting half of the pattern or adding artificial noise to the whole pattern? To address this question, we devised the “Bioelectric Robustness” experiments (Supplementary Material, Figure S6). We focused this experiment solely on organisms with indirect patterns, due to their lack of morphological robustness and the possibility that the key to their robustness lies in the nonlinear dependencies they learned from their bioelectric patterns. Also, the direct bioelectrical pattern encodes a direct mapping with the goal to reach in the morphospace and would not be expected to be robust to its alteration.

We evaluated the resilience of the indirect pattern organisms to bioelectric perturbations by either: i) resetting half of their bioelectrical patterns (by hyperpolarizing them to -100mV, similar to section 3.2) ii) introducing noise (by polarizing cells by a random choice of -20, -10, -5, 0, 5, 10, or 20mV). Each organism was evaluated on these two developmental cycles for their bioelectric robustness.

We found that there is large variance at population level, which has “Cut” robustness fitness scores of 31 (σ=31%) for indirect pattern “Half and Half”, 42 (σ=29%) for “Cross”’, 43 (σ=32%) for “Triangles” and 40 (σ=32%) for “LL” for an average fitness of 39 of fitness score for organisms with indirect pattern (Table 1). “Noise” population robustness fitness scores also had high variance, with a score of 47 (σ=26) for “Half and Half”, 43 (σ=24%) for “Cross”, 41 (σ=27) for “Triangles”, and 44 (σ=27%) for “LL” for an average fitness of 43.75. It appears that, the population of indirect organisms display a bimodal distribution for fitness scores under the “Cut” and “Noise” robustness experiments (Figure 4I and Figure 4J). We have two groups in each condition, one composed of organisms that didn’t survive the bioelectrical perturbation and another that presented emergent robustness during evolution.

The representative example of the survival group of the organisms with the “Triangles” indirect pattern (Figure 7B) achieved a fitness of 82.9 (σ=7.16%) and a fitness of 84.4 (σ=6.89%) over 20 runs for the cut and noise experiments, respectively. “Half and Half” organisms outperformed other indirect populations in Noise Robustness, but underperformed in Cut Robustness. The former may be due to the relatively low structural complexity of “Half and Half” compared to the other indirect patterns. The latter is easy to understand, because the “Half and Half” contains half of its information on each side, since both left and right halves have different polarization levels. Therefore, we conclude that indirect bioelectric patterns have emergent robustness to bioelectric perturbations.

4.5. An Emergent Generalizability Competency to New Bioelectrical Pattern for the Direct and Indirect Patterns Organisms

Animals can be anatomically altered by reprogramming their developmental bioelectric software [1], such as: regenerating planaria induced to grow two heads or no heads [100] or heads of other species [56,65], or X. laevis developing functional ectopic eyes in aberrant locations [101]. We wanted to know if we could reprogram organisms to reach new target morphologies by simply changing their native bioelectric pattern at initialization. Can they generalize past the relationship between their native bioelectrical pattern and target morphologies?

To answer that question, we devised the “Generalizability” experiment (Supplementary Material, Figure S7). For each organism, we created 3 new bioelectric patterns by rotating their native pattern by 90°, 180°, and 270° (relative to their initial orientation). Each organism was initialized with these novel patterns after they have been evolved with the standard parameters and patterns and then evaluated for their generalizability.

At the population level, direct pattern organisms reached a fitness score of 83 (σ=9%). While indirect organisms were less successful with generalization, they reached at population level a fitness score of 57 (σ=20%) for indirect pattern “Half and Half”, 27 (σ=7%) for “Cross”’, 24 (σ=2%) for “Triangles” and 23 (σ=3%) for “LL” for a general average of 32.75 of fitness score for organisms with indirect pattern (Table 1). This difference in performance is also reflected in the population level fitness distributions for generalizability, shown in Figure 4K and Figure 4L. Direct organisms all survived development (Figure 4K), while there seems to be two classes of indirect organisms: those that die (the majority of the organisms) and those that are moderate to high for generalizing and are successful at development, mainly with ‘Half and Half” pattern (Figure 4L).

The representative examples of the direct pattern organism and indirect pattern organism (utilizing a “Half and Half” pattern) shown in Figure 8A and 8B, have fitness scores of 92.8 (σ=6.93%) and of 79.4 (σ=16.1%) over 20 runs, respectively. The strategy that organisms followed was almost the same as their original strategies in standard conditions.

Thus, we observed that the organisms learned a larger mapping between the bioelectrical patterns and the target morphology. Through experiments comparing normal development to outcomes after bioelectric pattern rotation, we found that organisms with direct and indirect bioelectrical patterns adapted well and show emergent generalizability competency to new bioelectrical patterns.

4.6. An Emergent Repatterning Capabilities for Direct Pattern Organisms in Post-Developmental Phase

We showed in the previous section that artificial embryos could be reprogrammed to generalize to novel target morphologies that had never been learned in their developmental cycles. A similar biological phenomenon involving morphological remodeling is that of repatterning [102]. This involves the reorganization and reconfiguring of biological tissue to achieve a new outcome, such as the remodeling of adult epibranchials of metamorphosing species of lungless salamanders during development [103] or wound healing and blastema formation during Drosophila larvae imaginal disk regeneration [104]. In the context of our simulation, we define repatterning as the ability for an organism, after the conclusion of its 100 timestep developmental period, to remodel its cells to reach a new morphological state given a novel bioelectric pattern. Previous work [85] showed that artificial organisms are capable of morphological allostasis in this post-developmental period. This ability to maintain anatomical stability with dynamic responses after development was not directly programmed in by the modeler and is therefore emergent. If artificial organisms during the post-developmental period have the emergent ability to maintain their form, do they also have the emergent ability to create new forms (repatterning)? Can we induce them to repattern themselves into new morphologies by only changing their bioelectric patterns during the developmental phase?

To answer that question, we devised the “Repatterning” experiment for direct and indirect patterns organisms (Supplementary Material, Figure S8). This experiment is identical to the “Generalizability” experiment (Section 3.4) except each organism was allowed to develop normally until timestep 100 when these novel patterns replaced the native ones. The energy levels for these organisms were reset back to maximum at the start of the repatterning period as for most of the organisms were exhausted after development. They were then evaluated at timestep 200 for their repatterning ability.

We found that at populations level, for direct bioelectrical pattern organisms, their fitness was 22% (σ=33%), a stark drop from the 83 (σ=9%) population average for generalizability (Table 1). The high variance in outcomes for organisms repatterning means that some organisms start to repattern, while others cannot. Indirect organisms also had low repatterning capabilities, with an average fitness of 13 (σ=14%) for the entire population of Half/Half indirect pattern organisms, 7 (σ=11%) for “Cross”’, 1 (σ=3%) for “Triangles” and 0 (σ=0%) for “LL” for an general average of 5.25 of fitness score for all the organisms with indirect pattern (Table 1). Looking at the population level fitness histograms, we found that even if the entire direct pattern population tended to have low fitness scores, they contained a few organisms that showed some level of competency at repatterning (as displayed in the histogram, Figure 4M), such as in Figure 9A. However, overall, indirect organisms performed poorly with no organisms that can reach a fitness score over 50 after 100 steps, meaning a 0% survival rate (Figure 4N).

We observed that successful repatterning ability for direct pattern organisms follow the general behavior of previous experiments with direct pattern organisms typically more capable than the indirect pattern ones. The direct pattern organism and indirect pattern organism (utilizing a “Half and Half” pattern) examples are shown in Figure 9A and 9B, scoring a 200-timestep fitness of 83.0 (σ=7.02%) and of 54.3 (σ=2.96%) over 20 runs, respectively. Both patterns are rotated 90 degrees counterclockwise at time 100, and fitness was now evaluated at time 200 instead of 100.

Thus, we can conclude that artificial organisms with direct patterns have the emergent ability to repattern themselves after evolutionary training and alter their adult morphology in response to a changed bioelectric pattern. Conversely, those with indirect patterns don’t exhibit any repatterning capability (Figure 4N).

5. Selective Serotonin Reuptake Inhibitors (SSRI) Simulation and Experimental Results

Selective serotonin reuptake inhibitors (SSRIs) are a class of drugs prescribed for cases of depression and anxiety [105,106]. They mainly target a component of the serotonin signaling system, the serotonin transporter, SERT. Serotonin signaling has previously been implicated as a component of bioelectric circuits [105]. Some have suggested that SSRIs target bioelectricity in the context of developmental defects [105,107,108]. SSRIs act on the serotonin signaling system, which has been suggested to be involved in processing rewards and punishments in the brain [109,110]. SSRIs can also function to change membrane voltage of neural and non-neural cells [105,111,112].

As found above, artificial organisms learned a mapping between target morphology and their native bioelectric pattern that is robust and amenable to generalization and repatterning. Does that mapping start to break down under the influence of SSRI-like conditions? Does that lead to morphological defects in artificial organisms as seen in biological organisms [108]? To answer these questions, we simulated SSRI exposure to explore its effects on the development of artificial organisms and made predictions we verified experimentally.

5.1. Simulated SSRI Induced Loss of Regenerative Precision and a Bistable Morphogenetic Process

In order to simulate an SSRI exposure experiment, we started by evolving organisms that would pursue the French Flag morphology. These indirect pattern organisms have a homogenously depolarized bioelectric pattern (as shown in Figure 10) instead of the bioelectric patterns used for the other direct organisms (as shown in Figure 2B). This new pattern resembles the “Rest” pattern shown for binary trigger organisms (Figure 2B). Even so, these organisms are not binary trigger organisms because they are only exposed to one pattern for their entire developmental period (as opposed to binary trigger organisms which alternate between two patterns).

For those reasons, we selected a homogenously polarized bioelectric pattern (which is a type of indirect pattern). Next, we selected 20 of these organisms that had a fitness of 90 or more on the given task. At the population level, the fitness for this group was 92.4 (σ=2.1%) (Supplemental Material, Figure S9A). They also showed dependency on bioelectricity (Section 3.2). Then, we altered their cellular reward machinery by eliminating energy loss at each step similarly to SSRI exposure similarly to [109,110]. Thus, cells may no longer gain top-down instructive information on morphological error. We also perturbed the voltage of 1, 2, 5, and 10 random cells at once, by hyperpolarizing them anywhere from -5, -10, or -30 mV, simulating the weakening of cells’ abilities to correctly interpret the organismal bioelectric pattern under SSRI [105,111,112].

This simulated experiment led to two discoveries: global morphological degradation and a bistable development. Global morphological degradation, in a biological context, can be thought of as the inability for a cellular collective to determine how to differentiate to reach the target morphology. This state is characterized by a significant deviation from the target morphology and frequent random, undirected differentiation of the cells, or a loss of precision in organisms’ normally quite tight ability to reach the same, invariant pattern. Multistability in dynamical systems theory is “the ability of a system to achieve multiple steady states under the same external conditions” [113]. In a developmental biology context, this would equate to an organism that can become one of multiple target anatomies, under the same environment. There are many animals that exhibit this behavior under different contexts, such as planaria during regeneration [114].

5.1.1. SSRI Exposure Leads to Global Morphological

By simulating SSRI exposure as defined above with the alteration of the cellular reward mechanisms (via removing energy loss) and bioelectric states (via randomized voltage perturbations of cells), we found that the number of cellular voltage perturbations was the key factor that determined whether global morphological degradation (loss of morphogenetic precision) occurred or not (Figure 10). Hyperpolarizing only one cell did little to affect development. However, hyperpolarizing two cells was enough to begin to see a disruption of anatomical homeostasis. Development would typically proceed along a similar trajectory as normal, however there would be a slight increase in rapid “back-and-forth” cellular differentiations. Sometimes, the target morphology would be worse due to misplaced cells and incomplete structures. Moreover, the global morphological degradation became most apparent when five or more perturbations were introduced to the SSRI-exposed organism. Depending on their placement, a general weakening of anatomical homeostasis could occur (Figure 10C), or the birth of a bistable developmental process (Figure 10B).

A factor that impacted the occurrence of global morphological degradation was the location of perturbations. As discussed in depth in the next section, on rare occasions the specific cells that become hyperpolarized can be the difference between a benign perturbation that has seemingly no effect on morphogenesis, and an alternate developmental pathway. That factor usually impacts degradation in concert with another factor, the order of cellular activations during development. Interestingly, an insignificant factor was the magnitude of the cellular hyperpolarizations. For any given organism, hyperpolarizing an identical set of cells by -5, -10, or -30 mV had no distinguishably different effect. This could be because the organism and its ancestors were only exposed to one voltage value (0 mV) throughout their lifecycles, thus they never learned to distinguish between different voltage states.

Based off these observations, we designed a test to quantify morphological degradation by comparing the fitness of these 20 organisms before and after SSRI exposure with perturbations (Supplemental Material, Figure S9). In this experiment, we exposed each organism to SSRI and perturbed 5 random cells, as described above. We did this three times for each organism and averaged the fitness values to get a mean fitness score for each organism. At the population level, the fitness for these organisms dropped from 92.4 (σ=2.1%) to 54.8 (σ=13.3%) (Supplemental Material, Figure S9).

By altering cellular reward mechanisms (via removing energy loss) and bioelectric states (via randomized voltage perturbations of cells) similarly to SSRI exposure, we induced widespread morphological degradation in our model organisms similarly to results reported in biological experiments [108].

5.1.2. SSRI Exposure Leads to a Bistable Developmental Process

We also found that SSRI simulated exposure leads to leads to a multistable developmental process in a one representative organism (Figure 10). Under normal conditions, this organism reaches the French Flag morphology with a fitness of 74.5 (σ=15.1%), as shown in Figure 10A. When this organism was subjected to random perturbations of 1 cell at a time, no notable degradation or change in stability ever occurred. When multiple cells were perturbed at random, there was noticeable degradation but no multistability. However, we observed that when only the cells of the left stripe were perturbed, regardless of the amount of hyperpolarization for each cell or the number of cells that were hyperpolarized, the organism exhibited a bistable developmental process. Over 100 runs on the configuration shown in 10B, this organism had a 67% probability of developing into the wild type expected morphology, with the remaining 33% of times resulting in an entirely distinct form (Figure 10B). In this bistable mode, the organism had a fitness of 91.4 (σ=11.5%). Interestingly, when it solved the French Flag it always did so with near perfection, getting 98-100 fitness. When it pursued the mutant type morphology, the fitness was closer to 55-66. This new mutant type of morphology had no red stripe, a central 1-cell thick blue and red stripe at time 100, and a central red dot or dots at time 200 that replaced the blue stripe.

This bistable phenomenon quickly disappeared as soon as a cell was perturbed outside of the left stripe area (Figure 10C). Once again, the number of perturbations and degree of hyperpolarization does not matter – as soon as a cell outside of the blue stripe is perturbed on the bioelectrical pattern, there is morphological degradation. This leads to a French flag-like morphology, except it is typically missing its red stripe and with many actively differentiating cells in the red stripe area (Figure 10C). The fitness dropped to 63.6 (σ=0.82%) when perturbing outside of the left stripe. As indicated by the small deviation, this effect reliably reproduced the same morphologically decoherent morphology shown in Figure 10C.

We can conclude that this bistable phenomenon indicates that SSRI-like exposure paired with specific bioelectric modifications can induce a bifurcation in developmental pathways, leading to alternative stable morphologies. This connects to work done in tadpoles showing how depolarization of as little as three cells could cause an alternate, large-scale developmental pathway to emerge [115]. The complexity of bioelectric influences on development suggests that SSRI exposure by disrupting the top-down information related to the developmental reward may unlock latent developmental trajectories, resulting in divergent morphological outcomes under identical conditions. This finding is surprising as changing only a few cells via induced hyperpolarization causes a multi-stable developmental trajectory to emerge. Therefore, we decided to test experimentally this prediction.

5.2. SSRI Experimental Results: Regenerating Planaria in Fluoxetine or Sertraline

Our model made a prediction about the loss of morphogenetic precision resulting from loss of stress-induced motivation driving the anatomical homeostatic loop. We sought to test this in vivo, using a widely available class of drugs that are used for reducing anxiety in humans and other animals [116,117,118,119,120,121,122] and that has previously been tested in planaria in toxicology [123,124] and parasitology [125] contexts: the serotonin reuptake inhibitors Fluoxetine and Sertraline. Specifically, these SSRIs have been studied with respect to their effects on reducing behavioral motivation and the drive for specific outcomes [126,127,128]. Applying either drug to middle fragments of planarians (Figure 11A) resulted in a mix of 0-, 1-, and 2-headed regenerates (Figure 11B, quantified in 11C). Compared to exclusively 1-headed animals in control conditions (n>500), cohorts treated with Sertraline or Fluoxetine exhibited respectively 31% and 17% 0-headed animals, and 4% and 8% 2-headed animals (n=90, n=155 respectively). This effect was significant to p<0.01 (2.79×10^-36), from which we conclude that serotonin reuptake inhibitors reduce the precision (but do not induce one consistent anatomical fate) of regeneration.

6. Discussion

This study utilized evolutionary simulations to explore the role of bioelectric signals in guiding morphological outcomes during development. Although it is known that bioelectric patterns are crucial for pattern formation and can be classified into three types—direct encoding, indirect encoding, and binary trigger—there is a lack of computational tools to understand how these patterns evolve and confer advantages. To address these gaps, we developed a neural cellular automaton (NCA) and used evolutionary algorithms to optimize these models for reliable, regulative morphogenesis driven by the three types of bioelectric patterns. Our findings are as follows, we found that: (1) All types of bioelectrical codes allow the reaching of target morphologies (2) Morphogenetic behavior depends on the bioelectrical patterns and the duration of the binary trigger (3) An emergent morphological robustness for the direct pattern (4) A relative emergent bioelectric robustness for the organisms with indirect patterns to bioelectrical perturbation (5) An emergent generalizability competency to new bioelectrical pattern for the direct and indirect patterns organisms (6) Emergent repatterning competency for the direct pattern organisms in post-developmental-phase. We also developed a selective serotonin reuptake inhibitors simulation, which diminished the ability of artificial embryos to correctly interpret bioelectric patterns due to an altered organismal reward machinery, leading to higher variance of developmental outcomes (loss of precision), global morphological degradation, and induced in some organisms anatomical bistability. These computational findings have been validated by data collected from planaria SSRI exposure. Additionally, we simulated the effects of selective serotonin reuptake inhibitors (SSRIs) on artificial embryos, observing that these inhibitors impaired the correct interpretation of bioelectric patterns. This resulted in increased variance in developmental outcomes, global morphological degradation, and bistability in the target anatomies in some artificial organisms. This computational prediction was validated by experimental data from planaria exposed to SSRIs.

In the context of somatic pattern memories and bioelectric circuits underlying regeneration, propose that bistability enables stochastic outcomes in the restoration of correct large-scale anatomies from diverse initial states. This is achieved through the exploitation of bioelectric encoding of distributed goal states, akin to how the brain utilizes memory to achieve specific outcomes [114]. The concept of bistability is not only relevant in the context of regeneration and developmental biology but also extends to neuroscience. It suggests that bistability plays a role in cognitive processes and memory formation, drawing parallels between bioelectric circuits in regeneration and neural networks involved in memory, representation, and perceptual bistability [129]. This highlights the interconnectedness of bistability across different biological systems, emphasizing its significance in understanding complex biological phenomena. Moreover, the utilization of bistability in bioelectric circuits underscores the dynamic nature of cellular processes, similarly to neural processes, where stability and plasticity coexist to drive cellular differentiation and specialized functions in multicellular organisms [130].

We found several differences between organisms with direct and indirect bioelectrical pattern. Direct pattern organisms have a good morphological robustness to changes in the initial conditions while indirect pattern organisms in our simulation didn’t develop this capability fully. However, these organisms showed an emergent robustness to bioelectric perturbations while direct pattern organisms by definition don’t have this competency as the direct mapping between bioelectricity and the anatomy implies that any change in the bioelectrical pattern will result in a developmental defect. This is maybe why evolution co-opted indirect bioelectrical patterns. Indeed, organisms solely relying on direct pattern are disrupted more easily as organisms learned a direct mapping, any change in the goal-encoded direct pattern would lead to developmental defect, whereas indirect pattern may allow for the encoding of redundancy of the information necessary to reach the morphological goal. Evolutionarily, indirect pattern organisms may have a better fitness overall as they have such a robustness. In addition, learning an indirect mapping between morphology and an indirect bioelectrical pattern may allow the exploration of a larger morphospace than the one encoded in a direct one. It allows for more adaptability as an indirect mapping can represent more morphologies than a direct one.

Evolutionarily, organisms with indirect bioelectrical patterns likely possess greater overall fitness due to their robustness. Moreover, learning an indirect mapping between morphology and bioelectrical patterns enables the exploration of a larger morphospace compared to direct mappings. This greater adaptability stems from the ability of indirect mappings to represent a broader range of morphologies or compressed morphological information. Our simulated experimental comparisons between normal development and outcomes following bioelectric pattern rotation indicate that both direct and indirect pattern organisms adapt well and exhibit emergent generalizability to new bioelectrical patterns. This competency may facilitate the introduction of new patterns in evolution, enhancing adaptability and overall fitness.

Our model makes a specific prediction for loss-of-function targeting of motivation pathways in regenerative morphogenesis, which we decided to test in vivo, via the serotonergic route. The serotonin neurotransmitter system exists in planaria. An ortholog of tryptophan hydroxylase, which catalyzes the rate-limiting step in serotonin synthesis, has been identified in D. japonica [131], planarian serotonin receptors have also been characterized [132,133], and several studies have found that SSRI treatment impacts planaria in numerous ways including behavior, DNA damage, and regeneration [124,125,134,135].

We found that two different serotonin reuptake inhibitors induced a “randomizing” effect in planarian regeneration, in which the normally extremely reliable regeneration of 1-headed animals from fragments was driven to sometimes produce 2-headed and sometimes 0-headed animals. This is unusual because conventional treatments, such as inhibition of the Wnt signaling pathway, typically reliably produce inhibition of head character or excess anteriorization [136,137,138].

The ability to reduce morphogenetic precision, and thus create a stochastic phenotype where both extremes are represented, instead of pushing the system into one specific direction along the anterior-posterior decision axis, is notable. That, and the fact that the phenotypes featured apparently normal heads with no tissue-level defects, are consistent with our contention that the stress-based system is regulating global information processing and error minimization, not a specific molecular- or cell-level endpoint. Furthermore, the ability of reagents used in human patients to reduce anxiety to affect morphogenesis in a way predicted by our computational model further supports a research program in which tools and concepts from cognitive science are used beyond brains, to address the complex navigation of cell collectives in their anatomical problem space [7,139].

Our study had a number of limitations, which can be overcome in future work. First, because we wanted to isolate and understand the evolutionary properties of one specific mechanism, we focused on bioelectricity; numerous other developmental mechanisms – including biomechanics, chemical gradients, etc. – could be included to determine whether and how these influence the evolutionary process. Also, in our simulations, the indirect pattern organisms didn’t show the emergent competency to be robust to changes to initial morphological conditions (see ‘Morphological robustness’ experiment). This is in contrast with many systems in biology, which show high adaptability to different initial conditions The task of learning the indirect mapping is harder than learning a direct one, demonstrating here likely the limits of our leaning algorithm, even if we found in some cases some outliers with organisms that could show an appropriate development but only to specific changes in initial conditions. Our system also had a static bioelectrical pattern, which is a simplification of highly dynamical bioelectrical patterns for development that will be implemented in future iterations of our model.

One purpose of this kind of simulation is investigating the computational role and encoding of bioelectrical patterns in morphogenesis, with the aim to uncover novel mechanisms of development and evolution. This knowledge is fundamental, not only for basic evolutionary developmental biology but also for biomedicine and synthetic bioengineering by allowing the development of a rational control over large-scale growth and form and potentially paving the way for innovative therapeutic strategies in regenerative medicine and bioengineering targeting anatomical setpoints [35,140].