Closing the loop on morphogenesis: a mathematical model of morphogenesis by closed-loop reaction-diffusion

Morphogenesis, the establishment and repair of emergent complex anatomy by groups of cells, is a fascinating and biomedically-relevant problem. One of its most fascinating aspects is that a developing embryo can reliably recover from disturbances, such as splitting into twins. While this reliability implies some type of goal-seeking error minimization over a morphogenic field, there are many gaps with respect to detailed, constructive models of such a process. A common way to achieve reliability is negative feedback, which requires characterizing the existing body shape to create an error signal–but measuring properties of a shape may not be simple. We show how cells communicating in a wave-like pattern could analyze properties of the current body shape. We then describe a closed-loop negative-feedback system for creating reaction-diffusion (RD) patterns with high reliability. Specifically, we use a wave to count the number of peaks in a RD pattern, letting us use a negative-feedback controller to create a pattern with N repetitions, where N can be altered over a wide range. Furthermore, the individual repetitions of the RD pattern can be easily stretched or shrunk under genetic control to create, e.g., some morphological features larger than others. This work contributes to the exciting effort of understanding design principles of morphological computation, which can be used to understand evolved developmental mechanisms, manipulate them in regenerative-medicine settings, or engineer novel synthetic morphology constructs with desired robust behavior.

The generation of complex form during embryonic development, and its repair and remodeling during regeneration, highlight fundamental problems that range from cell and evolutionary biology to control theory and basal cognition (Pezzulo and Levin, 2015; Pezzulo and Levin, 2016; Harris, 2018; Pezzulo, 2020). How can collections of cells cooperate to reliably produce the same species-specific target morphology? Moreover, what mechanisms enable them to robustly do so despite various perturbations? For example, planarian flatworms regenerate their entire body from large or small fragments of any type (Cebrià et al., 2018), while amphibian embryos maintain the right proportions even when many cells are missing (Cooke, 1979; Cooke, 1981) or made too large (Fankhauser, 1945a; Fankhauser, 1945b). This homeostatic property of multicellular morphogenesis has numerous implications beyond basic science, as it represents an attractive target for regenerative medicine and synthetic bioengineering approaches that seek efficient methods for the control of growth and form. A number of mathematical frameworks have been developed to help understand, predict, and control the decision-making of cells in the morphogenetic problem space. Here, we first review several popular approaches to modeling this process, highlighting their positive features and limitations. We then propose a new model which has interesting and useful features for quantitative modeling of morphogenesis.

Negative feedback is a very common and effective means of achieving reliability in nature (Alon, 2019). We will consider negative-feedback systems to arrive at a target body shape. We consider, as an example, the specific question of how morphogenesis creates bodies with five toes rather than four or six. In mice, this has been shown to be accomplished with a five-peaked reaction-diffusion pattern (Raspopovic et al., 2014).

More generally, reaction-diffusion (RD) and positional information (PI) are perhaps the two best known hypotheses in the field of morphogenesis. Green (Green and Sharpe, 2015) gives an excellent summary of both hypotheses as well as contrasting the two. RD (Turing, 1953; Gierer and Meinhardt, 1972; Kondo and Miura, 2010; Marcon and Sharpe, 2012; Meinhardt, 2012; Green and Sharpe, 2015; Painter et al., 2021) was proposed by Turing in 1952. In its simplest form, it uses two chemical species, A and I. A (the “activator”) generates more of A and/or I via chemical reactions; I (the “inhibitor”) similarly reduces their concentrations. Surprisingly, combining these reactions with the diffusion of both A and I can, in many cases, amplify small random concentration gradients into definite and striking patterns (see (Kondo and Miura, 2010) for many examples).

Intuitively, the activator A promotes more of both A and I. Thus, any small excess of A at any location quickly grows by positive feedback. Of course, [I] also grows at the same location; but I is assumed to diffuse faster than A, so there is soon relatively little of I at this peak, and so the peak stays a peak. The I near the peak prevents new peaks from forming until you get far enough away for [I] to drop, at which point the pattern repeats. This concept, local self-activation with long-range inhibition, has been the basis of most RD systems [though new versions have also been discovered (Marcon et al., 2016; Landge et al., 2020)]. All of the variants have the basic ability to start with small, random concentration variations and amplify them into stable large-scale patterns.

Almost 20 years after Turing, Lewis Wolpert published his positional-information hypothesis (Wolpert, 1969). It is attractively simple. First, some unspecified process creates a gradient of a morphogen from, say, head to tail. Next, cells use a gene regulatory network (GRN) to determine their position by sampling the morphogen gradient, and then differentiate accordingly. PI gained rapid popularity. But it never per se explained where the initial gradient came from. Furthermore, most morphogen gradients exhibit exponential decay, which implies that much of the field will contain very low concentrations. This would make it difficult (Lander, 2007) for a GRN to determine spatial locations in those areas.

RD is not inherently scalable. RD patterns have a characteristic length λRD, typically given by λRD=DAKD,A (where D A is the diffusion constant of A in m2/sec and K D,A is the degradation constant of A in sec−1). In a field of length L, an RD pattern typically repeats L/λRD times. Thus, longer fields typically result in more pattern repetitions. This was originally seen as an argument against RD (Marcon and Sharpe, 2012); while it is reasonable for larger animals to have, say, more spots, we would not expect a larger embryo to have extra fingers or toes. This objection was eventually partially overcome. Gierer and Meinhardt first proposed (Gierer and Meinhardt, 1972) a scale-independent version of RD. The advent of modern molecular-biology techniques produced evidence (Raspopovic et al., 2014) that mouse digits are formed with an RD system that uses feed-forward techniques–a molecule that affects embryo size also feeds forwards to affect λRD, thus keeping λRD reasonably aligned with embryo size and tending to produce the correct number of digits.

Barkai later proposed expansion-repression (Ben-Zvi and Barkai, 2010; Ben-Zvi et al., 2011), which uses a morphogen A and “expander” species E. A is generated at one end of the field at x = 0, and then diffuses and decays everywhere, again with a characteristic length of λRD=DAKD,A. Thus [A] falls off as x>λRD. Because E is generated only when [A] is less than some repression threshold T rep, then E serves as a way to detect that λRD is shorter than L. They propose that E diffuses very quickly and causes λRD to increase everywhere (either by increasing D A or by decreasing K D,A). By using [E] to alter λRD, they robustly set λRD = L/2 and create exactly one repetition of an RD pattern, proportionally scaled to L (Figure 1). This is a negative-feedback system, where [E] essentially measures L/λRD and then feeds back to adjust λRD.

RD and PI were typically seen as competing theories, until experimental evidence mounted for each being used in different circumstances. Green (Green and Sharpe, 2015) suggested that RD and PI can work together in the same system, e.g., by having RD lay down a gradient that PI then uses; Tewary et al. (2017) is an example of this. So is expansion-repression; it lays down a scalable one-peak pattern, thus creating a morphogen gradient varying from low at one end of an organism to high at the other. Since it is scale independent, the coordinate system scales with the length of the organism.

The central problem is how to take an existing set of features, which may or may not be correct, and move in the correct direction in morphospace. Our understanding of the capability of RD to adapt to the field length L has clearly improved over time; from not at all in Turing’s original work (Turing, 1953), to the simple feed-forward hypothesis for mouse digits (Raspopovic et al., 2014), to negative feedback in expansion-repression (Ben-Zvi and Barkai, 2010).

But the negative feedback in expansion-repression only applies to RD patterns with a single peak. How might we get this level of reliability for, say, a five-peak pattern that creates toes? One way could be by adding reliability to RD by wrapping it in a negative-feedback controller. This would require counting the peaks for the current λRD, comparing it to five and adjusting λRD as needed. But how do we count the number of peaks? While this may seem simple, it is not as easy as it seems; here, we accomplish it using waves.

Why is not it easy? Computation outside of the nervous system is typically done with gene-regulatory networks (GRNs). A GRN is a system of promoters inside a cell that senses messenger signals, turns genes on or off to control protein levels, and in turn creates new messengers. Communication between cells then happens as signals leave cells, and travel by diffusion or gap junctions between neighboring cells. We thus have large numbers of small, independent computing entities, working essentially independently and driven largely by local communication–but unfortunately attempting to solve a global-scale problem.

This type of computation is reminiscent of a cellular automaton (CA). A CA is a large collection of identical processing units (called cells, by analogy to biological cells). Each cell in a CA communicates directly only with its neighbors and indirectly via a chain of cells (analogous to diffusion working most quickly at close range). In fact, our counting problem is similar to the classic majority-detection problem for a CA, where you must look at a field of cells that are either 0 or 1, and determine if there are more 0s or more 1s. While the problem may seem simple, it is not at all so (Gács and Levin, 1978; Fates, 2013).

CAs have been used to navigate morphospace; e.g., by Mordvintsev (Mordvintsev et al., 2020) to robustly create images. However, (Mordvintsev et al., 2020) uses an entire neural network in each cell, which is an unrealistic amount of compute power. Furthermore, like many deep neural networks, theirs is so complex that it is difficult to know if it is correct, let alone understand how to modify it to generate slightly different shapes. CAs have also been used to evolve spiking neural nets for the purpose of evolving artificial brains (de Garis et al., 1993). Most CAs use discrete (e.g., Boolean) state variables, though some CAs use analog state [e.g., the floating-point numbers in (Mordvintsev et al., 2020)].

There are several simulators that are capable of simulating networks of cells with GRNs, as well as reactions and diffusion between the cells. Garmen (Kaul et al., 2023) models GRN nodes as multilevel, with state variables either on, off or medium. PhysiBoss (Letort et al., 2019) is a combination of MaBoSS for Boolean models of the GRN and multicellular behaviour using agent-based modelling (PhysiCell). CompuCell 3D (Swat et al., 2012) is a well-known simulator based on the Cellular Potts model. It can simulate a wide variety of partial differential equations (PDEs). Any of the above simulators would probably have been suitable for this work. We use a simulator based on BETSE (Pietak and Levin, 2016; Pietak and Levin, 2017), largely for our convenience based on familiarity with it. Our simulator supports numerical integration of the PDEs for diffusion and for chemical reactions by efficient implicit techniques.

Waves of various types have been well researched in biological tissues (Deneke and Di Talia, 2018; Durston et al., 2018), regulating outcomes ranging across cell death (Cheng and Ferrell, 2018), proliferation (Anderson et al., 2017), and handedness (Anderson et al., 2017). Differentiation waves (Gordon and Gordon, 2019; Gordon and Stone, 2021) are mechanical waves that have been proposed as the driving force behind embryonic differentiation. It was hypothesized that a mechanical differentiation wave reaches a cell, which then decides (via a “cell-splitter” organelle) between one of two resultant states, and which results in one of two refinements of the cell state as well as potentially launching a new differentiation wave. The computational methods by which a differentiation wave might cause a state decision to be reached are unknown. By contrast, our wave is chemical rather than mechanical, and we are proposing a specific GRN to, specifically, count.

Chhabra et al. (2019) describes a fascinating example of waves during gastrulation of human embryonic stem cells in vitro. After seeding with BMP4, they see waves of WNT and NODAL signaling, with cell differentiation happening along with the waves. They show that the cell fates are not decided as a result of reading static concentration gradients left by RD. Rather, the waves are not part of a RD process but arise from some other method; and the cells decide their fate dynamically during the waves. In fact, the WNT and NODAL waves eventually die and concentrations become homogeneous, but the cell fates remain. They have thus used waves as an integral part of differentiating a gastrula into its three layers. By contrast, we will divide a tissue into an arbitrary number of layers, and easily make some layers larger or smaller, using waves that are part of a robust negative-feedback process.

Waves have been observed in many different biological systems using disparate signaling modalities. Bacillus subtulus biofilms produce systemic oscillations in response to nitrogen stress [Figure 2A (Chou et al., 2022)]. Colonies of facultative social Dictyostelium amoebae produce chemoattractive spiral waves of cAMP [Figure 2B (Fujimori et al., 2019; Singer et al., 2019; Ford et al., 2023)]. In cultured mammalian MDCK cells, Erk signaling waves manifest in response to wounding [Figure2C (Aoki et al., 2017; Hino et al., 2020)]. The Xenopus laevis embryo displays numerous waves from the fertilization-triggered calcium wave that initiates development (Busa and Nuccitelli, 1985) to the cell cycle trigger waves that coordinate early cell divisions (Chang and Ferrell, 2013; Anderson et al., 2017), to travelling gene expression waves that help establish the segmental pattern of the embryo (Jen et al., 1999).

Systems-level wave phenomena are likely under-reported in biological datasets as their detection requires long-term live imaging with faithful dynamic reporters. Many commonly used endpoint techniques, like in situ hybridization and immunohistochemistry, are poorly suited to detect dynamic patterns as they require the tissue to be fixed. Others, like RNA sequencing, remove both temporal and spatial information from the assayed tissue. Recent breakthroughs in live imaging of non-neural tissue signaling have revealed a diverse array of mesoscale temporospatial dynamic patterns.

Perhaps the most studied example of waves of gene expression is the segmentation clock, which helps establish the segmental pattern of the vertebrate axis [Figure 2D (Palmeirim et al., 1997; Soroldoni et al., 2014)]. It originated as the clock-and-wavefront model (Cooke and Zeeman, 1976; Tabin and Johnson, 2001) and typically requires a global clock; some versions (Cotterell et al., 2015) avoid this requirement and have an RD-like flavor. It is able, across species, to produce a variable number of vertebrae.

Furthermore, the segmentation clock is remarkably resistant to perturbation. Dividing cells add noise to the system, lose their phase but are re-synchronized to the phase of their neighbors via local interactions (Zhang et al., 2008). When segmentation is transiently chemically interrupted it will reliably recover, even when clock components are heterozygously mutated (Mara et al., 2007). When paraxial mesoderm is surgically removed and cultured ex vivo it will manifest waves and form segments in its new geometry [Figure 2E (Lauschke et al., 2013)], and even when completely dissociated and re-aggregated it will re-establish coordinated oscillations (Tsiairis and Aulehla, 2016; Hubaud et al., 2017).

Despite decades of research, many details of the segmentation clock–including reasons for its robustness–remain unknown. Some segmentation-clock errors do not seem to self-correct (Pourquie, 2022). Its waves are an integral part of the differentiation process. By contrast, our waves perform measurement, are part of an even more robust negative feedback system, and do not require a clock.

Our work, then, combines the two themes of robust morphogenesis and wavefronts. We introduce a GRN that, when replicated in cells and properly stimulated at one end, launches a wave. We will show that the wavefront counts the peaks in an RD pattern as it passes them. This then allows for a closed-loop negative-feedback system to robustly control the number of peaks. In our system, observable waves of morphogen concentration essentially occur as a side effect of computing the properties of a tissue.

Consider a simple RD pattern on a 1D field of cells. The number of replications of the basic RD pattern will depend on the length L of the field and the pattern’s intrinsic length λRD; we typically get roughly L/λRD pattern repetitions in a field of length L. Expansion-repression, as noted above, can alter λRD so that we always get exactly one peak at the source, independent of L. We go a step further: given an integer goal N, we will alter λRD so as to instead obtain exactly N peaks (Figure 1), again independent of L. To do this, we use a wave, based on an identical simple GRN in each cell, that serves to count the number of peaks. We then wrap the GRNs in a top-level control loop that iteratively adjusts λRD and launches a new wave until the wave counts exactly N peaks. Essentially, we have built a closed-loop negative-feedback goal-seeking machine for morphogenesis. It knows its target pattern shape and adjusts parameters iteratively until that goal is achieved.

We add one more capability. Green (Green and Sharpe, 2015) has proposed systems where RD acts downstream of PI, with a morphogen gradient inducing a gradual increase in λRD so that, e.g., digits in a mouse paw are wider at their distal end than their proximal end (Sheth et al., 2012). Meinhardt (2021) has proposed a similar mechanism in Hydra. Both of these serve to build an RD pattern where λRD varies from a small, tight pattern at one end of the field to a larger λRD at the other end. As our wave counts the peaks in an RD pattern, it leaves behind digital breadcrumbs such that each cell knows its exact ordinal position. A small amount of per-cell logic can then examine those signals and increase or decrease λRD in any given cell(s). As a result (Figure 3), we too can make λRD larger or smaller at different locations in the field; but we can do it arbitrarily, rather than only a simple monotonic increase from one end to the other as in (Meinhardt, 2012).

In order to have closed-loop negative feedback to create a RD pattern with the desired number of peaks, we must be able to count the current number of peaks. That is what our GRN does–but unfortunately, counting is more difficult than it may seem. To demonstrate why, we’ll start with a “strawman” solution whose failings will motivate the actual solution.

A human six-year-old, looking at the patterns in Figure 1, can easily tell that Figure 1A has one peak and Figure 1B has three. He would basically scan the picture from left to right, counting each peak as it occurs. But how can an organism, using only a simple GRN replicated in every cell, perform this task? We might start with a set of signals S0, S1, S2, etc., denoting the number of peaks to any cell’s left.

We are assuming that, for any given cell, [A]me is its own concentration of the RD activator A, and [A]left is the concentration of A in the cell on its left. And while we have expressed this logic in terms of Boolean AND gates and IF statements, it can easily be translated into a GRN (Alon, 2019). Figure 4A shows which signals should be expressed in which locations as per this strawman scheme.

However, multiple issues prevent the strawman GRN from actually producing Figure 4A. The first issue is directionality and loops. Our six-year-old human knows to count by sweeping his vision unidirectionally from left to right. Cells have no inherent concept of left and right; and signaling molecules, traveling by diffusion, move in all directions equally well. We would thus be susceptible to improper counting with a scenario like Figure 4B. How can we implement unidirectional counting when our signalling molecules diffuse equally well in all directions? Issues like this make CA algorithms challenging (Gács and Levin, 1978; Fates, 2013).

But the problems with this strawman solution are not over. We’ve talked about a cell measuring [A]left without explaining how that communication could happen. And even if it could happen, our system thus far is not robust to noisy signals. For example, if [A] had some noise that caused wiggling around the cell where [A] = .2, our GRN might miscount the noise as an extra peak (Figure 5). Biology is noisy (McAdams and Arkin, 1999), and real-life systems must be robust.