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

For convenience, in this section we recap all symbols that have been used globally throughout the paper. We do not include symbols that are used only once and defined at their point of use.

• λRD: the “intrinsic wavelength” of a RD pattern. When the pattern repeats itself, each repetition tends to have length λRD (and cannot be longer than that).

• N: the target number of repetitions of an RD pattern. E.g., for the fingers of a human hand, N = 5.

• [A]left: ditto, but in the cell to the left of the “me” cell. This is only used in the purposely-incorrect “strawman” GRN

We show multiple simulations of the system with different field length L and goal RD peaks N. With a constant cell radius and cell-to-cell spacing, the field length is proportional to the number of cells in the simulation.

We picked our parameters to show a sampling of the system’s capabilities. The first three simulations use a field 200 cells long. In the first simulation, our goal is two peaks: i.e., a morphogen profile of LHLH. We start out (Table 1) with λRD = 7.2 × 10−7, which yields four peaks rather than the desired two. The top-level controller then directs seven iterations of counting, noting that there are too many peaks, and increasing λRD by 10% on each iteration, eventually giving the desired two peaks. The second simulation (Table 2) starts with the same initial conditions but has a goal of N = 5 rather than N = 2. This time, the controller goes through four iterations of decreasing λRD, eventually creating the desired extra peak.

The third simulation (Table 3) shows an interesting quirk of some RD patterns. With the same 200-cell field, we started with a different initial λRD. As noted above, a single field length L can often support various numbers of peaks at the same λRD; it is a dynamic system with multiple stable points, and it is often difficult to know which stable point the system will travel to. We thus see iteration #1 setting λRD = 5.4 × 10−7, which could have given us five peaks but instead gave six. Iteration #3 hopped over five peaks directly to four. Eventually, though, the controller keeps adjusting λRD until it successfully reaches the goal. This ability to compensate for the unpredictability of RD patterns highlights a strength of our closed-loop feedback approach.

The fourth simulation (Table 4) uses a small 100-cell-long field. Unsurprisingly, the initial pattern had only two peaks rather than the four in Tables 1, 2. However, multiple iterations of decreasing λRD succeeded in reaching the target of four peaks.

The fifth simulation (Table 5) uses a long 300-cell field. It now starts with five peaks (rather than the four peaks of a 200-cell field). This time, multiple steps of increasing λRD enable it to reach the goal of three peaks.

We have shown that varying λRD allows us to control the number of RD peaks we create. Once we have achieved the desired peak count, we can then vary λRD locally to further control pattern shape. Figure 3A shows a pattern with 100 cells that was targeted to a three-peak pattern. The peaks are originally at cells 20, 60 and 100 (blue graph). A subsequent simulation then slightly modifies the GRN so that any cell expressing S0H increases its λRD by 60% (Figure 3A, orange graph). Figure 3B is quite similar, except in this case we have modified the GRN so that any cell expressing either S0H or S1L increases its λRD by 60%. In each case, the appropriate segment(s) of the RD pattern increase in length at the expense of their immediate neighbors.

This capability is our digital equivalent of what Green (Green and Sharpe, 2015) calls “RD acting downstream of PI.” In his version, PI first lays down a coordinate system and then RD uses it to affect λRD. In ours, our computation wave first lays down one or more coordinate systems and we can then stretch or shrink any subset of them almost arbitrarily.

Robustness is one of the great challenges in morphogenesis (Lander, 2007), and many strategies have evolved to achieve it. Many of them use negative feedback, which is an extremely common motif in nature (Alon, 2019). For example, mice use RD to reliably pattern fingers and toes (Raspopovic et al., 2014). To then prevent large embryos from having extra toes, mice use fibroblast growth factor (FGF), which affects embryo growth, to also increase λRD. If FGF were the only factor affecting embryo growth, this strategy would be completely robust. However, in addition to the inherent unreliability in RD, any variations in embryo size from sources other than FGF are not controlled for, thus occasionally resulting in four- or six-toed mice.

In error-correcting control systems, you get what you measure. FGF concentration is serving as a proxy (albeit an imperfect one) for the field length L, and being used to control for the fact that increasing L would normally increase the number of RD pattern replications. In other words, the proxy for L is being fed forward to control λRD. However, arguably the most appropriate goal is to preserve the number of toes–and the most reliable way to do that is not to measure L at all, but to directly measure the number of toes as a feedback mechanism.

Expansion-repression basically generates A at one end of the field at some rate G A, measures [A] at the other end, and increases (decreases) λRD when the distal [A] is too small(large). It is thus using the distal [A] as a proxy for L. If L doubles, then expansion-repression will double λRD (e.g., by quadrupling D A). The combination of doubling L and quadrupling D A can easily be shown to exactly restore the original [A] profile. Since [A] affects λRD, which then affects [A], this is indeed negative feedback.

If, instead of doubling L, we double G A, then the distal [A] will originally double. It is also easy to show that if we then double both D A and K D,A, then the profile of [A] will be restored. However, in both these cases, if we restore the distal [A] by any other combination of changing D A and K D,A, then [A] at the source will change, as well as the exact profile. The issue is that we are attempting to preserve an exact profile shape but using [A] at a single point to serve as a proxy for that profile.

Our top-level loop, measuring the number of peaks and altering λRD accordingly, is clearly a negative-feedback system. Specifically, our feedback variable is the number of pattern peaks, which is exactly the final variable most important to control. Thus, as long as our feedback system itself is operational, we will be immune from changes in other, non-feedback variables.

This is particularly important since RD patterns are not fully predictable. As shown by our simulation #3 and also noted in (Werner et al., 2015), while a pattern with characteristic length λRD will be stable on a field of length L, it may also be stable on any field of length longer than L. Similarly, increasing λRD such that (as noted above) a six-peak pattern is no longer viable may lead instead to a four-peak pattern, even though five peaks would be feasible. A field can thus often stably sustain a choice of multiple RD peak counts. Each choice has its own stability region around it, where unstable initializations will flow to that particular stable point. The stability regions are often difficult to predict.

The basic building blocks of biology are almost always noisy and difficult to predict (McAdams and Arkin, 1999). Building systems that nonetheless work reliably is thus often difficult, and our case is no exception. While our basic RD patterns can be difficult to use, the most effective feedback system–one where we close the loop by directly monitoring the variable we care most about–is quite effective. This is exactly the system we have built.

While it is clear that morphogen gradients exist in nearly all complex organisms, it is less clear exactly how they are used. The gradients may be quite small (Warmflash et al., 2014; Tewary et al., 2017; Kaul et al., 2023). Furthermore, it is not known whether organisms read the morphogen concentration, its concentration gradient, a fold-change between neighboring cells, or something else (Simsek and Ozbudak, 2022). Since it is impractical to read a morphogen gradient with more than perhaps 10 levels (Simsek and Ozbudak, 2022), a simple gradient is clearly not enough to drive differentiation of a human body, which has over 3 × 1013 cells (Bianconi et al., 2013). A more plausible strategy is to divide and conquer, where the body first divides into (e.g.,) organs, which then form and differentiate independently from each other, each using its own smaller coordinate system. More complex organisms might create their form using even more levels of hierarchy. The task of forming a foot might involve subdivision into five smaller pieces, each with its own coordinate system, and then using a common “toe routine” to further develop each identical piece.

Our closed-loop RD system can do this easily, partitioning the partially-grown field of cells into smaller pieces, each with its own coordinate system that can be used by PI. We may further want different toes to take up different amounts of space, with a big toe typically being wider than the others, which would use our capability for unequal subdivision.

Our system of leaving digital breadcrumbs, via the Pre* signals, is a powerful means of letting each cell know which segment it is part of. It is somewhat reminiscent of Drosophila embryos, where maternal genes cause expression of the gap genes Kni, Hb, Kr and Gt; the embryo cells decode their position by decoding the combination of these four morphogens. It is believed (Petkova et al., 2019) that the embryo decodes these signals near optimally, though it is not exactly clear how. Since our Pre* signals are digital, it is quite easy to decode them; but it takes more morphogens than the Drosophila scheme (which uses analog values and thus requires fewer signals, but is less noise resistant).

Our work is complementary to Chhabra (Chhabra et al., 2019) in several ways. While we both use wavefronts as part of morphogenesis, the usage is somewhat different. Their wavefront is proposed as controlling differentiation; it is part and parcel of the cells’ fate decision. Ours, however, is more of a pure analysis wave, used simply to characterize an existing morphogen pattern as part of overall negative feedback.

Chhabra experimentally shows the occurrence of wavefronts, and that they are linked to cells’ fate decisions. They note that these decisions seem tied to the time between the WNT and NODAL wavefronts. However, they did not know the exact mechanism for this, nor how cellular fate decisions were remembered after the wave had passed. Our work could suggest a specific implementation mechanism for both of these unknowns; we detect and initiate computation upon detecting the WNT wavefront, and could easily measure the time between that and the subsequent NODAL wavefront.

Our work assumes that cells communicate by diffusion, with the diffused signals being potentially regenerated at each cell. However, Chhabra (Chhabra et al., 2019) shows that their signals are not regenerated, and are instead presumably traveling solely by diffusion. Unaided diffusion is an order-N 2 process, while their results show the waves traveling at constant speed. They thus suggest that some form of active transport may be involved. Others (Reddien, 2018; Bischof et al., 2020) have also seen a role for active transport of morphagens. If indeed this is the case, it would greatly simplify the GRN in our work–since the majority of our complexity is to implement effectively-unidirectional signaling, while active transport is often naturally unidirectional.

Another way to reduce the complexity of our GRN might be the imposition of an electric field. This would cause charged signaling ions to travel mostly unidirectionally, though diffusion still occurs–it would be interesting to see if this could simplify our model. Electrical waves have been observed in conjunction with potassium concentration in bacterial colonies (Prindle et al., 2015) and in Xenopus embryos (Vandenberg et al., 2011); our hypothesis may turn out to be a reason for those waves.

Tewary et al. (2017) describes a system where an RD system lays down one repetition of a pattern, which is then interpreted by PI to create the layers of the human gastrula. They show experimentally that increasing the size of the system can result in inadvertent repetitions of the RD pattern. Since larger human embryos do not have multiple gastrulas, they conclude that a higher-level system is probably preventing that. Our work fits nicely with this view; we have proposed a candidate high-level system.

A limitation of our approach is that we have posed our problem as starting with an existing field of cells and subdividing it. However, most embryos are growing at the same time as cells are differentiating, and cells are migrating as well. Our simpler case is not uncommon, though. Mammalian and avian blastoderms can divide in two to create identical twins; each of the two new embryos then reforms itself, differentiating anew to alter each of the existing embryos. Planarian morphallaxis is another interesting case; when an adult planarian is cut into fragments, each fragment can regrow into a full new worm. However, since a fragment may be missing a mouth or indeed an entire digestive system, the fragment cannot increase its mass until it has the capability of eating. It thus undergoes morphallaxis (Reddien, 2018; Reddien, 2021), where the fragment reforms itself into a fully formed but small-scale planarian, and then eats and grows.

We have similarly discussed the problem as occurring in discrete steps; the top-level controller measures a morphological property via a computation wave, then alters parameters, waits for them to take effect, and repeats until it converges to the target. Others (Friston et al., 2015; Friston, 2019) have suggested that the process is more asynchronous and continuous; that each cell continuously monitors if its expectations for the cells in its neighborhood match what it currently senses, and migrates or differentiates so as to bring the two into accord. The large collection of independent agents would march towards convergence by a free-energy-minimization process. We do not know which hypothesis is correct; it could be that a combination of both happens, and computation waves are used to calculate their form of free energy. In the case of RD patterns, however, they tend to have to converge to a fixed pattern before their shape makes any sense at all, which would seem to fit our approach.

The use of discrete signals in negative-feedback controllers in nature is not uncommon. Bacterial chemotaxis uses an integral-based feedback system (Alon, 2019) to decide how often to tumble (i.e., to try a new direction in a search for food). A receptor can be methylated in any of five locations; more methylation implies a larger frequency of tumbles. The controller, like ours, uses a discrete variable with a small number of legal values to serve its purpose.

This work is purely in silico; we have yet no evidence that this particular negative-feedback system exists in nature. However, it seems fairly clear that some sort of negative-feedback system must exist in morphogenesis. The ability of a mammalian embryo to successfully recover from disturbances as varied as being split in two (e.g., for identical twins) and transient mRNA interference to the early embryo (Vandenberg et al., 2012) would be difficult to explain otherwise. Regardless of whether evolution found exactly this scheme, it can now be used in synthetic biology approaches to engineer novel patterning systems (Tiwary et al., 2015; Velazquez et al., 2018; Santorelli et al., 2019; Silva-Dias and Lopez-Castillo, 2022).