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

We can resolve our second and third issues with a simple trick that is very common in human-designed noise-filtering schemes, called a Schmidt trigger (Holdsworth and Woods, 2002). A Schmidt trigger does not try to detect edges directly, but instead uses hysteresis. Instead of trying to detect an edge at 0.2, it uses two thresholds; e.g., .1 and .3. When the signal falls below .1, it is noted as low; when it rises above .3, it is noted as high; and an edge is counted when we rise from low to high. Any spurious noise between .1 and .3 has been made irrelevant.

What governs the choice of, in our case, .1 and .3 as our new thresholds? The further apart the two thresholds are, the more noise immunity we have. On the other hand, if we make (e.g.,) the .3 too high, we risk having some peaks that never reach that threshold and are not counted.

We have now fixed our robustness-to-noise issue, and in fact also no longer need a cell to measure any concentration other than its own. We have incurred the cost of now needing two signals to count each peak: our new signals S0L, S0H, S1L, S1H, etc., are now generated as per Figure 6.

We must still deal with the problem of directionality, loops and double counting (Figure 4B). Our trick is to take advantage of having a prior process break symmetry and give us a known head and tail, thus enabling a global unidirectional sweep to simply count. In this sense, our model is a contribution to the classic problem of leveraging large-scale morphogenetic order from molecular symmetry breaking (Brown and Wolpert, 1990; Vandenberg and Levin, 2009; Vandenberg and Levin, 2010; Naganathan et al., 2016).

We implement unidirectional signaling by tying computation to a wavefront. Cell #0 (at the left) initiates the wave by generating S0L. When cell #1 sees the wavefront, it looks at the local [A] and decides whether to regenerate S0L or to instead generate S0H. Importantly, it then freezes its decision until the next wavefront (if any) happens. It finally sources the appropriate signaling molecule (S0L or S0H) that it wants to travel to cell #2. Of course, these molecules travel to the left towards cell #0 as well–but it is moot, since all cells on the left have already seen the wavefront, made their decision, and will not change it until another wave comes.

With this system, even though our signaling molecules diffuse equally in all directions, the computation wavefront proceeds only from left to right. The key is that any path from cell #0 to a cell C via a loop will be longer than the simple direct path from cell #0 to C, and thus will arrive at C after the direct signal has arrived, and thus cannot affect the action that cell C takes. As the wave hits any cell, that cell examines the local environment (specifically, [A]) and the accumulated state arriving in the wave (i.e., whether the wave is sending S0L, S0H, etc.) to decide what signal to send out in the exiting wave to the next cell to the right.

Here is the final robust GRN in each cell (where the top-level controller now seeds cell #0 with S0L): Pre0L=Pre0L nothing S0L & !very_high& !done Pre0H=Pre0H S0L & very_high & !done S0H& !very_low& !done Pre1L=Pre1L S0H& very_low & !done S1L & !very_high& !done Pre1H=Pre1H S1L & very_high & !done S1H& !very_low & !done … done=Pre0L | Pre0H Pre1H Pre1L Pre2H … S0L=Pre0L S0H=Pre0H S1L=Pre1L S1H=Pre1H

The combination of Schmidt triggers with computation wavefronts counts peaks quite robustly; we will discuss the limits of robustness shortly. Again, while we have expressed this logic in terms of simple Boolean AND, OR and NOT gates, it can easily be translated into a GRN (Alon, 2019), where each signal in the GRN becomes either a protein or a messenger induced by a protein. Finally, assume a top-level controller forcing the leftmost cell to express S0_L.

Cells communicate with each other with the S* signals. Any cell participates in the computation wave by waiting to receive an S* signal, then deciding (based on the local [A]) which S* signal to relay onwards. The Pre* and done signals are local (i.e., do not leave the cell that generates them) and assure that each cell computes only at the wavefront. Any cell, once it receives an S* signal, makes its decision by driving one of the Pre* signals. These signals stay within the cell and are purely for internal calculation. Once a cell drives any of its Pre* signals high, its done (which also remains in the cell) also goes high. This then feeds back to the first set of equations and serves to cut off the Pre* signals from looking at any incoming S* signal any longer. At this point, the self-loop in each Pre* equation takes over, so that whichever Pre* signal is asserted will stay asserted.

Finally, the appropriate S* signal gets driven out of the cell, and stays asserted until new_wave comes in and breaks the Pre* self-loops. New_wave is a global (i.e., widely-diffusing) signal that operates by substantially increasing the degradation rate of the Pre* signals (e.g., by adding a degradation tag). This breaks the self-loop, and thus turns off the Pre* and then done signals in each cell. Morphogenesis is a dynamic process; cells respond to cues throughout morphogenesis, and there will thus be frequent computation waves. Each is preceded by a new_wave signal, which is issued by the top-level controller; while new_wave must reach a high enough level everywhere to act as a global signal, it does not need to reach a level of full spatial homogeneity.

Given that each S* signal is merely a buffer of its corresponding Pre* signal, why bother with the extra signals? The self-loop on each Pre* signal implements our memory of the decision a cell takes; if we tried to put that self-loop on the S* signal instead, then each cell would latch the incoming S* signal before making a decision of its own.

In our GRN, the signals Pre*, done, very_high and very_low are proteins, expressed by genes given the appropriate transcription factors. However, the S* signals must travel between cells and are thus unlikely to be proteins and are thus not the direct output of genes; they may, e.g., be created from the gene products Pre* by chemical reactions.

The next piece in our system is the top-level controller. The controller takes a goal N. It uses the computation wave to count peaks; then compares the count to N and adjusts λRD as needed; and iterates this sequence until we have exactly N peaks. It is conceptually quite simple (Figure 7). The parenthesized numbers below correspond to the flowchart steps in Figure 7.

The system starts with an initial λRD (1) and settles to the resulting RD pattern. This is done with an initial simulation run (2). At that point, the controller launches a new computation wave to count the number of peaks in the RD pattern (3). If, fortuitously, we have exactly the desired number of peaks, then we are done. If not, then the controller implements negative-feedback control. If there are more peaks than our goal, the controller will slightly increase λRD (4). If, on the other hand, there are too few, it decreases λRD (5) and re-seeds (6). In either case, it then re-simulates (2), and the loop iterates until we have attained our goal.

What does the “re-seed the pattern” step mean? It has been known since the 1970s (Bard and Lauder, 1974) that the number of repetitions of a Turing pattern is sensitive to initial conditions. Werner has noted (Werner et al., 2015) that while L/λRD is an upper bound for the number pattern replications in a field of length L, the lower bound can sometimes be 1 (Figure 8). In other words, while we cannot fit (e.g.,) four RD pattern repetitions in a space only large enough for three, it is possible [albeit unlikely (Werner et al., 2015)] for one single pattern repetition to stretch/scale itself up to an almost arbitrarily large field.

The top-level code above thus includes a small trick. When we are increasing λRD, we merely continue the simulation with the larger value. But when we decrease λRD, this may not succeed at creating extra peaks. Instead, it turns out that setting [I] = 0, with [A] rising linearly from 0 at tail to 1 at the head is reasonably reliable at seeding the maximum number of peaks in a given field size. As discussed below, it is not 100% reliable–but our closed-loop controller can successfully work around the unreliable building block.

The top-level controller must be located where it can see the result of the computation wave; for a tail-to-head wave it must live near the head. The controller is small and simple enough that it could easily be implemented as a GRN. There must be only one top-level controller; if it is located in a cell, then there must be a mechanism for it to only be active in one location. For simplicity, we have merely left it as software in our simulations.

Here, again, is the detailed GRN that implements our cellular automaton: Pre0L=Pre0L nothing S0L & !very_high& !done Pre0H=Pre0H S0L & very_high & !done S0H& !very_low& !done Pre1L=Pre1L S0H& very_low & !done S1L & !very_high& !done Pre1H=Pre1H S1L & very_high & !done S1H& !very_low & !done … done=Pre0L | Pre0H Pre1H Pre1L Pre2H … S0L=Pre0L S0H=Pre0H S1L=Pre1L S1H=Pre1H

We implement the logic equations, as is often done, as Hill functions. The Pre* gates are the most complex: e.g., for Pre0H (Eq. 1): Pre0H=kvPre0H.531+Pre0H.53+S0L/.331+S0L/.33∙A/.331+A/.33∙11+done/.53+S0H/.331+S0H/.33∙A/.131+A/.13∙11+done/.53 (1)

Robust operation of the computation wave places some constraints on system parameters. For example, the S* signals are meant to travel by diffusion to their nearest neighbors. A molecule S generated at a constant rate G S moles/s at x 0 and diffusing freely will have its concentration given by Gx=GSKD,Se−λSx−x0 , where K D,S is the decay rate for S, λS=DSKD,S and D S is the diffusion rate for S. Clearly λS must be large enough for the signal to reach its nearest-neighbor cells, which is hopefully easy. However, it must also be small enough that the S* signals do not travel further than one half cycle of the pattern. So, e.g., S1H will first be generated at the cell C 0 where S1L is seen and [A]>.3; it will correctly be generated at cells further to the right until we reach a cell C 1 where [A]<.1. At that point, S2L will correctly be generated, and regenerated at successive cells until we reach a cell C 2 where [A]>.3 again. But what if S1H (traveling by diffusion from cell C 1) reaches cell C 2 before S2L (being regenerated at each cell from C 1 to C 2) does? In that case, cell C 2 will incorrectly see S1H as the incoming wavefront, and will express Pre1H immediately, and will be locked into that decision before it sees S2L and tries to express S2H. As a result, the count will be too low by one.

How do we avoid this issue? Diffusion is an order (distance squared) process and thus slow over long distances. Our trick then is simply to be sure that the half-pattern-length distance is always more than just one or two cells. This imposes a minimum size on λRD.

In most of the RD literature, the reaction and diffusion can occur anywhere in a homogeneous fixed-length field. Simulating cells that interact with RD then requires interfacing this continuous RD field with cells that are at discretized locations. Often (e.g., the sims in (Kaul et al., 2023) that explain the results from (Tewary et al., 2017; Chhabra et al., 2019)), this is done by stepwise alternation between two simulators; first simulate RD, then each cell reads the concentrations at its discrete location to give to its GRN simulation. Since it is uncommon for RD systems to allow an analytic solution, numerical simulation is used for RD, which discretizes the simulation area in any case. The RD activator A and inhibitor I might then be small molecules that could diffuse through a cell membrane and then interact anywhere in the field.

Our simulations are done with BITSEY, a faster but less powerful version of BETSE (Pietak and Levin, 2016; Pietak and Levin, 2017). All of the code to reproduce the simulations is publically available at (repository location to be added after the blind-review process https://gitlab.com/grodstein/bitsey/-/tree/master/RD). BITSEY treats the world as a collection of cells interconnected by gap junctions (GJs). Since the proportion of cell surface area covered by GJs is typically small, reactions in a cell tend to happen on a faster scale than inter-cell communication (Pietak and Levin, 2017): BITSEY thus models molecular concentrations as being constant within a cell rather than subdividing a single cell into multiple finite volumes, leading to fast simulation. While BITSEY models ion channels via the Goldman-Hodgkins-Katz (Keener and Sneyd, 2009) model, our current work does not use them. Communication between cells through GJs is modeled by a simple diffusive flow, discretizing Ficke’s Law.

It is easy to show that, as long as the RD pattern length is substantially longer than the diameter of a cell, BITSEY’s model is mathematically equivalent to a simple discretization of the diffusion equations in 1D, where BITSEY essentially chooses the discretization length of numerical integration to be the cell diameter. This, of course, may not be numerically optimal; if the RD pattern size is far larger than a cell, then BITSEY’s model is unduly compute expensive, and if the RD pattern size is on the order of the cell diameter then BITSEY’s model is numerically inaccurate. However, the model does give us one feature–the exact same simulator core handles both RD simulation and GRN simulation, and both could happen simultaneously if desired.

Each cell holds one GRN. Rather than modeling the GRN in each cell by Boolean or multi-level logic [as, e.g., (Kaul et al., 2023)], BITSEY models the differential equations as per a Hill model (Alon, 2019). Thus, the discretization in space is a single cell; the discretization in time is determined by the speed of change in cell-concentration changes, both from the GRNs and from GJs.

There are countless varieties of RD equations in the literature. Most would work equally well in our system, but we had to choose one. We used a very simple RD pattern taken from (Werner et al., 2015), where (Eq. 2) ∂A∂t=gAGA,I−kD,AA+DA∂2A∂x2,∂I∂t=gIGA,I−kD,II+DI∂2I∂x2 and GA,I=11+I/Ah (2)

There are numerous more-complex RD systems in the literature, including entire categories of new systems that do not place stringent requirements on diffusivity ratios (Marcon and Sharpe, 2012; Landge et al., 2020). The choice of RD system is immaterial to this work, as long as it can be manipulated by λRD and it forms an oscillating pattern of some species that the per-cell GRN can see.

The feedback mechanism in our system operates by repeatedly adjusting λRD. Since λ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), this can be done by either adjusting D A or k D,A. Nature has access to multiple means of controlling both of these (Lander, 2007): degradation tags, competing reactions to bind a morphogen, and even lipid modification to affect diffusivity. We choose to alter D A by changing GJ density. Altering the density of GJs between cells changes the effective diffusion rate of molecules passing through those GJs. While this choice is largely for convenience, there is evidence that such a mechanism does exist (Kleber and Jin, 2021; Tripathi, 2021). Furthermore, there is substantial evidence (Iovine et al., 2005; Watanabe et al., 2006; Watanabe and Kondo, 2012) in fish models of RD molecules traveling through GJs, so this seemed to be a reasonable choice. Once more, we choose to use this mechanism in our simulations for compatibility with an existing simulator, and it is not central to our results.

Our simulation methodology imposes limits on the values of N (the target number of RD pattern repetitions) and L (the field length). As N gets larger, we are asking for each RD pattern repetition to occur across a smaller number of cells. At some point, quantization errors make the RD patterns less stable. This tends only to be a problem for simulations that (e.g., in order to improve simulation time) use a biologically unrealistically-small field length.

The top-level controller, as noted, is implemented as Python code for convenience. It simply iterates between calling the main simulation engine to let the RD pattern stabilize, simulating again to kick off a new computation wave that counts pattern peaks, and adjusting λRD accordingly. The simulation engine treats the RD simulation no differently than the GRN simulation and in fact simulates both concurrently. The simulation run that simulates the RD pattern to steady state is also simulating the GRN–but since the top-level controller has not seeded cell #0 with S0L, the GRNs never see a wave and are inactive.