In the following, we consider a class of developmental models based on the combination of (at least) two different transcriptional modules. Biologically, those two modules correspond to two sequential developmental phases. The main assumptions are that those transcriptional modules are globally regulated for multiple genes at the same time (which could be done for instance through chromatin global regulations) and that there is a continuous transition from one to the other. Here, we focus on metazoan segmentation and regionalization, but the formalism might be applicable to other patterning processes where both an initiation phase and a maintenance phase have been described.

We use ordinary differential equations to model our system. Calling $P$ a vector encoding the state of all proteins in any given cell (typically $P$ corresponds to concentrations of proteins), a generic single-cell equation describing all models presented in the following is:

In Equation 1, variable $g$ encodes an external control parameter of the developmental transition. For example, $g$ could be an external morphogen concentration driving patterning, but more complex situations with feedback are possible, where $g$ could also be part of the system (e.g. the phase difference between oscillators [Beaupeux and François, 2016; Sonnen et al., 2018]). For simplicity, we rescale variables so that $g$ is constrained between 0 and 1. The terms $D\left(P\right)$ and $S\left(P\right)$ correspond to different sets of modules, their influence on the dynamics being weighted by functions ${\theta }_D\left(g\right)$ and ${\theta }_S\left(g\right)$, respectively. The term $\eta \left(g,P\right)$ encodes the noise. Finally, $C\left(P\right)$ represents dynamical terms that are independent of the transcriptional modules, such as protein degradation.

We focus here on the simplest two-module case, where $S\left(P\right)$ encodes a multistable system (i.e. with multiple fixed points at steady state) and $D\left(P\right)$ a dynamic system (i.e. oscillatory). In this situation we will assume ${\theta }_S\left(0\right)=1$, ${\theta }_S\left(1\right)=0$, ${\theta }_D\left(0\right)=0$, and ${\theta }_D\left(1\right)=1$, meaning that for $g=1$ the network is in a pure dynamic phase, while for $g=0$ the network is multistable. Details on the specific forms of $D\left(P\right)$, $S\left(P\right)$, ${\theta }_D\left(g\right)$ and ${\theta }_S\left(g\right)$ are given in the following and in the Appendix. We study two types of models: gene-network like models where $D\left(P\right)$ and $S\left(P\right)$ explicitly model biochemical interactions between genes (such as transcriptional repression), and geometric models where $D\left(P\right)$ and $S\left(P\right)$ directly encode flows in an abstract 2D phase space, similarly to Corson and Siggia, 2017.

We model an embryo as a line of cells, corresponding to the antero-posterior axis. The dynamics within each cell (position $x$) is described by Equation 1. The only difference between cells is that the dynamics of $g$ is a prescribed function of $x$, for example we assume that there is a function $g\left(x,t\right)$ describing the dynamics of a morphogen. We focus on the transition between the two regimes as $g$ continuously changes from 1 to 0 in different cells as a function of time. We will typically consider a sliding morphogen gradient moving along the antero-posterior axis with speed $v$, described by $H\left(s\left(x-vt\right)\right)$ where the function $H$ encodes the shape of the morphogen, and parameter $s$ is a measure of the gradient’s spatial steepness.

We also include noise in the system with the help of an additive Gaussian white noise. For gene networks, we follow an approach similar to the τ-leaping algorithm (Gillespie, 2001), where the variance of the noise corresponds to the sum of the production and the degradation terms (approximating independent Poisson noises). A multiplicative noise intensity term ${\displaystyle \sqrt{1/\mathrm{\Omega }}}$ is introduced, where Ω can be interpreted as the typical concentration of the proteins in the system, so that bigger Ω corresponds to lower noise. In addition, we add diffusion coupling the cells in the stochastic gene network models. For the geometric model, the variance of the noise is held independent of the position $x$. A more detailed description of the noise and diffusion terms is provided in the Appendix.

All source codes and data used for this paper are available at: https://github.com/laurentjutrasdube/Dual-Regime_Geometry_for_Embryonic_Patterning (copy archived at https://github.com/elifesciences-publications/Dual-Regime_Geometry_for_Embryonic_Patterning).

The model described above aims at being generic, but one can relate it to existing developmental networks. A two modules dichotomy was initially proposed in the Drosophila context, where two enhancers (early and late) were observed for Krüppel and knirps (El-Sherif and Levine, 2016). An extension of this model has been proposed in Zhu et al., 2017 for a gap gene cascade under control of the maternal gene Caudal, which would thus play the role of $g$ (a detailed model is reproduced in the Supplement, see Figure 1—figure supplement 2). The dynamic module corresponds to a genetic cascade comprising hunchback, Krüppel, mille-pates, and giant. Each gene in this sequence activates the next one, and later genes repress earlier ones. For the static module, each gene self-activates, and hunchback and Krüppel repress one another. The situation is less clear for vertebrates, given the plethora of oscillating genes and possible redundancy in three different pathways, specifically Notch, Wnt and FGF (Dequéant et al., 2006). The two-module system we consider assumes both the dynamic and static regimes are realized by the same set of genetic components (genes and/or signalling pathways). Notch signaling is an ideal candidate to be a component of both a dynamic and a static regime that might mediate vertebrate somitogenesis, since Notch is implicated in both the core segmentation clock of vertebrates (e.g. via Lfng [Dale et al., 2003], or the hes/her family in zebrafish [Lewis, 2003]) and oscillation stabilization (Jiang et al., 2000; Oginuma et al., 2010). Importantly, several genes of the Notch signalling pathways (e.g. DeltaC in zebrafish) are first expressed in an oscillatory manner then stabilize in striped patterns, as expected in our model (Giudicelli and Lewis, 2004; Wright et al., 2011). There are also multiple genetic interactions between members of the Notch pathway, in particular again Delta genes (Schwendinger-Schreck et al., 2014), with different roles and changes of regulations in the dynamic vs static phase (see e.g. [Wright et al., 2011]).The oscillation itself could be mediated through one or several negative feedback loops in this pathway (Lewis, 2003), and stabilization could be realized through one of the multiple Notch signaling interactions (possibly via cell coupling, similarly to what is observed in other systems [Corson and Siggia, 2017]).

In Zhu et al., 2017, it was suggested that the transition from a ‘wave-like’ behaviour to a static pattern during Tribolium segmentation was mediated by a smooth transition from one set of modules (corresponding to the oscillatory phase) toward another (corresponding to the fixed pattern). This explained the ‘speed-gradient’ mechanism where the typical time-scale of the dynamical system depends strongly on an external gradient (in this case, the concentration of the transcription factor Caudal). In the Appendix, we further study the associated bifurcation, and observe that new fixed points corresponding to the stabilization of gap gene expressions appear on the dynamical trajectory of those gap genes (Figure 1—figure supplement 2F). In simple words, the gap gene expression pattern slowly ‘freezes’ without any clear discontinuity in its behaviour from the dynamic to the static phase, which is reminiscent of the ‘infinite-period’ scenario displayed on Figure 1.

We first aim to generalize this observed property. A simple way to generate waves of gene expressions (as in the gap-gene system described above) is to consider an oscillatory process, so that each wave of the oscillation corresponds to a wave of gap genes. We are not saying here that the gap-gene system is an oscillator, but rather that its dynamics can be encompassed into a bigger oscillator (Verd et al., 2018). The other advantage of considering oscillators is that we can better leverage dynamical systems theory to identify and study the bifurcations. Furthermore, it allows for a better connection with oscillatory segmentation processes in vertebrates and arthropods.

We thus start with an idealized gene regulatory network with three genes under the control of two regulatory modules (Figure 2). In the dynamic phase $D\left(P\right)$, we assume that the three genes are oscillating with a repressilator dynamics (Elowitz and Leibler, 2000), so that the system keeps a reference dynamical attractor and an associated period. In the static phase $S\left(P\right)$, we assume that the module encodes a tristable system via mutual repression (Figure 2A).

Figure 2 with 2 supplements see all

Download asset Open asset

We study the dynamics in a simulated embryo under the control of a regressing front of $g$ (Figure 2B). Transition from the dynamic module to the static module is expected to form a pattern by translating the phase of the oscillator into different fates, implementing a clock and wavefront process similar in spirit to the one in François et al., 2007. We compare two versions of this model, presenting the two different behaviors that we found. In Model 1 (Figure 2C–H), the weights of the two modules are non-linear in $g$: ${\theta }_D\left(g\right)=g^2$ and ${\theta }_S\left(g\right)={\left(1-g\right)}^2$ (Figure 2C). In Model 2 (Figure 2I–N), the weights of the two modules are linear in $g$: ${\theta }_D\left(g\right)=g$ and ${\theta }_S\left(g\right)=1-g$ (Figure 2I). We note that the initial and final attractors of both models are identical. Importantly, only the transition from one set of modules (and thus one type of dynamics) to the other is different. This two-module system thus offers a convenient way to compare the performance of different modes of developmental transition while keeping the same ‘boundary conditions’ (i.e. the same initial and final attractors).

Figure 2E and Figure 2K show the kymographs for both models without noise, with behaviors of individual cells in Figure 2D and Figure 2J. While the final patterns of both models are the same (Figure 2F and Figure 2L), giving rise to a repeated sequence of three different fates, it is visually clear that the pattern formed with Model 2 is more precise and sharper along the entire dynamical trajectory than the one formed with Model 1, which goes through a ‘blurry’ transitory phase (compare mid-range values of $g$ on Figure 2E and Figure 2K).

To better understand this result, we plot the bifurcation diagram of both models as a function of $g$ in Figure 2G and Figure 2M. As $g$ decreases, Model 1 is the standard case of a local Hopf bifurcation (Strogatz, 2015), which happens at $g=0.72$. Three simultaneous saddle-node bifurcations appear for lower values of $g$, corresponding to the appearance of the fixed points defining the three regions of the pattern. The behaviour of Model 2 is very different: the fixed points form on the dynamical trajectory, via three simultaneous Saddle Node on Invariant Cycle (or SNIC) bifurcations (Strogatz, 2015). Both models display waves corresponding to the slowing down of the oscillators, leading to a static regime. In Model 1, the time-scale disappears with a finite value because of the Hopf bifurcation (Figure 2H). For Model 2, it diverges because of the SNIC (Figure 2N), suggesting an explicit mechanism for the infinite-period scenario of Figure 1.

To further quantify the differences of performance between the two models, we introduce noise (encoded with variable Ω, see the Model section and the Appendix) and diffusion (Figure 3A–D). We also define a mutual information metric measuring how precisely the phase of the oscillator is read to form the final pattern (Figure 3E, see the Appendix for details), consistent with the experimental observation in vertebrate segmentation that oscillatory phases and pattern are functionally connected (Oginuma et al., 2010). Intuitively, this metric quantifies in a continuous way the number of fates encoded by the system at steady state. Ideal mutual information for the three mutually exclusive genes of Models 1 and 2 gives ${\displaystyle lo{g}_2\left(3\right)\sim 1.6}$ bits of mutual information, meaning that the pattern deterministically encodes the phase of the cycle into three static fates with equal weights. While addition of noise decreases this mutual information as expected (Figure 3E), Model 2 (black curves) always outperforms Model 1 (red curves). For a reasonable level of noise corresponding to a few thousands of proteins in the system, Model 2 can encode 2^1.3 ∼ 2.5 fates, close to the optimum 3. Furthermore, for a given diffusion constant, Model 1 requires a ten times smaller noise level than Model 2 to encode the same amount of mutual information, which thus indicates much better noise resistance for Model 2.

Figure 3 with 1 supplement see all

Download asset Open asset

Those observations suggest that appearance of stable fixed points through SNIC bifurcations rather than through Hopf bifurcations generates a more robust pattern. The superiority of Model 2 can be rationalized in the following way: when there is a Hopf bifurcation, only one fixed point exists for a range of $g$ values, so that all trajectories are attracted towards it. This corresponds to the ‘blurred’ zone in the kymographs of Figure 2 and Figure 3. In presence of noise, the effect is to partially erase the memory of the phase of the oscillation when only one fixed point is present for the dynamics. Conversely, a SNIC bifurcation directly translates the phase of the oscillation into fixed points, without any erasure of phase memory, ensuring higher information transfer from the dynamic to the static phase, and therefore more precise patterning.

The Hopf bifurcation of Model 1 occurs when the weights of the dynamic and static modules become small compared to the degradation term, which generates an ‘intermediate regime’ with one single fixed point after the oscillations of the dynamic module and before the multistability of the static module. The specific form of the weights is not what determines the bifurcation, but rather the presence or absence of an intermediate regime. We confirm this observation with similar 3-gene models that used Hill functions for the weights ${\theta }_D$ and ${\theta }_S$ (Figure 2—figure supplement 1 and Figure 3—figure supplement 1). Interestingly, we get both Hopf and SNIC bifurcations with the same shape for the two weights; the Hopf is obtained by shifting the weight of the dynamic term toward larger values of the control parameter. This effectively generates the required intermediate regime where both weights are small compared to the degradation term.

Hopf and saddle-node bifurcations are ‘local’ bifurcations, in the sense that changes of the flow in phase space are confined to an arbitrarily small region of phase space as the bifurcation is approached. They do not in principle require complex changes of the flow or fine-tuning of the parameters to happen. As such, they are the most standard cases in many natural phenomena and in most theoretical studies. Conversely, SNIC bifurcations are ‘global’ bifurcations (Strogatz, 2015; Ermentrout, 2008): they are associated to changes of the flow in large regions of phase space (e.g. when a limit cycle disappears with a non-zero amplitude) and usually require some special symmetries or parameter adjustments to occur (e.g. to ensure that a saddle-node collides with a cycle).

It is therefore a surprise that SNIC bifurcations spontaneously appear in the models considered here. To better understand how this is possible and if this is a generic phenomenon, we follow ideas first proposed by Corson and Siggia, 2012 and consider geometric (or gene-free) systems. We aim to see if: 1. SNIC bifurcations are generically observed, and 2. a model undergoing a SNIC bifurcation is in general more robust to perturbations than a model undergoing a Hopf bifurcation, with initial and final attractors being held fixed. We thus build 2D geometric versions of the system (variables $y$ and $z$). The dynamic module $D\left(P\right)$ is defined by a non-harmonic oscillator on the unit circle, while the static module $S\left(P\right)$ is defined by two stable fixed points, at $y=\pm 1$, $z=0$ (see Figure 4A, and the Appendix for the equations). Like previously, we build a linear interpolation between the two systems as a function of $g$ and explore the consequence on the bifurcations (Figure 4B–H). Since the flow in the system is 2D, we can also easily visualize it (Figure 4I and Figure 4—video 2).

Figure 4 with 9 supplements see all

Download asset Open asset

In brief, this geometric approach confirms all the observations made on the gene network model of the previous section, and further clarifies the origin of the SNIC bifurcation. Because of the smooth transition between modules, the entire flow in 2D needs to interpolate from a cycle to a bistable situation. When both modules have close to equal weights (around $g=0.5$), the flow and associated cycle concentrate around two future fixed points. This appears in retrospect as the most natural way to interpolate between the two situations since both types of attractor (stable limit cycle, and multiple stable fixed points) are effectively present at the same time around $g=0.5$. For this reason, the oscillations are also more similar to relaxation oscillations, rapidly jumping between two values corresponding to the future fixed points. When $g$ is further lowered, the weight of the static module dominates and ‘tears apart’ the cycle, forming two fixed points.

This situation is so generic that in fact, to obtain a Hopf bifurcation, we need to mathematically reinforce the fixed point at $y=0$ for intermediate values of $g$. More precisely, three things are required to get a Hopf bifurcation. First, we need to add an extra term, the ‘intermediate module’, characterized by a single fixed point at $y=0$ and $z=0$. Without this module, we consistently get a SNIC bifurcation, for all forms of coupling tested, including linear and non-linear couplings (Figure 4—figure supplement 1A–D). Second, we need to set the weight of the intermediate module to 0 for $g=1$ and $g=0$, so that we get the oscillations of the dynamic module at the beginning of the simulation and the bistability of the static module at the end of the simulation. We achieve this by setting the weight of the intermediate term to $g\left(1-g\right)$, which is of second order in $g$. Third, we need to make the weights of the dynamic and static modules smaller than the weight of the intermediate module for intermediate values of the control parameter, that is for $g$ around 0.5. This is achieved by using cubic weights for the dynamic and static modules (Figure 4—figure supplement 1G and Figure 4—figure supplement 2). Weights of the fourth order in $g$ also lead to a Hopf bifurcation (Figure 4—figure supplement 1H), while linear weights lead to a SNIC bifurcation (Figure 4—figure supplement 1E). Interestingly, quadratic weights lead to simultaneous supercritical Hopf and pitchfork bifurcations (Figure 4—figure supplement 1F). The non-linearity of the coupling helps reinforce the relative weight of the intermediate module for values of $g$ around 0.5, but the exact shape of the non-linearity is not crucial. Furthermore, our mutual information metric confirms that the pattern is more robustly encoded when the system goes through a SNIC bifurcation rather than through a sequence of Hopf and pitchfork bifurcations (Figure 4—figure supplement 5F).

Detailed simulations of the gene-free model with the intermediate module and with cubic weights for the dynamic and static modules are shown on Figure 4—figure supplement 2. As expected for a Hopf bifurcation, the flow first concentrates on the central fixed point at $y=0$, before re-emerging in a bistable pattern for lower $g$ (Figure 4—figure supplement 2I and Figure 4—video 1). There are technically two types of Hopf bifurcations, depending on the stability of the limit cycle. A Hopf bifurcation is supercritical (resp. subcritical) if a stable (resp. unstable) limit cycle becomes a stable (resp. unstable) fixed point. All Hopf bifurcations discussed previously are supercritical. However, by changing slightly the dynamic module of the gene-free model (and including the intermediate module as well as cubic weights for the dynamic and static modules) we can also obtain a subcritical Hopf bifurcation (Figure 4—figure supplement 3). During this bifurcation, an unstable limit cycle forms around the origin. This unstable limit cycle coexists with the stable limit cycle of the dynamic module for some range of $g$ values before they annihilate each other during a so-called ‘saddle-node of limit cycles’ bifurcation. Finally, the two fixed points of the static module form during two simultaneous saddle-node (of fixed points) bifurcations (Figure 4—figure supplement 3I and Figure 4—video 3). Again, our mutual information metric confirms that the pattern is more precise when the system goes through a SNIC bifurcation (Model 2), rather than through supercritical or subcritical Hopf bifurcations (Models 1 and 3, resp.) (Figure 4—figure supplement 5E). Taken together, these results suggest that keeping the static and dynamic attractors fixed, patterning is both more generic and more robustly encoded through a SNIC bifurcation than through a Hopf bifurcation, at least in simple low-dimension models.

A concern with the results of the previous section might be that those mathematical models are in fact fine-tuned and too symmetrical, so that in particular when the transition occurs, both new fixed points appear for the same value of the control parameter. Furthermore, real biological networks have no reason to be perfectly symmetrical (although evolution itself might select for more symmetrical dynamics if needed). We thus relax our hypotheses to study a system where parameters and trajectories are not symmetrical (Figures 5 and 6).

Figure 5 with 2 supplements see all

Download asset Open asset

Figure 6 with 1 supplement see all

Download asset Open asset

Going back first to the gene network model, we induce an asymmetry between the fixed points by changing thresholds of repression in the static phase (Figure 5A). The bifurcation diagrams of Figure 5B–C indicate that the asymmetry of the fixed points indeed breaks the simultaneity of appearance of all fixed points in both scenarios. We nevertheless notice that for those changes of parameters, all bifurcations still happen in a very narrow range of $g$ for the SNIC model.

Asymmetry of the fixed points might therefore destroy the advantage of SNIC vs Hopf by creating a transient zone where one of the fixed points always dominates. We thus perform a comparison between Models 1 and 2 with the same asymmetric static enhancers (Figure 5, see also Figure 5—figure supplements 1 and 2, and the Appendix for details). To compare the two cases, we consider different time-scales of the morphogen gradient. The reasoning is that the slower the decay of $g$, the more time the system spends in a region of parameter space without all three final fixed points, allowing the system to relax and ‘lose’ phase information. Conversely, a faster decay of $g$ means that less time is spent in a region with few fixed points, and therefore the patterns are expected to be more robust.

We first decrease the thresholds of repression of gene A by both genes B and C (Figure 5A). Results of these simulations are shown in Figure 5: Model 2 with a SNIC bifurcation still outperforms Model 1 with Hopf and saddle-node bifurcations. In particular, it is again visually clear on kymographs that Model 2 produces a robust and well-defined pattern at any time point of the simulations, while Model 1 gives rise to a much ‘fuzzier’ pattern before the transition. Model 1 produces a robust static pattern only for a steep gradient (allowing to quickly move through the ‘fuzzy’ phase) and a weak asymmetry in the static module (Figure 5E). It is brittle to any change of the dynamics of $g$ (Figure 5H) or to stronger asymmetry in the static module (Figure 5—figure supplement 1E,H). Conversely, Model 2 is robust to different shapes of the morphogen (Figure 5F,I). Only for a strong asymmetry does the system lose one fixed point (Figure 5—figure supplement 1I), but even in this case transitions through a SNIC bifurcation appear superior to transitions through a Hopf bifurcation.

The fragility of the Hopf bifurcation to asymmetries in the parameters can be understood as follows. In the asymmetric versions of Model 1, one of the fixed points of the static term forms at the Hopf bifurcation, way before the two other fixed points form. It is therefore the only attractor available for a large range of $g$ values. However, in Model 2 the same asymmetry only favors one of the fixed points for a small range of $g$ values, generating a robust pattern. Again, we can use the mutual information metric defined above to quantify the robustness of the pattern and confirm the superiority of Model 2 (Figure 5—figure supplement 2J). We also confirmed these results for the case of random modifications of the repression thresholds of all interactions in the static term (Figure 5—figure supplement 2).

The asymmetry introduced in Figure 5 changes the shapes of the basins of attraction and the positions of the fixed points. The geometric model allows to change those features independently. The most generic way to introduce an asymmetry in the system is to fix the positions of the fixed points of the static regime and change only the positions of the basins of attraction (the reason is that the future fates depend on the position of the separatrix between different regimes [Corson and Siggia, 2012]). To replicate this situation in the 2D gene-free models, we move the unstable fixed point of the static term along the $y$ axis. Results of this procedure are shown on Figure 6 and confirm our results on the gene-network based models: Model 2 bifurcates via a SNIC and is always more robust than Model 1. When we change the positions of the fixed points in the static regime to move them away from the limit cycle (still in an asymmetric way), interestingly both Models 1 and 2 now bifurcate via SNICs (Figure 6—figure supplement 1). Furthermore, we see that for Model 1, the amplitude of the limit cycle decreases before the bifurcation, while for Model 2, the amplitude increases (Figure 6—figure supplement 1E).

We conclude from all those numerical perturbations that even with asymmetric basins of attraction and asymmetric parameters, transitions based on SNIC bifurcations are both more generic and more robust than the ones based on Hopf bifurcations, at least in simple low-dimension models.

Our theoretical study suggest that SNIC based transitions are both more robust and more generic than Hopf/saddle-nodes ones. We thus now examine ways to distinguish between Hopf and SNIC based transitions experimentally. The natural method would be to modify the value of the control parameter of the bifurcations (in our case parameter $g$) and check how the attractors of the system change. A recent experimental example can be found in the auditory hair-bundle context where it has been suggested that Hopf bifurcations can be distinguished from SNIC bifurcations by tuning the external driving force (Salvi et al., 2016). However, the situation in development is different from any sensory system since the actual control parameters are not known and are likely combinations of various signaling systems (such as FGF or Wnt). There could also be many compensatory developmental mechanisms, making it difficult to fully take control of the system, as well known by developmental geneticists.

In the absence of a direct control, we have to rely on indirect measurements. The most obvious choice is to use the antero-posterior position along the embryo (or along the PSM in the case of somite formation) as a proxy for the control parameter. In situs are informative of what locally happens in insects (El-Sherif et al., 2012; Zhu et al., 2017). In vertebrates, it is possible to monitor in real-time segmentation oscillators in embryos (Aulehla et al., 2008; Webb et al., 2016; Delaune et al., 2012) and in tissue cultures (Lauschke et al., 2013; Matsuda et al., 2020; Diaz-Cuadros et al., 2020), so that many properties of the cycles can be inferred.

A first relevant metric (used in Salvi et al., 2016) is the shape of the limit cycle as a function of the control parameter. There is a clear contrast between Hopf and SNIC bifurcations on our simulated kymographs: for Hopf bifurcations, as mentioned above the transition zone is ‘blurred’ because the damped oscillations relax toward the unique fixed point, while for SNIC bifurcations the oscillations localize close to the fixed points, and as a consequence the pattern is gradually reinforced (Figures 2, 3, 4, 5; damped oscillations are especially visible in the geometric model). To better visualize the ‘blurry’ transition zone characteristic of Hopf bifurcations in the gene-free model (resp. in the 3-gene model), we compute the distribution of the values of the geometric variables (resp. of the gene concentrations) for different values of parameter $g$ (Figure 7—figure supplement 2). In models with Hopf bifurcations, the distribution becomes very narrow for intermediate values of $g$, which is a direct consequence of the damped oscillations past the Hopf bifurcations. This suggests that measuring experimentally the distribution of gene concentrations across developmental time at fixed antero-posterior positions can help distinguish SNIC from Hopf bifurcations during actual pattern formation processes. The experimental picture, both from in situs and live imaging, does not support a narrowing of the distribution: it shows clear increases of amplitude that seems more consistent with a gradual reinforcement of the pattern, similar to a SNIC (see e.g. [Delaune et al., 2012]).

Another metric that could help identify the type of bifurcation is the time evolution of the period (resp. of the frequency), which presents a discontinuity for both types of Hopf bifurcations, but continuously goes to infinity (resp. zero) for SNIC bifurcations. The local spatial wavelength of the pattern provides a continuous measurement related to the local period (as first used and derived in Giudicelli et al., 2007). Calling $S$ the somite size (or the wavelength of anterior/posterior markers within somites after clock stopping), $T\left(x\right)$ and $S\left(x\right)$ the respective period and local wavelength of the pattern at position $x$ (with $x=0$ being the tail and $x=1$ the front), and assuming the cells move towards the anterior with constant speed, we have

In the posterior $T\left(x\right)=T\left(0\right),$ so that $S\left(x\right)$ is infinite: all cells are synchronized locally. As cells move towards the anterior, $T\left(x\right)$ decays and $S\left(x\right)$ decreases but stays bigger than $S$. If there is a Hopf bifurcation, one expects that the clock will stop with a final period $T\left(1\right)$, corresponding to a critical wavelength of the wave $S\left(1\right)=S \frac{T\left(1\right)}{T\left(1\right)-T\left(0\right)}$, which is strictly greater than $S$. There is therefore in principle a discontinuity between the wavelength of the propagating wave in the PSM and the wavelength of the static pattern. Conversely, if $T\left(x\right)$ goes continuously to infinity, $S\left(x\right)$ decays continuously to converge towards $S$ (as measured in Giudicelli et al., 2007). This is a simple and intuitive experimental outcome: in this situation the wavelength of the propagating wave in the PSM decreases until it exactly matches the pattern of anterior markers. Therefore, a signature of the finite vs infinite period bifurcation might be visible by comparing the local wavelength of propagating waves in the PSM to the wavelength of anterior markers in the pattern. Practically, monitoring the distance between the peaks of the oscillations does not give a clear difference between Hopf and SNIC, the reason being that such a discrete measurement can smooth out an abrupt change in the period and that past the Hopf bifurcation, damped oscillations actually give rise to a ‘ghost’ frequency which can decrease quite significantly (see Figure 2—figure supplement 2). Therefore, a measurement of the local wavelength might be less conclusive than a measurement of the amplitude. Nevertheless, we notice that a SNIC scenario gives rise to a combination of an increase of the amplitude with a continuous shrinking of the wavelength, meaning that many refining waves (corresponding to broad variations of periods) will simultaneously propagate.

Finally, SNIC bifurcations are also accompanied by characteristic changes of the shape of spatial wave profiles that might be observable in experiments. We further compare a SNIC-based model to a phase model where an infinite-period behaviour is explicitly assumed, namely the model of a collection of coupled oscillators from Morelli et al., 2009. A kymograph of the spatio-temporal profile of the frequency imposed on the oscillators is shown in Figure 7A, and the dynamics of the resulting pattern formation process is shown on the kymograph of Figure 7B, with the final pattern on Figure 7C. The most striking difference is observed on the shape of the spatial wave profile as it moves towards the region where the pattern stabilizes. In the infinite-period scenario of Morelli et al., 2009, the phase profile is by construction symmetric (albeit stretched in the posterior compared to the anterior, see Figure 7E). In the SNIC scenario, we see a clear asymmetry in the peaks of the wave profile (as shown schematically in Figure 7D). Indeed, following the peaks of oscillations spatially from anterior to posterior (left to right), we see that the transition from positive to negative values of $y$ occurs in two steps. During the first step, the system stays close to $y=1$ for some time, with $y$ values decreasing slowly, while in the second step the system goes rapidly towards negative $y$ values and reaches $y=-1$ fast (Figure 7F). This phenomenon, which gives rise to a ‘sawtooth’ pattern of propagating waves, is observed in all our versions of Model 2 (and is notably absent from all our versions of Model 1, see Figure 7—figure supplement 1). Those spatial asymmetries are likely due to the asymmetries of the limit cycles of relaxation oscillators, where a system jumps between two (or more) steady states in an asymmetric way (which can also be observed in systems close to criticality, see Tufcea and François, 2015). Our model thus offers a simple explanation of wave asymmetry, solving the long-standing problem of the asymmetry of AP vs PA transitions, which is possibly crucial for segment polarity as first suggested by Meinhardt, 1982.

Figure 7 with 3 supplements see all

Download asset Open asset

In absence of a direct handle on the control parameter, comparison of SNIC and Hopf bifurcations in simulated systems reveal three possible signatures of SNIC bifurcations in the propagating antero-posterior waves of genetic expression: 1. oscillations gradually reinforcing the pattern for SNIC (instead of damped oscillations for Hopf) 2. relaxation-like oscillations leading to some asymmetry in the wave pattern, and 3. a continuous shrinking of the spatial wavelength during the dynamic phase to match the static pattern for SNIC, with many well-defined propagating waves.

For prediction 1, decades of in situs of oscillating genes and of experimental monitoring of Notch signalling pathways argue against damped oscillations in somitogenesis. Weak waves gradually refine into well-defined stripes in the anterior PSM, in most species (see for example the comparison made in Gomez et al., 2008), suggesting an increase of amplitude in oscillations. Consistent with this, visualising segmentation clock oscillations in live vertebrate embryos suggests an increase of oscillation amplitude as cells get more anterior in the PSM (Delaune et al., 2012; Lauschke et al., 2013). A similar observation was also made for gap genes during short-germ segmentation in Tribolium (Zhu et al., 2017).

For prediction 2, in somitogenesis, there is an asymmetry in the wave pattern before stabilization. The transition from the high to low phase within one somite is shallower than the transition between those two phases from one somite to the other (i.e. posterior of one somite to anterior of the next), as detailed in Shih et al., 2015. A possible read-out is the downstream pattern of Cadherin controlling segment boundary formation, which presents a sawtooth profile (McMillen et al., 2016), consistent with the observed asymmetry in the wave pattern. It has been proposed that the main reason for the clock control was precisely to generate such periodic sawtooth pattern, and a SNIC bifurcation offers a plausible mechanism. In insects, it has been proposed that one of the roles of traveling waves of segmentation genes (e.g. eve genes shift) is indeed to provide segment polarity (Rothschild et al., 2016; Clark, 2017; Clark and Akam, 2016), which is consistent with this proposal, although precise quantification of wave asymmetry has not been done in this context to our knowledge. In situs of gap genes in short germ insects show asymmetry between the anterior and posterior borders of the propagating gene expression waves as well (Zhu et al., 2017).

For prediction 3, the slowdown of oscillations during vertebrate segmentation appears quite variable between species. Mouse segmentation clock seems at first more consistent with a relatively sudden Hopf scenario. There is roughly a 2π phase shift between the posterior and the front, consistent with only a moderate change of the period/wavelength (Lauschke et al., 2013). Other species have several propagating waves within the PSM, seemingly more consistent with a broader period regime. In zebrafish, in vivo real-time imaging suggests a wavelength of propagating waves in the anterior PSM of roughly twice the segment length, meaning that the period of the clock at the front is at least half of the period in the tail bud (Shih et al., 2015). However, this measurement likely overestimates the wavelength and thus underestimates the period at the front since it is based on a peak to peak measurement between the last two waves. Thus, the period of the clock might further increase in the anterior, consistent with a period divergence and with measurements of Giudicelli et al., 2007 on in situs, who proposed that the period indeed diverges (Giudicelli et al., 2007). Importantly, as said above, while the period increases the amplitude increases too Shih et al., 2015, which is more consistent with a SNIC bifurcation. In situs in snake PSM clearly show a decrease of the wavelength that seems to continuously match the pattern in the anterior, possibly more consistent with an infinite-period bifurcation as well (Gomez et al., 2008).

A direct experimental measure of the wavelength of the pattern for the last oscillations can be found with Mesp2 (Takahashi et al., 2000; Saga and Takeda, 2001). The consensus on Mesp2 is that it is expressed in the last few waves of oscillations, and that its pattern continuously regresses to reach exactly the wavelength of the anterior pattern. This is consistent with an infinite-period bifurcation. In Oginuma et al., 2010 it is argued that Mesp helps setting the segment boundary in mouse by precisely reading the Notch wave, explaining both its biological role and why it is a convenient marker of the wavelength. We notice that the model used to explain the experiments in Oginuma et al., 2010 is Julian Lewis’ model from the appendix of Palmeirim et al., 1997, thus presenting an infinite-period divergence.

One cannot completely exclude from all those experimental observations a Hopf scenario where waves would not be very damped (or would be amplified for some reason) past the bifurcation, while the period of damped oscillations would become long until the saddle-node bifurcation is reached. A way to experimentally falsify this scenario is to induce a smoother transition and see if damped oscillations appear. Interestingly, both in vertebrates and in short germ insects, changes of Wnt signalling indeed give rise to smoother transitions from the dynamic to the static regime (suggesting that Wnt could be related to the control parameter $g$). In mouse, beta-catenin gain of function mutants give rise to considerably extended PSM toward the anterior (Aulehla et al., 2008), with up to five waves of oscillating Lfng (compared to only one in WT). Thus, this extended PSM qualitatively looks much more like a zebrafish or a snake PSM, with reinforcement of Lfng in situs from anterior to posterior (suggesting a refinement of the pattern and not damped oscillations, see Figure 4g,h of Aulehla et al., 2008). In the anterior PSM of such mutants, the wavelength of the oscillations decreases to match exactly the Mesp2 wavelength (Aulehla et al., 2008, see e.g. Figure 3), and the very last waves move extremely slowly towards the anterior (Alexander Aulehla, private communication), consistent with a divergence of the time-scale. In Tribolium, axin RNAi similarly changes the wave pattern, where smoother and more spatially extended propagation of gap gene expressions are observed in those mutants, with reduced wavelengths compared to the WT (see Figure 4 of Zhu et al., 2017). Those phenotypes are exactly what is expected from an infinite-period bifurcation, again consistent with the SNIC scenario.

Lastly, damped oscillations would also be absent from Hopf scenarios where the bifurcation that destroys the limit cycle occurs at a value of the control parameter that is almost exactly the value at which the system becomes multistable. This happens for example when the supercritical Hopf bifurcation that kills the limit cycle happens at the same value of parameter $g$ than the pitchfork bifurcation that generates bistability (Figure 4—figure supplement 1F). A similar phenomenon occurs in the model presented on Figure 4—figure supplement 4: the two simultaneous subcritical Hopf bifurcations that generate the bistability and the saddle-node of limit cycles that destroys the stable limit cycle happen in a very narrow range of $g$ values. We cannot fully exclude those cases from the experimental data, but we notice they are intrinsically less generic since multiple bifurcations have to occur at the same time, which happens only in our models because of specific parameter choices, symmetries in the equations or specific choices of non-linearity for the weights of the modules.

The most straightforward prediction of the model proposed here is the presence of several global transcriptional modules between strongly interacting genes, directly controlling the smooth changes of developmental time-scale (in a similar way to the ‘speed-gradient’ model in Zhu et al., 2017). Many developmental genes are regulated by multiple ‘shadow’ enhancers (Cannavò et al., 2016). A smooth transition between different enhancers has even been observed for gap genes in Drosophila (El-Sherif and Levine, 2016). Global regulation of transcriptional modules could be biologically achieved through ‘super enhancers’ regulating many genes at the same time (Hnisz et al., 2017). A non-trivial prediction of our model is that the intrinsic time-scale of the system is a function of the relative balance of transcriptional activities of the modules. The transcriptional control described here naturally allows for infinite-period bifurcations, an implicit mechanism in several models of metazoan segmentation. This is to be contrasted with classical models of negative feedback oscillators such as the Goodwin model, where the time-scale is entirely controlled by degradation and is independent from transcription and translation rates (Forger, 2011), and with models of delayed oscillators, where the time-scale is essentially controlled by transcriptional delays (Lewis, 2003).

Our model is controlled by an external parameter $g$. The most natural hypothesis would be that $g$ corresponds to an actual morphogen gradient. As said above, Wnt is a natural candidate, but feedbacks clearly exist with Caudal in Tribolium (Zhu et al., 2017) and FGF in vertebrates. However, in the spirit of the initial wavefront proposal by Cooke and Zeeman, $g$ could also be in some context a temporal variable, for example an effective timer. Recent works on somitogenesis have suggested that the segmentation front could also be coupled to the slowing down of two oscillators (Lauschke et al., 2013), possibly one corresponding to Notch/FGF and the other to Wnt (Sonnen et al., 2018), so that the oscillation could feedback on itself to define $g$. One could also imagine that $g$ is related to more biophysical variables (density, elasticity) (Hubaud et al., 2017). It is important to point out that in our framework the nature of the bifurcation does not depend on the nature of $g$. While it might be difficult to experimentally disentangle feedbacks between the bifurcations and the control parameter from actual properties of the bifurcations themselves, our predictions on the nature of the bifurcations would not change.

An assumption of this two-module framework is that the same genes interact to control the system in both the dynamic and static regimes, giving rise to a smooth dynamical transition during development. This is consistent with what is observed for gap gene dynamics in short germ insects (Zhu et al., 2017). For vertebrate segmentation, we do not know yet mechanistically how both regimes are controlled, but the Notch signalling pathway is implicated in all steps of somitogenesis and in particular is known to gate information from the oscillatory to the segmented regime (Oginuma et al., 2010). An opposite view would be that the transition from dynamic to static regime is de facto sudden (even if it appears smooth for other reasons). Such scenario could be realized in different ways. For instance, different enhancers could regulate completely different sets of genes in the dynamic vs static regimes. The ‘static’ genes would then interact with the ‘dynamic’ genes only briefly during development, ensuring transmission of positional information between the static and dynamic regimes in a very localized region of time and space. In somitogenesis, as said above specific genes are indeed expressed at the so-called ‘front’ (such as Mesp2) and could act like gating processes transferring the information from the clock to an independent patterning system. In this case, we would be back to a sequential point of view where different regimes of development live in different regions of phase space, and the local bifurcation scenario would then be more plausible.

Evolutionary simulations for the evolution of patterning have favored a scenario based on Hopf and saddle-node bifurcations (François et al., 2007). Those simulations did not include multiple enhancers like here, and all transcriptional regulations had essentially to evolve from scratch, possibly suggesting a ‘Ur-segmentation clock’ based on Hopf and saddle-node bifurcations. This scenario is not excluded by our model: in particular there would be no difference between Hopf and SNIC bifurcations under the control of steep gradients of $g$, which shortcut the region where only one fixed point is present. So Hopf/saddle node bifurcations under control of a steep gradient are likely the easiest solution found by (simulated) evolution. However, adding a more complex, enhancer-based evolutionary grammar might allow for more combinatorial use of dynamical attractors, and associated robustness to external perturbations (such as the shape of the $g$ gradient). A SNIC bifurcation might then plausibly evolve from the modularization (West-Eberhard, 2003) of a system based on Hopf and saddle-node bifurcations (both in vivo and in silico). First, since the SNIC bifurcation is the generic scenario that we observe in our framework, it should be easy to discover by evolution once a multi-enhancer system combining an oscillator and a bistable system has evolved. Second, SNIC-controlled segmentation is plastic in the sense that changes of dynamics of the control parameter would change the transient phenotypes (such as the number of oscillating waves) but would still generate a pattern. Modularizations leading to developmental plasticity has been suggested to be an important engine for evo-devo (West-Eberhard, 2003), since it allows for intra-specific variability without impinging on the most important phenotypes (here, segmentation).

Importantly, such plasticity is indeed observed experimentally both for short germ and vertebrate segmentation. For instance, in Tribolium one can considerably modify Caudal gradient dynamics and still see proper patterning (Zhu et al., 2017). It could thus explain how and why there is so much quantitative variability in segmentation mechanisms, such as short vs intermediate germ band segmentation (as suggested in Zhu et al., 2017). In somitogenesis, there is a lot of interspecies quantitative variability in the numbers of waves and rescaled periods (Gomez et al., 2008) while the qualitative dynamics itself appears very conserved see for example (Krol et al., 2011). In other words, a two-module mechanism makes the dynamics both more robust – a generic bifurcation scenario gives precise phase encoding – and more evolvable – one can vary many features of the system (e.g. basins of attractions, dynamics of the control parameter) and still get proper patterning, again a hallmark of developmental plasticity.

In brief, we have discussed a two-module based model of the smooth transition in development from a dynamical regime to a static one. This model explains time-scale divergence, as well as robustness to changes of morphogen dynamics (Zhu et al., 2017) and to noise. It provides a possible explanation for smooth robust transitions in metazoan segmentation, with a non-trivial (global) bifurcation. Further experimental and theoretical studies are required to assess the importance of smooth transitions for encoding dynamic information into spatial patterns of genetic expressions.

In Zhu et al., 2017, two of us proposed a model of Tribolium segmentation relying on the interplay of two sets of enhancers. In short, two sets of enhancers (static $S(P)$, dynamic $D(P)$) were used, where the role of parameter g is played by morphogen Caudal (cad) (El-Sherif et al., 2012; Figure 1—figure supplement 2). $S(P)$ encodes a multistable system and $D(P)$ a sequential cascade of genetic expression of gap genes ($hb$, $Kr$, $mlpt$, $gt$). This system was found to implement a ‘speed gradient’ model, where the speed of traveling waves of gap genes from posterior to anterior depended on the level of cad concentration (Figure 1—figure supplement 2B–D). This led to robust patterning of the embryo (Figure 1—figure supplement 2E) but the mathematical origin of the speed gradient was not explained.

To better understand the underlying dynamics of the system, we consider the time courses of multiple cells at different positions and thus with different final fates. Figure 1—figure supplement 2F shows the projection of the cells’ dynamics on a 2D plane corresponding to the first two genes expressed in the cascade ($Kr$ and $hb$), as well as a typical flow for different values of $cad$ while keeping the other genes ($mlpt$, $gt$ and $X$) at zero. Importantly, for $cad=0.13$, we see the appearance of a new fixed point (green disk on Figure 1—figure supplement 2F).

We make four observations:

• The flow of the system is canalized. The trajectories of the cells stay very close to one another in phase space.

• As cad is lowered, the new fixed point appears very close to the common trajectory of all cells (Figure 1—figure supplement 2F, top row), and clearly separates the trajectories of cells ending up at different fates (Figure 1—figure supplement 2F, bottom row).

• When cad further decreases, the new fixed point moves in the high hb, low Kr region, corresponding to the eventual fate of Cell 1.

• When the new fixed point appears, the flow of cells past this fixed point is slowed down (Figure 1—figure supplement 2G).

These four observations offer a concise explanation to the ‘speed gradient’ model: as the system gets closer to the bifurcation happening at $cad=0.13$, the system is slowing down because of the future fixed points appearing on the trajectory. Intuitively, this is due to the fact that a fixed point corresponds to a frozen state, and thus to an infinite time-scale (static). When cad varies, the system has to interplay between a non-zero time-scale (dynamics) and such infinite time-scale, and it thus makes sense a priori that in between, the time-scale of the system diverges. This mechanism is close in principle to the critical timing proposed in Tufcea and François, 2015.

For all models except the model for Tribolium segmentation and the infinite-period scenario of Figure 7A–E, the following function is used to describe the spatio-temporal profile of the input g, which is treated either as the concentration of a morphogen in the gene network models, or as an abstract control parameter in the gene-free models:

where parameter s controls the steepness of the gradient and v represents the speed at which the gradient moves along the antero-posterior axis. Parameter $x_{\text{osc}}$ allows to generate a few oscillations inside the first simulated cell before g starts decreasing. Note that the position vector x is normalized in all our simulations, such that positions are constrained from 0 to 1. Appendix 1—table 7 lists the values of the parameters used for the gradients of all models (except the model for Tribolium segmentation): the gradients of the infinite-period scenario and of the Hopf scenario used to generate the results of Figure 1B–F and Figure 1G–K, respectively, the shallow gradient used in the 3-gene models of Figure 2, Figure 2—figure supplement 1 and 2, Figure 3, Figure 3—figure supplement 1, Figure 5, Figure 5–figure supplements 1 and 2, and Figure 7—figure supplement 2, the steep gradient used in the 3-gene models of Figure 5, Figure 5–figure supplements 1 and 2 and , and the gradients used in the gene-free models of Figure 4, Figure 4–figure supplements 2, 3, 4 and 5, Figure 6, Figure 6—figure supplement 1, Appendix 1—table 7, Figure 7—figure supplements 1 and 2 and Figure 7—video 1.

In the model for Tribolium segmentation, the role of input g is played by the maternal gene $cad$. The dynamics of $cad$ is modelled with a Hill function:

where the time dependencies of parameters $x^{*}(t)$ and $n(t)$ encode respectively the regression of the morphogen gradient along the antero-posterior axis, and the gradual increase in the steepness of the morphogen gradient:

In deterministic simulations, the initial phase of the genetic oscillation inside a given cell determines in which part of the pattern this cell will end up. This is not necessarily the case in stochastic simulations. To quantify the robustness to noise of a given model for specific values of parameter $\mathrm{\Omega }$ (and of the diffusion constant $D$ in the case of the 3-gene models) it is required to define a metric that measures the accuracy with which the initial phase of the genetic oscillations inside a cell predicts the region of the pattern in which this cell will end up. The specific metric used in Figure 3, Figure 3—figure supplement 1, Figure 4—figure supplement 5 and Figure 5—figure supplement 2 is the mutual information between the initial phase of the oscillator and the final protein expression state of the simulated cells. The mutual information $I(x,y)$ between two discrete variables x and y is given by the following expression:

where $X$ and $Y$ are the sets of possible values for x and y, respectively. Intuitively, the mutual information between two variables quantifies the amount of information obtained on the value of the first variable by knowing the value of the second variable (and vice versa). If the logarithm is in base 2, the units of the mutual information are bits. To measure how precisely the phase of the oscillator is read to form the final pattern, variable x is set to the phase of the oscillation in protein expression at the beginning of the simulation ${\varphi }_i$, and variable y is set to the protein expression state at the end of the simulation ${\displaystyle \begin{array}{}\text{(31)}& & I({\varphi }_i,P_f)={\displaystyle \sum _{P_f}\sum _{{\varphi }_i}p({\varphi }_i,P_f)\log \left(\frac{p({\varphi }_i,P_f)}{p({\varphi }_i)p(P_f)}\right)}& \text{(32)}& & {\displaystyle =\sum _{P_f}\sum _{{\varphi }_i}p({\varphi }_i,P_f)\log \left(\frac{p({\varphi }_i,P_f)}{p({\varphi }_i)p(P_f)}\right)}& \text{(33)}& & {\displaystyle =\sum _{P_f}\sum _{{\varphi }_i}p(P_f|{\varphi }_i)p({\varphi }_i)\log \left(\frac{p(P_f|{\varphi }_i)}{\sum _{{\varphi }_i}p(P_f|{\varphi }_i)p({\varphi }_i)}\right)}& \end{array}}$

To get the second equality, the fact that ${\displaystyle p(x,y)=p(x|y)p(y)}$ for any two variables x and y was used to get rid of the joint probability ${\displaystyle p(\varphi ,R_i)}$, which is not straightforward to evaluate directly. Similarly, the fact that ${\displaystyle p(y)=\sum _{x\in X}p(x,y)=\sum _{x\in X}p(x|y)p(y)}$ for any two variables x and y was used to get rid of $p(P_f)$, which is less easy to compute than $p({\varphi }_i)$. Indeed, ${\varphi }_i$ is sampled uniformly in the simulations of the 3-gene and gene-free models, since the speed of regression of the input g is constant throughout the simulations. In the 3-gene models, the different phases ${\varphi }_i$ are defined as the different states of protein expression along the oscillation cycle generated by the dynamic module ($g=1$). A uniform sample of ${\varphi }_i$ is obtained by sampling this oscillation cycle at constant time intervals for a total time length of one period. In the gene-free models, the different phases ${\varphi }_i$ are defined as the different sets of (y, z) values along the oscillation cycle generated by the dynamic module ($g=1$). Since the oscillations are on the unit circle (centered at the origin) and have a constant speed along the cycle, sampling uniformly the angles from the positive y axis (starting at 0 and stopping at $2\pi$) generates a uniform sample of ${\varphi }_i$.

All codes are written in the python3 programming language (except for two Mathematica notebooks). Commented jupyter notebooks can be found on Github at the following address: https://github.com/laurentjutrasdube/Dual-Regime_Geometry_for_Embryonic_Patterning. This repository also contains folders with the source data files, as well as the source codes used to generate the data files.

• 3-gene_det.ipynb This notebook performs deterministic simulations of the symmetric 3-gene Models 1 and 2, and deterministic simulations of the 3-gene Models 1 and 2 with Hill functions for the weights of the dynamic and static modules. It also performs a bifurcation analysis of these models using the data found in the XPPAUTO_data folder, which also contains the .ode files used to generate the data with the XPP AUTO software (Ermentrout, 2008). Figure 2, Figure 2—figure supplements 1 and 2, and Figure 7—figure supplement 2 show the results obtained with this notebook.

• 3-gene_stoch.ipynb This notebook performs stochastic simulations of the symmetric 3-gene Models 1 and 2, and stochastic simulations of the 3-gene Models 1 and 2 with Hill functions for the weights of the dynamic and static modules. It also generates plots of the mutual information using the data found in the Mutual_info_data folder, which also contains the python codes used to generate the data. Figure 3 and Figure 3—figure supplement 1 show the results obtained with this notebook.

• 3-gene_asym.ipynb This notebook performs deterministic simulations of the asymmetric 3-gene Models 1 and 2. It also performs a bifurcation analysis of these models and generates plot of the mutual information using the data found in the XPPAUTO_data and Mutual_info_data folders, respectively. Figure 5 and Figure 5—figure supplements 1 and 2 show the results obtained with this notebook.

• Gene-free_det.ipynb This notebook performs deterministic simulations of the symmetic gene-free Models 1, 2, 3 and 4. It also performs a bifurcation analysis of these models and generates flow plots using the data found in the XPPAUTO_data and Mathematica_data folders, respectively. Figure 4, Figure 4—figure supplements 1, 2, 3 and 4, Figure 4—videos 1, 2, 3 and 4, and Figure 7—figure supplement 2 show the results obtained with this notebook.

• Gene-free_stoch.ipynb This notebook performs stochastic simulations of the symmetric gene-free Models 1, 2, 3 and 4. It also generates the mutual information plots using the data found in the Mutual_info_data folder. Figure 4—figure supplement 5 shows the results obtained with this notebook.

• Gene-free_asym.ipynb This notebook performs deterministic simulations of the asymmetic gene-free Models 1 and 2. It also performs a bifurcation analysis of these models using the data found in the XPPAUTO_data folder. Moreover, it generates plots of the flow and of the spatial wave profiles. Figure 6, Figure 6—figure supplement 1, Figure 7 and Figure 7—figure supplement 1 show the results obtained with this notebook.

• Hopf_scenario_Fig1.ipynb This notebook performs deterministic simulations of the gene network model evolved in silico in François et al., 2007. Results are shown on Figure 1. It also performs a bifurcation analysis of this model, shown on Figure 1—figure supplement 1.

• Infinite-period_scenario_Fig1.ipynb This notebook performs deterministic simulations of the infinite-period model of Figure 1, which is a simplified version of the model in the appendix of Palmeirim et al., 1997.

• Infinite-period_scenario_Fig7.ipynb This notebook performs deterministic simulations of the infinite-period model of Figure 7, which is adapted from Morelli et al., 2009.

• Tribolium_model.ipynb This notebook performs deterministic simulations of the model for Tribolium segmentation from Zhu et al., 2017. It also generates flow plots and computes the speed of the cells in phase space. Figure 1—figure supplement 2 shows the results obtained with this notebook.