From genes to patterns: a framework for modeling the emergence of embryonic development from transcriptional regulation

Understanding embryonic patterning, the process by which groups of cells are partitioned into distinct identities defined by gene expression, is a central challenge in developmental biology. This complex phenomenon is driven by precise spatial and temporal regulation of gene expression across many cells, resulting in the emergence of highly organized tissue structures. While similar emergent behavior is well understood in other fields, such as statistical mechanics, the regulation of gene expression in development remains less clear, particularly regarding how molecular-level gene interactions lead to the large-scale patterns observed in embryos. In this study, we present a modeling framework that bridges the gap between molecular gene regulation and tissue-level embryonic patterning. Beginning with basic chemical reaction models of transcription at the single-gene level, we progress to model gene regulatory networks (GRNs) that mediate specific cellular functions. We then introduce phenomenological models of pattern formation, including the French Flag and Temporal Patterning/Speed Regulation models, and integrate them with molecular/GRN realizations. To facilitate understanding and application of our models, we accompany our mathematical framework with computer simulations, providing intuitive and simple code for each model. A key feature of our framework is the explicit articulation of underlying assumptions at each level of the model, from transcriptional regulation to tissue patterning. By making these assumptions clear, we provide a foundation for future experimental and theoretical work to critically examine and challenge them, thereby improving the accuracy and relevance of gene regulatory models in developmental biology. As a case study, we explore how different strategies for integrating enhancer activity affect the robustness and evolvability of GRNs that govern embryonic pattern formation. Our simulations suggest that a two-step regulation strategy, enhancer activation followed by competitive integration at the promoter, ensures more standardized integration of new enhancers into developmental GRNs, highlighting the adaptability of eukaryotic transcription. These findings shed new light on the transcriptional mechanisms underlying embryonic patterning, while the overall modeling framework serves as a foundation for future experimental and theoretical investigations.

1 Introduction

A fundamental challenge in developmental biology is understanding how a group of cells is partitioned into distinct identities, where each group is defined by the expression of one or a few specific genes (Negrete and Oates, 2021; Wolpert, 1969; Green and Sharpe, 2015; Müller and El-Sherif, 2020). This process of cell diversification is further subdivided recursively, creating increasingly specialized subgroups, ultimately forming the complex and highly organized structure of the embryo with its diverse cell types. This phenomenon is known as embryonic patterning (Wolpert, 1969; Green and Sharpe, 2015; Müller and El-Sherif, 2020; Wolpert, 1989; Landge et al., 2020; Meinhardt et al., 1991; Roth, 2011; Briscoe and Small, 2015). Embryonic patterning relies on the precise regulation of gene expression—spatially and temporally—ensuring the right genes are activated in the right cells at the right time (Negrete and Oates, 2021; Wolpert, 1969; Spitz and Furlong, 2012). This regulation controls how cells acquire their identities and organize into complex structures. Embryonic patterning, therefore, is an emergent phenomenon, resulting from the intricate regulation of many interacting genes across numerous cells. This concept of emergence parallels statistical mechanics, where complex behaviors in macroscopic systems, like gases or liquids, arise from interactions among many microscopic particles (Anderson, 1972). Although the underlying laws for these particles, such as Newton’s laws or quantum mechanics, are relatively simple, their collective behavior leads to emergent large-scale phenomena. Similarly, in embryonic development, gene interactions across many cells produce the complex patterns guiding embryonic development (Davidson, 2010; Levine and Davidson, 2005). However, unlike well-understood models in statistical mechanics, our knowledge of gene regulation remains incomplete, particularly regarding how individual genes and their interactions drive these developmental patterns.

In emergent phenomena, not all features of the underlying microscopic elements are critical to the behavior of the larger macroscopic system. Similarly, in embryonic development, not all details of transcription at the single-gene level, which are often intricate (Fuda et al., 2009), are essential contributors to the emergent patterns of gene expression that guide development. Therefore, a key step in understanding how gene regulation mediates embryonic patterning is to develop abstract, simplified models (Alon, 2006; Alon, 2007) that focus on the essential mechanisms at the single-gene level while still capturing the emergent gene expression patterns. However, we currently lack such models.

On one hand, reductionist approaches that delve into the molecular details of transcriptional machinery have uncovered vast amounts of information (Spitz and Furlong, 2012; Cramer, 2019; Furlong and Levine, 2018; Fukaya et al., 2016), but they often fail to clarify how these details contribute to the overarching process of pattern formation. On the other hand, existing models of gene regulation in embryonic pattern formation often rely on simplified systems, such as those found in bacteria (Alon, 2006; Alon, 2007; Kalir et al., 2005; Bintu et al., 2005a; Bintu et al., 2005b; Shea and Ackers, 1985). While these models provide insights, they fall short in capturing the much greater complexity of transcription in eukaryotes, especially in metazoans (Kim et al., 2022; Vincent et al., 2016; Samee et al., 2017; Park et al., 2019). This raises the question of which aspects of eukaryotic transcription are essential for mediating pattern formation. Identifying these key aspects, even at a theoretical level, would enable more targeted experimental work aimed at understanding the transcriptional mechanisms that drive embryonic patterning. To address this challenge, we need to establish a modeling framework that links gene regulation at the molecular level to the higher-level process of pattern formation across embryonic tissues. Such a framework should not simply aim to reproduce experimental data, but rather expose the underlying assumptions about how gene regulation occurs, providing new perspectives on the problem. Moreover, to enhance understanding, it is crucial to supplement these models with computer simulations, offering intuitive and straightforward code that researchers can easily use and modify.

In this paper, we propose such a framework by developing simple toy models of transcription and gene regulatory interactions at the single-gene level. Building on these models, we construct basic pattern formation models that generate realistic yet simplified gene expression patterns. Even if some of these models are based on outdated experiments in bacteria, this effort will clarify the assumptions driving our understanding of gene regulation, highlight areas for future investigation, and point to key differences between bacterial and eukaryotic transcription. By making these assumptions explicit, we aim to shed light on the connection between the molecular processes of transcription and the large-scale phenomenon of embryonic pattern formation, opening new avenues for exploration in developmental biology and developmental genetics. To enhance accessibility and encourage further research, we provide accompanying computer simulations with intuitive and simple code for each model.

We begin our modeling effort by focusing on the transcription process of single genes, using basic chemical reaction models. From there, we model how multiple genes interact to form gene regulatory networks (GRNs) that mediate specific functions within single cells. Next, we introduce phenomenological models of embryonic pattern formation, including various versions of the French Flag model and the Temporal Patterning/Speed Regulation models. These models are then extended with molecular or GRN realizations, representing a comprehensive modeling approach that connects transcriptional-level processes to tissue-level pattern formation. As a case study, we then explore the impact of different strategies for integrating the activities of multiple enhancers regulating a single gene on the robustness and evolvability of embryonic patterning. Our simulations suggest that the two-step gene regulation strategy observed in eukaryotes—first enhancer activation, followed by competitive integration of enhancer activities at the promoter—provides a more standardized approach for incorporating newly evolved enhancers into developmental GRNs. This insight emphasizes the evolutionary adaptability of eukaryotic transcriptional regulation and its role in shaping robust embryonic development.

2 Materials and methods

2.1 Modeling the regulation of a single gene

We begin our modeling of transcription during embryonic patterning by focusing on the kinetics of a single gene regulated by a single activator. Transcription initiation is typically controlled by a non-coding DNA sequence associated with the gene. In prokaryotes and simple eukaryotes, this regulatory region is a short DNA sequence located directly upstream of the coding sequence, near or overlapping the promoter, sometimes called the “operator”. Transcription factors (TFs) that regulate the gene bind to the operator, either aiding in the recruitment of RNA polymerase (in the case of activator TFs) or blocking it (in the case of repressor TFs), thereby controlling transcription (Chen et al., 2021).

In higher eukaryotes, transcription initiation is more complex (Fuda et al., 2009) and often occurs in two stages. Regulatory TFs first bind to a non-coding region associated with the gene, known as an enhancer (Spitz and Furlong, 2012), which can be located at variable distances from the promoter, including far upstream, downstream, or within introns. These enhancers influence the recruitment of general TFs to the promoter, where transcription is initiated. For simplicity, we will start our effort of modeling gene regulation by modeling the binding of a single activator to a single regulatory region. This region may be considered a functional combination of both the promoter and an enhancer, abstracting away the details of their interaction. Later in the paper, we will extend this model to explore the case of interactions between one or more enhancers and the gene promoter, capturing the more complex regulatory dynamics of transcription.

In this simple model of gene activation, an activator TF, denoted as X, binds to the regulatory region of a gene Y (Figure 1A). We represent the unbound regulatory region as $D_0$ and the region bound to X as $D_1$. Transcription factors bind and unbind the regulatory region stochastically, with binding and unbinding rates denoted as $k_1$ and $k_{-1}$, respectively (Figure 1A), described by the equation:

$D_0+X\underset{k_{-1}}{\overset{k_1}{\hspace{0.17em}\rightleftharpoons \hspace{0.17em}}}{D}_1(1.1)$

Figure 1

Figure 1. Modeling the regulation of a single gene. (A) Schematic representation of gene activation by a single activator transcription factor. The transcription factor X binds to the unbound regulatory region (D₀) of gene Y, transitioning it to the bound state (D₁). Binding and unbinding occur at rates k₁ and k_-1, respectively. Transcripts are degraded at rate λ. As an illustration, transcript degradation is shown here but will be omitted in subsequent schematics. (B) Simulation of steady-state transcription levels as a function of activator concentration. The graph depicts the steady-state concentration of gene Y transcripts (Y_ss) increasing monotonically with activator X concentration. Saturation is observed when the concentration of X surpasses the dissociation constant (K_x), resulting in a maximum transcription rate of k_t/λ. This simulation corresponds to Computer Simulation 5 (Supplementary Text S2). (C) Schematic representation of gene repression by a single repressor transcription factor. Repressor X binds to the regulatory region of gene Y, inhibiting transcription. The unbound state D ₀ transitions to the bound state D ₁ upon binding of X, with rates k ₁ and k _-1 indicating binding and unbinding dynamics, respectively. (D) Simulation of steady-state transcription levels as a function of repressor concentration. The graph depicts the steady-state concentration of gene Y transcripts (Y_ss) decreasing as the concentration of repressor X increases. Maximum transcription occurs when X is zero, diminishing progressively with higher X concentrations. This simulation corresponds to Computer Simulation 6 (Supplementary Text S2). (E) Modeling cooperativity in gene activation. Schematic representation of the cooperative binding of two copies (n = 2) of activator X to the regulatory region of gene Y. Both copies must bind simultaneously to activate transcription, representing extreme cooperativity. (F) Effect of cooperative binding on gene activation (n = 2). The simulation demonstrates how cooperative binding of two activator molecules leads to a switch-like increase in gene expression. The response curve shows a sharper transition compared to non-cooperative binding. This simulation corresponds to Computer Simulation 7 with n = 2 (Supplementary Text S2). (G) Effect of cooperative binding on gene activation (n = 5). The simulation illustrates that with higher cooperativity (n = 5), the gene’s response to activator concentration becomes more digital, exhibiting an even sharper switch-like behavior. This simulation corresponds to Computer Simulation 7 with n = 5 (Supplementary Text S2). Some elements of this figure was created by BioRender.com.

When the regulatory region is unbound ($D_0$), transcription is inactive, whereas when X is bound ($D_1$), the gene Y is transcribed at a rate $k_t$, producing transcripts of Y, as described by the equation:

$D_1\stackrel{k_t}{\hspace{0.17em}\to \hspace{0.17em}}{D}_1+Y(1.2)$

Gene transcripts (or gene products more broadly) degrade at a rate $\lambda$, as modeled by the equation:

$Y\hspace{0.17em}\stackrel{\lambda }{\to }\hspace{0.17em}\varnothing (1.3)$

For simplicity, we assume that no transcription occurs in the absence of activator binding (i.e., no basal transcription).

The transcriptional activation of gene Y by TF X, described by Equations 1.1–1.3, can be translated into differential equations using the mass action law, as follows [For a detailed explanation of this conversion process, refer to Supplementary Text S1, Supplementary Figure S1, and Computer Simulations 1 and 2 (Supplementary Text S2)].

$\frac{d{D}_1}{dt}=k_1\left(1-D_1\right)X-k_{-1}{D}_1(2.1)$

$\frac{dX}{dt}=-\frac{d{D}_1}{dt}(2.2)$

$\frac{dY}{dt}=k_t{D}_1-\lambda Y(2.3)$

Here, $D_1$ represents the fraction of DNA bound by X (with $D_1+D_0=1$). These three differential equations can be solved to determine $Y\left(t\right)$, the concentration of gene Y transcripts over time (Supplementary Figure S2A; Computer Simulation S3 in Supplementary Text S2). However, because TF binding to DNA occurs on the millisecond timescale, while transcription proceeds on the order of seconds (Lammers et al., 2020), we can apply the quasi-steady-state assumption. This assumption implies that the steady-state levels of X and $D_1$ are established much more rapidly than that of Y. Consequently, on the timescale at which Y varies, we can reasonably assume that:

$\frac{d{D}_1}{dt}=0(3)$

By substituting Equation 3 into Equations 2.1–2.3, we reduce the system to a single equation, which takes the form of the Michaelis-Menten equation:

$\frac{dY}{dt}=k_t\frac{\frac{X}{K_X}}{1+\frac{X}{K_X}}-\lambda Y(4.1)$

Here, $K_X$ is the dissociation constant of TF X, representing the rate at which TF X detaches from DNA relative to its binding rate, and serves as a measure of the binding strength of X to DNA (lower values of $K_X$ indicate stronger binding):

$K_X=\frac{k_{-1}}{k_1}(4.2)$

As demonstrated by Computer Simulations 3, and 4 (Supplementary Text S2), the Michaelis-Menten equation (Equation 4.1) yields results consistent with the three original equations when $k_1,k_{-1}\gg k_t$ (compare the solution of the full equations describing gene activation, Equations 2.1−2.3, vs. Michaelis-Menten equation, Equation 4.1; Supplementary Figures S2B, C, respectively). To calculate the steady-state concentration of Y as a function of changing X, we set the transcription rate $\frac{dY}{dt}=0$, yielding the steady-state solution $Y_{ss}$:

$Y_{ss}=\frac{k_t}{\lambda }\frac{\frac{X}{K_X}}{1+\frac{X}{K_X}}(4.3)$

This results in a monotonic increase in steady-state transcription as the concentration of TF X increases, with saturation occurring at a maximum transcription rate of $\frac{k_t}{\lambda }$ as X exceeds the dissociation constant $K_X$ (see Computer Simulation 5 (Supplementary Text S2), and simulation results in Figure 1B).

The proposed model of gene activation significantly oversimplifies the transcription process, particularly in eukaryotes, by neglecting critical aspects such as chromatin remodeling, nucleosome clearance, transcription pausing and release, transcriptional bursting, and other important details (Fuda et al., 2009; Lammers et al., 2020; Yokoshi and Fukaya, 2019; Adelman and Lis, 2012; Gaertner and Zeitlinger, 2014; Green, 2005; Nechaev and Adelman, 2011). Additionally, the quasi-steady-state (or equilibrium) assumption, along with the simplification of TF binding and unbinding as fast dynamics, overlooks how these factors might influence the slower overall dynamics of transcription (Fuda et al., 2009). While this model’s simplicity has made it useful for gaining insights into gene regulation during development, assuming that slow dynamic processes have no effect on gene regulation may be misguided. As we will discuss in relation to combinatorial gene regulation, slow dynamics can indeed impact the overall transcriptional logic (Scholes et al., 2017).

2.1.1 Thermodynamic state ensemble (TSE) modelling

As shown above, deriving the Michaelis-Menten relationship from reaction equations can be quite tedious, particularly for more complex scenarios involving multiple TFs binding to the cis-regulatory region. A more streamlined approach for modeling transcriptional regulation that produces comparable results is the “Thermodynamic State Ensemble” (TSE) modeling (Bintu et al., 2005a; Bintu et al., 2005b; Sherman and Cohen, 2012). In this method, changes in transcript levels are influenced by the transcription rate ($T_y$) and transcript degradation:

$\frac{dY}{dt}=T_y-\lambda Y(5.1)$

The transcription rate, $T_y$, is expressed as the proportion of TF, DNA, and RNA polymerase configurations conducive to transcription relative to all possible configurations. Implicit in this approach the assumption that all processes involved in transcription are at equilibrium, and so it shares the quasi-steady state assumption as our modeling of gene activation presented above. Below, we outline how to use the TSE strategy to model gene regulation; for a more in-depth and rigorous explanation, readers are referred to the detailed exposition in (Sherman and Cohen, 2012).

As an example, we will use the TSE strategy to re-derive the rate equation governing the binding of a single activator X to the cis-regulatory region of gene Y. In this case, the denominator of $T_y$ accounts for two possible states: unbound DNA (represented by 1) and DNA bound by X (represented by $\frac{X}{K_X}$, where $K_X$ is the dissociation constant of X). The numerator includes only the states conducive to transcription—DNA bound by X, represented by $\frac{X}{K_X}$. The overall expression is then multiplied by the transcription rate $k_t$, which occurs once RNA polymerase binds to the DNA. This yields:

$T_y=k_t\frac{\frac{X}{K_X}}{1+\frac{X}{K_X}}(5.2)$

Substituting this into the differential equation for transcript dynamics gives the familiar Michaelis-Menten expression for a single activator (Equation 4.1):

$\frac{dY}{dt}=k_t\frac{\frac{X}{K_X}}{1+\frac{X}{K_X}}-\lambda Y(5.3)$

Next, we apply the TSE strategy to model gene regulation by a single repressor X (Figure 1C). In this case, the denominator of $T_y$ remains similar to the activation case, but the numerator reflects only the state conducive to transcription—unbound DNA, represented by 1 (i.e., transcription occurs in the absence of repressor X). This yields the following equation for the transcription rate:

$\frac{dY}{dt}=k_t\frac{1}{1+\frac{X}{K_X}}-\lambda Y(6)$

When plotting the steady-state value of Y as a function of X, the transcription rate reaches a maximum at zero concentration of the repressor X, then progressively decreases as X concentration increases (see Computer Simulation 6 (Supplementary Text S2), which results are depicted in Figure 1D).

2.1.2 Modeling cooperativity

Often, multiple copies of the same TF are required to bind cooperatively in order to effectively recruit RNA polymerases. In this section, we will use the TSE approach to model the cooperative binding of n copies of activator TF X for the activation of gene Y. We will assume a case of extreme cooperativity, meaning that all n copies of X must bind simultaneously to the regulatory region of Y, whereas fewer than n copies are unable to bind.

In this scenario, the denominator of $T_y$ represents all possible binding events, which include unbound DNA and the binding of n copies of X to the cis-regulatory elements. The numerator, however, includes only the state conducive to gene activation, which is the simultaneous binding of n copies of X to the regulatory region of Y (see Figure 1E). This yields the rate equation:

$T_y=k_t\frac{{\left(\frac{X}{K_X}\right)}^n}{1+{\left(\frac{X}{K_X}\right)}^n}(7.1)$

This type of equation is commonly referred to as the Hill equation.

Cooperativity introduces a nonlinear, switch-like behavior in the gene’s response to activator binding. The higher the value of n, the more pronounced this switch-like behavior becomes, with a sharper transition from low to high input. As n increases further, the gene’s response approaches digital behavior, as demonstrated by comparing the cases of $n=2$ and $n=5$ in Figures 1F, G, respectively (see Computer Simulation 7 (Supplementary Text S2)). The point at which the gene switches from low to high output is determined by the dissociation constant $K_X$, which characterizes the binding affinity of TF X (as shown in Equation 7.1).

Similarly, the Hill equation for the cooperative binding of n repressors follows a comparable form:

$T_y=k_t\frac{1}{1+{\left(\frac{X}{K_X}\right)}^n}(7.2)$

In both cases, the cooperative binding of TFs leads to a sharp, switch-like regulation of gene expression, reflecting the impact of cooperativity on transcriptional control.

2.1.3 Modeling multiple TFs binding cis-regulatory module

So far, we have considered the binding of one or multiple copies of the same transcription factor (either activator or repressor) to the regulatory region of a gene. In this section, we will extend this by considering the binding of multiple different TFs. Specifically, we will examine the binding of 2 TFs, A and B, to the regulatory region of gene Y. This can be extended to cases involving a greater number of different TFs with straightforward modifications.

The effect of binding multiple TFs depends on the logic that governs their regulation of the gene, often referred to as the “transcription logic” of the regulatory module. Below, we will explore several common transcriptional logic models, though more elaborate scenarios are certainly possible.

2.1.3.1 Two ORed activators

In this model, the binding of either A, B, or both activates the gene (Figure 2A). Using the TSE approach, and assuming there is no cooperative binding between A and B (i.e., the binding of A does not influence the binding of B and vice versa), we obtain the following equation:

$\frac{dY}{dt}=\frac{\frac{A}{K_A}+\frac{B}{K_B}+2\frac{A}{K_A}\frac{B}{K_B}}{1+\frac{A}{K_A}+\frac{B}{K_B}+\frac{A}{K_A}\frac{B}{K_B}}-\lambda Y=\frac{\frac{A}{K_A}}{1+\frac{A}{K_A}}+\frac{\frac{B}{K_B}}{1+\frac{B}{K_B}}-\lambda Y(8.1)$

Figure 2

Figure 2. Modeling multiple transcription factors binding to a cis-regulatory module. (A) Schematic of gene activation by two ORed activators. Transcription factors A and B independently bind to the regulatory region of gene Y. The presence of either A or B is sufficient to activate transcription, illustrating an OR logic gate. (B) Simulation of transcription levels with two ORed activators. The graph shows how the steady-state concentration of gene Y transcripts varies with different concentrations of A and B, reflecting the OR logic. This simulation corresponds to Computer Simulation 8 (Supplementary Text S2). (C) Schematic of gene activation by two competitive ORed activators. Activators A and B compete for the same binding site on gene Y’s regulatory region. Binding of either A or B activates transcription, but they cannot bind simultaneously due to competition. (D) Simulation of transcription levels with competitive ORed activators. The graph illustrates how the competition between A and B affects gene Y’s expression, with transcription levels depending on the relative concentrations of (A, B). This simulation corresponds to Computer Simulation 9 (Supplementary Text S2). (E) Schematic of gene activation by two ANDed activators. Both transcription factors A and B must simultaneously bind to gene Y’s regulatory region to initiate transcription, representing an AND logic gate. (F) Simulation of transcription levels with ANDed activators. This simulation corresponds to Computer Simulation 10 (Supplementary Text S2). Some elements of this figure was created by BioRender.com.

In cases where multiple copies of A and multiple copies of B, each bind cooperatively (with cooperativity values $n_A$ and $n_B$, respectively), but without cooperative interactions between A and B, we derive the following Hill equation

$\frac{dY}{dt}=\frac{{\left(\frac{A}{K_A}\right)}^{n_A}}{1+{\left(\frac{A}{K_A}\right)}^{n_A}}+\frac{{\left(\frac{B}{K_B}\right)}^{n_B}}{1+{\left(\frac{B}{K_B}\right)}^{n_B}}-\lambda Y(8.2)$

As shown in Figure 2B, simulation result of this relationship matches the OR logic as expected [see Computer Simulation 8 (Supplementary Text S2)]. Moving forward, we will focus on the binding of single copies of A and B, as the extension to cooperative binding is straightforward.

2.1.3.2 Two competitive ORed activators

Here, we consider the case where both A and B are activators that bind the same site on the regulatory region of gene Y, leading to competitive binding (Figure 2C). In this case, the denominator includes unbound DNA, DNA-bound A, and DNA-bound B. However, the simultaneous binding of both A and B is excluded because they compete for the same binding site and cannot bind at the same time. The numerator includes the states where A or B binds independently, resulting in the reaction equation:

$\frac{dY}{dt}=\frac{\frac{A}{K_A}+\frac{B}{K_B}}{1+\frac{A}{K_A}+\frac{B}{K_B}}-\lambda Y(9)$

As expected, simulation result of the Hill-function extension of this relationship matches the competitive OR logic as described above [Figure 2D; see Computer Simulation 9 (Supplementary Text S2)].

2.1.3.3 Two ANDed activators

In this scenario, both A and B must bind simultaneously to activate gene expression (Figure 2E). This logic results in the following equation:

$\frac{dY}{dt}=\frac{\frac{A}{K_A}·\frac{B}{K_B}}{1+\frac{A}{K_A}+\frac{B}{K_B}+\frac{A}{K_A}\frac{B}{K_B}}-\lambda Y=\frac{\frac{A}{K_A}}{1+\frac{A}{K_A}}\times \frac{\frac{B}{K_B}}{1+\frac{B}{K_B}}-\lambda Y(10)$

As expected, simulation result of the Hill-function extension of this relationship matches the competitive AND logic [Figure 2F; see Computer Simulation 10 (Supplementary Text S2)].

2.1.3.4 Two ANDed repressors

So far, we have considered activators, but the extension to models involving both activators and repressors is straightforward. As an example, we will examine the case of two ANDed repressors. In this case, only unbound DNA allows for transcription, and the transcription rate $T_y$ is the product of the individual repression effects of A and B:

$\frac{dY}{dt}=\frac{1}{1+\frac{A}{K_A}+\frac{B}{K_B}+\frac{A}{K_A}\frac{B}{K_B}}-\lambda Y=\frac{1}{1+\frac{A}{K_A}}\times \frac{1}{1+\frac{B}{K_B}}-\lambda Y(11)$

The transcriptional logics presented above represent relatively simple cases. However, more complex scenarios could be modeled within the TSE framework by incorporating cooperative binding between different TFs (Kim et al., 2022). However, the modeling approach presented here shares the same limitations as our earlier approach using rate equations (Equations 1.1−1.3 and 2.1−2.3), as it neglects many details of the transcription cycle in eukaryotes. Additionally, TSE models assume that binding and unbinding events are at equilibrium, similar to the quasi-steady state assumption used in the rate equations (Equations 3 and 4.1−4.3). However, the rate equation framework (Equations 1.1−1.3 and 2.1−2.3) is less constrained by assumptions and can, in principle, accommodate more transcriptional details. Employing rate equations to capture these complexities has proven effective in providing insights into gene regulation in eukaryotes, revealing a range of potential mechanisms for achieving combinatorial transcriptional logic beyond those typically assumed within TSE models (Scholes et al., 2017). Nevertheless, for simplicity, we will continue with the TSE modeling results in this study.

2.2 Modeling gene regulatory networks

So far, we have considered the regulation of single genes by one or more TFs. During development, terminal and housekeeping genes may be passively regulated by other transcription factors. However, many TFs regulate other TFs, forming gene regulatory networks (GRNs) (Figure 3A) (Levine and Davidson, 2005). In this section, we will explore examples of GRNs that act within a single cell and how to model them. Later, we will extend this to GRNs that function across multiple cells within embryonic tissues.

Figure 3

Figure 3. Modeling gene regulatory networks (GRNs). (A) Schematic representation of a gene regulatory network (GRN). Multiple genes usually interact to regulate each other’s expression, forming complex networks that mediate specific cellular functions. Shown is a positive feedback loop between two genes as an example. (B) Positive feedback loop between two activator genes. Genes X and Y activate each other’s expression. This mutual activation can lead to bistable states where both genes are either on or off. (B′, B″) Simulation of gene expression states in the positive feedback loop. The graph shows stable states achieved in the system, with both genes being not expressed or highly expressed, depending on initial conditions. This simulation corresponds to Computer Simulation 11 (Supplementary Text S2). (C) Positive feedback loop between two repressor genes. Genes X and Y repress each other’s expression. This mutual repression results in a system where one gene is active while the other is repressed. (C′, C′′) Simulation of gene expression states in the mutual repression loop. The graph demonstrates the system stabilizing with one gene active and the other inactive, based on initial expression levels. This simulation corresponds to Computer Simulation 12 (Supplementary Text S2). (D) Negative feedback loop between an activator and a repressor gene. Gene X activates gene Y, while gene Y represses gene X. This configuration promotes homeostasis, maintaining stable gene expression levels. (D′) Simulation of gene expression in the negative feedback loop. The graph shows how the system maintains stable expression levels of genes X and Y over time. This simulation corresponds to Computer Simulation 13 (Supplementary Text S2). (E) Oscillatory GRN with two activators and one repressor. Gene X activates gene Y, gene Y activates gene Z, and gene Z represses gene X. This network can produce oscillations in gene expression over time. (E′) Simulation of oscillatory behavior in the gene network. The graph displays cyclic fluctuations in the expression levels of genes X, Y, and Z, characteristic of an oscillatory system. This simulation corresponds to Computer Simulation 14 (Supplementary Text S2). (F) Repressilator network composed of three repressor genes. Genes X, Y, and Z repress one another in a cyclic manner, forming an oscillatory system known as a repressilator. (F′) Simulation of oscillations in the repressilator network. The graph shows the periodic expression patterns of genes X, Y, and Z over time. This simulation corresponds to Computer Simulation 15 (Supplementary Text S2). (G) Multi-stable GRN with mutually repressing genes. Multiple genes repress each other strongly, leading to stable states where only one gene is expressed while the others are repressed. (G′) Simulation of gene expression states in the multi-stable network. The graph shows that depending on initial conditions, one gene remains active while others are inactive. This simulation corresponds to Computer Simulation 16 (Supplementary Text S2). (H) Genetic cascade illustrating sequential gene activation. Each gene in the cascade represses the previous gene and weakly represses the next one. Activation of the first gene triggers a domino effect of gene expression. (H′) Simulation of gene expression in the genetic cascade. The graph depicts the sequential activation and repression of genes over time, following the cascade logic. This simulation corresponds to Computer Simulation 17 (Supplementary Text S2). (I) Oscillatory gene network in which a series of genes repress each other in sequence without the last gene repressing the first, resulting in continuous oscillations of gene expression. (I′) Simulation of oscillations in the gene network. The graph shows the cyclic expression patterns of the genes over time. This simulation corresponds to Computer Simulation 18 (Supplementary Text S2). Some elements of this figure was created by BioRender.com

2.2.1 Positive feedback loop

We begin by considering feedback loops involving two genes. If activation is positive and repression is negative, a positive feedback loop results when the product of all regulatory links’ signs is positive (Alon, 2006; Alon, 2007; Longabaugh et al., 2005). Positive feedback loops can function as toggle switches or memory devices—if one of the genes is activated, it remains active.

A positive feedback loop can be realized as either two genes activating each other (Equations 12.1, 12.2; Figure 3B) or two genes repressing each other (Equations 13.1, 13.2; Figure 3C). In a loop of two activators, both genes can be either high (“on”) or low (“off”) (Figures 3B’, 3B”; see Computer Simulation 11 (Supplementary Text S2)). In a loop of two repressors, the GRN stabilizes in a state where one gene is active and the other is repressed (Figures 3B’, 3C”; see Computer Simulation 12 (Supplementary Text S2)). In these equations, $K_{XY}$ denotes the dissociation constant of TF Y binding to the regulatory region of gene X.

$\frac{dX}{dt}=\frac{{\left(\frac{Y}{K_{XY}}\right)}^n}{1+{\left(\frac{Y}{K_{XY}}\right)}^n}-\lambda X(12.1)$

$\frac{dY}{dt}=\frac{{\left(\frac{X}{K_{YX}}\right)}^n}{1+{\left(\frac{X}{K_{YX}}\right)}^n}-\lambda Y(12.2)$

For a repression-based positive feedback loop, the equations become:

$\frac{dX}{dt}=\frac{1}{1+{\left(\frac{Y}{K_{XY}}\right)}^n}-\lambda X(13.1)$

$\frac{dY}{dt}=\frac{1}{1+{\left(\frac{X}{K_{YX}}\right)}^n}-\lambda Y(13.2)$

2.2.2 Negative feedback loop

A two-gene negative feedback loop (Figure 3D) can help maintain homeostasis, ensuring that gene expression remains stable around set values. For instance, the GRN modeled by Equations 14.1, 14.2 maintains gene Y’s transcriptional activity around the activation threshold for TF X binding to its regulatory region ($K_{YX}$), and gene X’s activity around the deactivation threshold for TF Y binding to its regulatory region ($K_{XY}$) (Figure 3D’; see Computer Simulation 13 (Supplementary Text S2)).

$\frac{dX}{dt}=\frac{1}{1+{\left(\frac{Y}{K_{XY}}\right)}^n}-\lambda X(14.1)$

$\frac{dY}{dt}=\frac{{\left(\frac{X}{K_{YX}}\right)}^n}{1+{\left(\frac{X}{K_{YX}}\right)}^n}-\lambda Y(14.2)$

2.2.3 Oscillatory GRNs

While two-gene negative feedback loops typically maintain homeostasis, more complex three-gene negative feedback loops can result in oscillations, depending on network parameters. For example, a GRN with two activators and one repressor (Figure 3E, E’; see Computer Simulation 14 (Supplementary Text S2)) can behave as an oscillator (Equations 15.1–15.3). A GRN composed of three repressors, also known as a repressilator (Elowitz and Leibler, 2000), can also function as an oscillator [Figure 3F, F’; see Computer Simulation 15 (Supplementary Text S1)].

$\frac{dX}{dt}=\frac{1}{1+{\left(\frac{Z}{K_{XZ}}\right)}^n}-\lambda X(15.1)$

$\frac{dY}{dt}=\frac{{\left(\frac{X}{K_{YX}}\right)}^n}{1+{\left(\frac{X}{K_{YX}}\right)}^n}-\lambda Y(15.2)$

$\frac{dZ}{dt}=\frac{{\left(\frac{Y}{K_{ZY}}\right)}^n}{1+{\left(\frac{Y}{K_{ZY}}\right)}^n}-\lambda Z(15.3)$

2.2.4 Developmentally inspired GRNs

So far, we have presented example GRNs with specific wiring capable of executing various functions within the cell. The structures we’ve discussed are not the only possible configurations, as many different GRN architectures with varying parameter values can accomplish the same function. However, GRNs involved in development, particularly those involved in embryonic pattern formation, have distinct characteristics. One notable feature is the minimal reliance on activators; most TFs involved in embryonic patterning are repressors, where activation is achieved through the repression of a repressor, a process called de-repression (Averbukh et al., 2018; Tufcea and François, 2015). For example, most of the patterning genes involved in the patterning of the Anterior-Posterior (AP) and Dorsoventral (DV) axes in the early Drosophila embryo are transcriptional repressors (Clark et al., 2019; Diaz-Cuadros et al., 2021; Moussian and Roth, 2005).

We will now explore some common GRN architectures that perform specific tasks during embryonic pattern formation, relying on repression and de-repression to mediate regulatory relationships between genes. Another important feature of these GRNs is the extensive cross-regulation between the genes, where oftentimes all constituent genes repress each other. Despite this complexity, variations in the strength of repression among the genes can produce meaningful and distinct functions, as illustrated by the GRNs discussed below.

2.2.4.1 Multi-stable GRNs

A multi-stable GRN (Tufcea and François, 2015) functions similarly to the positive feedback loop or genetic toggle switch discussed earlier. However, developmental GRNs typically consist of more than two genes. The simplest extension of the toggle switch to more than two genes is a multi-stable GRN, composed of M mutually repressing genes. This repression is characterized by small dissociation constants (${\theta }_s$), indicating strong repression (Equation 16; Figure 3G). In this configuration, depending on the initial conditions, one gene remains active while the others are inactive [Figure 3G’; see Computer Simulation 16 (Supplementary Text S2)]. Other variations of this realization are possible, where subsets of genes are co-expressed due to weakened repressive links between them (Zhu et al., 2017).

$\frac{d{G}_m}{dt}={\displaystyle \prod _{\forall i}^{i\ne m}}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}-\lambda G_m(16)$

2.2.4.2 Genetic cascade

Another developmentally inspired GRN is the genetic cascade, where a sequence of genes is activated one after another. A simple realization of this cascade involves one gene activating the next, which in turn represses the previous gene, and so on. In a GRN based on de-repression logic (Averbukh et al., 2018; Zhu et al., 2017; Jaeger et al., 2004), all genes are naturally active but kept repressed by other genes in the cascade (with small dissociation constants, ${\theta }_s$, indicating strong repression; Equations 17.1, 17.2). The only exception is that a gene weakly represses the gene preceding it in the cascade (with larger dissociation constants, ${\theta }_w$, indicating weaker repression) (Figure 3H). When the cascade is initialized by activating the first gene, a domino effect occurs where the next gene represses the previous one, leading to the de-repression of the subsequent gene, and so forth [Figure 3H’; see Computer Simulation 17 (Supplementary Text S2)].

For $m=1$:

$\frac{d{G}_1}{dt}={\displaystyle \prod _{\forall i}^{i\ne 1}}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}-\lambda G_1(17.1)$

For $m\ne 1$:

$\frac{d{G}_m}{dt}=\frac{1}{1+{\left(\frac{G_{m-1}}{{\theta }_w}\right)}^n}\times {\displaystyle \prod _{\forall i}^{i\ne m,m-1}}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}-\lambda G_m(17.2)$

2.2.4.3 Oscillators

The basic logic of an oscillator, which is based on mutually repressive links and de-repressive interactions, is similar to that of a genetic cascade. However, in this case, the last gene in the sequence does not repress the first gene, allowing the system to cycle continuously (Equations 18.1, 18.2 [where N is the total number of genes in the GRN); Figure 3I, I’; see Computer Simulation 18 (Supplementary Text S2)].

For $m=1$:

$\frac{d{G}_1}{dt}={\displaystyle \prod _{\forall i}^{i\ne 1,N}}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}-\lambda G_1(18.1)$

For $m\ne 1$:

$\frac{d{G}_m}{dt}=\frac{1}{1+{\left(\frac{G_{m-1}}{{\theta }_w}\right)}^n}\times {\displaystyle \prod _{\forall i}^{i\ne m,m-1}}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}-\lambda G_m(18.2)$

In our modeling framework of GRNs presented above, the cis-regulatory regions of constituent genes primarily encoded AND and/or OR regulatory logics of activators and repressors. However, developmental genes might exhibit more complex regulatory logic (Kim et al., 2022; Setty et al., 2003; Mayo et al., 2006; Yuh et al., 1998), with TFs acting as activators in some conditions and repressors in others (Hanna-Rose et al., 1997; Zuo et al., 1991; Ilsley et al., 2013). Given this complexity, one might question the value of the GRN models presented here—as well as many GRN implementations found in the literature—especially that their functionality often relies heavily on parameters that are experimentally challenging to determine. Some modeling approaches sidestep explicit GRN interactions altogether, focusing instead on the overall dynamic behavior of the network using simple equations to characterize overall functionality, like homeostasis or oscillations (Jaeger and Monk, 2014; Jutras-Dubé et al., 2020; Corson and Siggia, 2012), an approach that is often called “geometric” modeling. While such approaches have their merits, GRN modeling remains valuable despite its unverified assumptions. Although the detailed wiring may not be entirely accurate, the overall regulatory logic of a GRN can capture essential system dynamics. For instance, as discussed, developmental genetic cascades can follow different overall regulatory logics—either a relay logic, where each gene activates the next, or a logic that depends of de-repressions, where downregulation of one gene releases the next gene within the cascade. These two realizations may perform similarly under normal conditions but differ in robustness and responses to genetic perturbations (Averbukh et al., 2018). Pure dynamical modeling cannot capture these differences, whereas traditional GRN models, despite their limitations, offer insights into critical aspects of the network’s regulatory logic.

2.3 Cellular differentiation of GRNs: space enters the scene

So far, we have discussed GRNs that function within a single cell. In development, although cells within a tissue share the same genes and GRNs (due to having the same genome), they are regulated differently, allowing each cell to express a unique set of TFs and follow distinct developmental trajectories. This manifests as gene expression patterns across the tissue, where different subsets of cells express different genes.

While such expression patterns could, in theory, take any form, they usually adopt specific, stereotyped patterns. Two major classes of patterns typically observed during early developmental stages are periodic and non-periodic patterns (Figures 4A, A’; shown for simplicity as a one-dimensional row of cells) (Negrete and Oates, 2021; Diaz-Cuadros et al., 2021). Periodic patterns mediate the division of tissue into serial structures, such as demarcating the segmented structures along the AP axis in arthropods and vertebrates. Non-periodic patterns, which are more common, mediate the division of tissue into distinct cell fates (regionalization). This differential expression of genes across a tissue results from regulating the same GRN differently in various cells (Figure 4B).

Figure 4

Figure 4. Introduction of spatial patterns in gene regulatory networks. (A, A′) Illustration of non-periodic and periodic patterns in a tissue. The schematic shows a linear arrangement of cells exhibiting non-periodic (distinct regions) and periodic (repeating units) gene expression patterns. (B) Conceptual representation of cells within a tissue sharing the same GRN but being regulated differently. Different regulatory inputs lead to distinct expression patterns despite identical GRNs in each cell. (C) Formation of a morphogen gradient by an organizer. An organizer at one end of the tissue secretes a ligand that diffuses to form a gradient, providing positional information to cells. (C′) Patterning followed by tissue elongation. The diagram illustrates how a tissue is patterned by a morphogen gradient before undergoing elongation, resulting in the expansion of patterned regions. (D) Concurrent tissue elongation and patterning with a retracting gradient. A moving morphogen gradient, exposing cells to changing morphogen concentrations over time, patterns the tissue as it elongates. Some elements of this figure was created by BioRender.com.

To regulate cells differently across a tissue, a morphogen gradient is often employed (Briscoe and Small, 2015; Rogers and Schier, 2011; Simsek and Özbudak, 2022). An organizer (usually positioned at one extremity of the tissue to be patterned) secretes a ligand that diffuses and forms a morphogen gradient (Figure 4C). Different concentrations of the gradient set distinct developmental paths for cells along the tissue. However, tissue elongation—whether due to growth or convergent extension—often accompanies embryonic patterning. In some cases, tissue elongation occurs after patterning (Figure 4C’), while in others, it takes place concurrently with patterning (Figure 4D).

Patterning a tissue during elongation is mechanistically distinct from patterning a non-elongating tissue (or a tissue patterned after elongation), as it often involves a dynamically retracting morphogen gradient. As a result, cells are exposed to changing concentrations of the morphogen over time (Figure 4D). In summary, embryonic patterning can be categorized based on the nature of the pattern (periodic vs non-periodic) and whether the tissue is elongating or non-elongating (Diaz-Cuadros et al., 2021). These developmental processes are explored further in the Results section.

The mechanisms that drive morphogen gradient formation range from passive diffusion to highly regulated ligand transport across a tissue, which are covered extensively in other reviews (Rogers and Schier, 2011; Zhu et al., 2020; Bollenbach et al., 2008; Stapornwongkul and Vincent, 2021; Lord et al., 2021) and are beyond the scope of this study. Here, we use a simplified formulation for the gradient, independent of the mechanism that generates it. For the case of a non-retracting gradient in a non-elongating tissue, we use the following sigmoid function, where $x_0$ specifies the inflection point of the gradient, and $k$ denotes the gradient slope:

$M\left(x\right)=\frac{1}{1+e^{-k\left(x-x_0\right)}}(19.1)$

For the case of a retracting gradient, we employ the following moving sigmoid function, where $v$ is the velocity of gradient retraction:

$M\left(x,t\right)=\frac{1}{1+e^{-k\left[x-\left(x_0\left(0\right)-vt\right)\right]}}(19.2)$

3 Results

In this section, we will discuss two of the most prominent models of embryonic pattern formation: the French Flag model (Wolpert, 1969; Green and Sharpe, 2015) and the Temporal Patterning (or Speed Regulation) model (Diaz-Cuadros et al., 2021; Zhu et al., 2017; Boos et al., 2018; Rudolf et al., 2020; El-Sherif et al., 2012; El-Sherif et al., 2014; Ebisuya and Briscoe, 2018; Dessaud et al., 2007). These models are phenomenological in nature, focusing on descriptive aspects and leaving out the molecular details of their genetic or molecular realization. We will explore four variations of each model, addressing the four common scenarios of embryonic patterning (Negrete and Oates, 2021): a non-periodic pattern in a non-elongating tissue (Wolpert, 1969), a non-periodic pattern in an elongating tissue (Green and Sharpe, 2015), a periodic pattern in a non-elongating tissue, and (Müller and El-Sherif, 2020) a periodic pattern in an elongating tissue.

Following this, we will apply the gene regulation models previously discussed to investigate different molecular and genetic realizations of these two phenomenological models of embryonic pattern formation. Although the reaction-diffusion model is another significant framework for understanding pattern formation, little is known about its regulation at the molecular level, and it will not be covered in this study.

3.1 The french flag model

The core mechanism of the French Flag model (Wolpert, 1969; Green and Sharpe, 2015; Diaz-Cuadros et al., 2021) is that different concentrations of a morphogen gradient activate different genes or cellular states (Figure 5A). This mechanism is ideal for mediating non-periodic patterns in non-elongating tissue (Figure 5A, left panel), where certain ranges of the morphogen gradient activate specific genes. However, a more challenging scenario arises when applying the French Flag model to pattern an elongating tissue during elongation, as the model assumes a stable gradient to set precise boundaries between gene expression domains. To address this, a retracting gradient (M’; black in Figure 5A, right panel; Equation 20.1) activates a gene (M; grey in Figure 5A, right panel; Equation 20.2) with a slow decay rate (ideally $\lambda =0$), forming a stable long-range morphogen gradient (grey in Figure 5A) from a retracting short-range gradient (black in Figure 5A). This stable gradient then mediates the division of the tissue into distinct gene expression domains. Extending the French Flag model to periodic patterns (Figure 5B) is conceptually straightforward: alternating ranges of morphogen gradient concentrations turn a gene on and off. However, achieving this at the molecular level is more complex as will be discussed later.

$M^{\prime }=\frac{1}{1+e^{-k\left[x-\left(x_0\left(0\right)-vt\right)\right]}}(20.1)$

$\frac{dM}{dt}=\frac{M^{\prime }}{1+M^{\prime }}(20.2)$

Figure 5

Figure 5. The French Flag Model and its gene regulatory network realizations. (A) French Flag model for non-periodic patterning. Left panel: A stable morphogen gradient divides the tissue into distinct gene expression domains (blue, red, green). Right panel: In an elongating tissue, a retracting gradient (black) induces a stable gradient with slow decay (gray) to pattern the tissue. (B) French Flag model applied to periodic patterning. Alternating concentrations of a morphogen gradient lead to the activation and repression of genes, creating repeating patterns across the tissue. (C) A GRN realization of the French Flag, where patterning genes (shown G1-G4) are activated or repressed by different thresholds of the morphogen gradient M (gray), resulting in specific expression domains. (C′) Simulation of GRN realization shown in (C). This simulation corresponds to Computer Simulation 19 (Supplementary Text S2). (D) Gene expression boundaries defined by two morphogen gradients. Gradients M1 (dark gray) and M2 (light gray) overlap to establish precise expression domains for genes G1 (blue) and G2 (red). (D′) Simulation of GRN realization shown in (D). This simulation corresponds to Computer Simulation 20 (Supplementary Text S2). (E) GRN configuration involving cross-regulatory interactions (nested feedforward loops). The morphogen gradient M activates genes G1, G2, and G3 at different thresholds, while cross-repressions between them refine their expression domains. (E′) Simulation of GRN realization shown in (E). This simulation corresponds to Computer Simulation 21 (Supplementary Text S2). (F) The periodic expression of the gene eve in Drosophila is the result of the additive effect of multiple stripe-specific enhancers, each mediating one or two stripes. The logic of each enhancer is mediated by the multiplicative effect of several TF factors. Shown is the regulatory logic of enhancer 3 + 7 (mediating the expression of eve stripes 3 and 7), and enhancer 4 + 6 (mediating the expression of eve stripes 4 and 6). Both enhancers mediate the formation of stripes by integrating inputs from gap gene morphogens (mainly hb and kni). (G) A closer look at the regulatory logic of the eve 3 + 7 and 4 + 6 enhancers in Drosophila given the gene expression patterns of gap genes hb and kni. Harboring strong binding sites (thick line) for kni (green) but weak bind sites (thin line) for hb (blue), enhancer 3 + 7 place stripes 3 and 7 away from kni expression and close to abutting hb expressions. Harboring weak binding sites for kni but strong bind sites for hb, enhancer 4 + 6 place stripes 4 and 6 closer to central kni expression and away from abutting hb expressions. (H) Additive effect of separate enhancers on gene expression. Gene G is expressed in two domains, each regulated by a distinct enhancer (E1 and E2). The overall expression pattern is the sum of both enhancer activities. The same basic principle can be extended to form periodic pattern of arbitrary repetitions. (H′) Simulation of GRN realization shown in (H). This simulation corresponds to Computer Simulation 22 (Supplementary Text S2). Some elements of this figure was created by BioRender.com.

It is important to note that the French Flag model is a purely phenomenological, descriptive model and does not specify how it is realized at the molecular or genetic level. Here, we will explore potential molecular and GRN realizations of this model. Despite the simplicity of the French Flag model, realizing it at the GRN level is challenging, as some genes must be both activated and repressed by the same morphogen gradient. For instance, in Figure 5A, the red gene is activated by moderately high levels of the morphogen gradient (grey) but repressed at very high concentrations. This implies that the red gene responds to two different thresholds, T1 and T2, of the same morphogen gradient. None of the transcription models discussed so far can accommodate this feature for a single gene. While this phenomenon has been experimentally observed (Zuo et al., 1991), it has not yet been fully mechanistically elucidated. For now, we will limit ourselves to simpler models where each gene responds to a single threshold of the morphogen.

We will start by reviewing potential GRN realizations of the French Flag model through a basic gene expression pattern where genes are expressed in bands extending to the tissue ends (Figures 5C, C’). This simple pattern can be realized by adjusting the activation and repression thresholds of patterning genes. For example, in Figure 5C, the blue gene (G1) is strongly activated by the morphogen gradient [M, grey in Figures 5C, C’; see Computer Simulation 19 (Supplementary Text S2)], extending its expression towards the lower end of the gradient. The red gene (G2), on the other hand, is weakly activated by M. This can be achieved by setting G1 to have a low activation threshold (low $K_{G1M}$, indicating strong binding sites for M at G1’s Cis-Regulatory Module (CRM)) and G2 to have a high activation threshold (high $K_{G2M}$, indicating weaker binding sites for M at G2’s CRM). The expression domains of the green and orange genes (G3 and G4, respectively) are formed by having them repressed by M—strongly in the case of G3 (low $K_{G3M}$) and weakly in the case of G4 (high $K_{G4M}$) (Equation 22).

${\left.\frac{d{G}_i}{dt}\right|}_{i=1,2}=\frac{{\left(\frac{M}{K_{GiM}}\right)}^n}{1+{\left(\frac{M}{K_{GiM}}\right)}^n}-\lambda G_i,{\left.\frac{d{G}_i}{dt}\right|}_{i=3,4}=\frac{1}{1+{\left(\frac{M}{K_{GiM}}\right)}^n}-\lambda G_i(22)$

To create gene expression domains that do not extend to the boundaries of the morphogen gradient, two morphogen gradients can be employed. An example of this is shown in Figure 5D, D’ [see Computer Simulation 20; (Supplementary Text S2)], where two gradients (M1 and M2, represented by light and dark grey, respectively) together define the boundaries of the blue and red genes (G1 and G2). In this scenario, G1 is strongly repressed by M1 but weakly repressed by M2 (i.e., small $K_{G1M1}$ and large $K_{G1M2}$), while G2 is weakly repressed by M1 and strongly repressed by M2 (i.e., large $K_{G2M1}$ and small $K_{G2M2}$) (Equation 23). This approach can be easily generalized to multiple genes by adjusting the repression thresholds of the two gradients for each gene.

$\frac{d{G}_i}{dt}=\frac{1}{1+{\left(\frac{M1}{K_{GiM1}}\right)}^n}·\frac{1}{1+{\left(\frac{M2}{K_{GiM2}}\right)}^n}-\lambda G_i(23)$

To generate multiple gene expression domains from a single morphogen gradient, cross-regulatory interactions between genes are required. The specific GRN wiring can vary, but one straightforward realization is shown in Figure 5E. In this configuration, the morphogen gradient M activates patterning genes (G1, G2, and G3) at different thresholds. Without cross-regulatory interactions, this would result in broad expression domains, extending from the high concentration of M to the respective activation thresholds. Repression from other genes is needed to refine these domains. For instance, G1 is repressed by both G2 and G3, restricting its expression to the lower end of the gradient (see Computer Simulation 21 and its simulation in Figure 5E, E’; Equations 24.1–24.3; Note that the circuit depicted in Figure 5E is for a 3 genes realization of this regulatory scheme, whereas Simulation 21; Figure 5E’ are for a 5-genes realization)

$\frac{d{G}_i}{dt}=\frac{{\left(\frac{M}{K_{GiM}}\right)}^n}{1+{\left(\frac{M}{K_{GiM}}\right)}^n}·\frac{1}{1+{\left(\frac{G_{i+1}}{{\theta }_w}\right)}^n}·{\displaystyle \prod _{j=i+2}^N}\frac{1}{1+{\left(\frac{G_j}{{\theta }_s}\right)}^n}-\lambda G_i(24.1)$

For $i=N-1$:

$\frac{d{G}_i}{dt}=\frac{{\left(\frac{M}{K_{GiM}}\right)}^n}{1+{\left(\frac{M}{K_{GiM}}\right)}^n}·\frac{1}{1+{\left(\frac{G_{i+1}}{{\theta }_w}\right)}^n}-\lambda G_i(24.2)$

For $i=N$:

$\frac{d{G}_i}{dt}=\frac{{\left(\frac{M}{K_{GiM}}\right)}^n}{1+{\left(\frac{M}{K_{GiM}}\right)}^n}-\lambda G_i(24.3)$

This GRN configuration can be seen as a nested feedforward loop motifs (Mangan and Alon, 2003), where the morphogen gradient regulates each gene via two parallel paths: one direct, and one indirect through an intermediate gene. Such configuration was shown to be effective in generating single or multiple gene expression domains (Ishihara et al., 2005). Other GRN configurations are also possible. For instance, in the AC-DC GRN motif (Balaskas et al., 2012; Perez-Carrasco et al., 2018; Panovska-Griffiths et al., 2013), the morphogen gradient thresholds emerge from cross-regulation between genes, rather than being dictated by the gradient itself.

So far, we have discussed the formation of non-periodic patterns, noting that their GRN realization is more challenging than the simplicity of the phenomenological French Flag model suggests. This difficulty primarily stems from the challenge of making a single gene respond to two different thresholds of the same morphogen gradient, which we addressed using either two morphogen gradients or cross-regulatory interactions. Applying the French Flag model to periodic patterns is even more complex, as genes must respond to multiple thresholds, with the morphogen gradient alternately activating and repressing the same gene. One solution is similar to the non-periodic case: decomposing the periodic pattern into several non-periodic patterns, each mediated separately by two morphogen gradients. The final periodic pattern is the sum of all sub-patterns.

This method of generating periodic patterns has been extensively studied in the context of AP patterning in the early Drosophila embryo (Figures 5F, G), providing insights into the organization of cis-regulatory modules during development (Schroeder et al., 2011). A notable example is the expression of the even-skipped (eve) gene, which forms seven distinct stripes in the early Drosophila embryo (Schroeder et al., 2011; Goto et al., 1989). Rather than arising from a periodic process like an oscillator or reaction-diffusion mechanism (Landge et al., 2020; Meinhardt et al., 1991), this pattern results from decomposing the periodic pattern into several non-periodic sub-patterns, each composed of a single stripe or pair of stripes (Clark et al., 2019; Diaz-Cuadros et al., 2021; Akam, 1989; Lynch et al., 2012). Each sub-pattern is mediated by a dedicated enhancer, and the final expression of eve is the sum of all enhancer activities. The regulatory logic of each enhancer is primarily determined by the multiplicative (AND) interactions of repressing TFs known as gap genes, which form (localized) morphogen gradients in the early Drosophila embryo. Figures 5F, G illustrates the regulatory logic of two enhancers (Clyde et al., 2003): the 3 + 7 enhancer (which mediates the formation of the 3rd and 7th stripes) and the 4 + 6 enhancer (which mediates the formation of the 4th and 6th stripes). The positioning of these stripes is determined by the binding strengths of the gap gene morphogens hunchback (hb) and knirps (kni) to the 3 + 7 and 4 + 6 enhancers along the AP axis.

This example, along with others, demonstrates that an additive OR logic is mediated by using separate enhancers spaced far enough to avoid interference between their activities, while multiplicative AND logic is mediated by placing TF binding sites close to each other within the same enhancer.

A simple example of modeling this mechanism is shown in Figure 5H, H’, where gene G is expressed in two domains, each mediated by the additive effects of two separate enhancers: Enhancer 1 (E1) and Enhancer 2 (E2). E1 is strongly repressed by morphogen 1 (M1) and weakly repressed by morphogen 2 (M2), while E2 is weakly repressed by M1 and strongly repressed by M2 (see Computer Simulation 22 and its simulation in Figure 5H’; Equations 25.1, 25.2).

$E_i=\frac{{\left(\frac{M1}{K_{EiM1}}\right)}^n}{1+{\left(\frac{M1}{K_{EiM1}}\right)}^n}·\frac{{\left(\frac{M2}{K_{EiM2}}\right)}^n}{1+{\left(\frac{M2}{K_{EiM2}}\right)}^n}-\lambda G_i(25.1)$

$G=E_1+E_2-\lambda G(25.2)$

In the GRN realizations of the French Flag model presented above, we considered only the case of non-elongating tissues. The extension to the case to elongating tissues, where both M and M′ gradients are used (right panels in Figures 5A, B) is straight forward (Boos et al., 2018).

3.2 Temporal patterning or the speed regulation model

In temporal patterning, a temporal sequence is translated into a spatial pattern rather than directly mediating a spatial pattern (Diaz-Cuadros et al., 2021; Durston et al., 2012; Palmeirim et al., 1997). This mode of patterning, although less intuitive, is a standard mechanism during neurogenesis (Kohwi and Doe, 2013; Naidu et al., 2020; Filippopoulou and Konstantinides, 2024) in many animals and at various stages of development. Recent studies have shown it is involved in many other tissues as well, such as segmentation and regionalization of the AP axis in short-germ insects (El-Sherif et al., 2012; Sarrazin et al., 2012) and vertebrates (Palmeirim et al., 1997; Oates et al., 2012; Pourquié, 2003), ventral neural tube patterning (Dessaud et al., 2007), and limb bud development (McQueen and Towers, 2020; Roensch et al., 2013). This raises the possibility that temporal patterning may be the default mechanism for embryonic pattern formation. While the reason for this preference is not fully understood, it has been suggested that this mechanism may offer greater robustness compared to alternatives, such as the French Flag model (Jutras-Dubé et al., 2020; El-Sherif et al., 2014; Balaskas et al., 2012). It has also been hypothesized that since animals evolved from single-celled organisms, many evolutionarily conserved GRNs are temporal in nature, and multicellularity may have evolved by translating these conserved temporal patterns into spatial ones (Filippopoulou and Konstantinides, 2024).

The Speed Regulation (SR) model (Zhu et al., 2017) was recently proposed as a unifying mechanism for various types of temporal patterning observed in different tissues. The SR model synthesizes several patterning schemes that all share the core feature of a morphogen gradient modulating the speed of a temporal sequence, whether periodic or non-periodic. In the non-periodic version of the SR model, each cell within a tissue progresses through successive states (depicted in different colors in Figure 6A), with each state defined by the expression of one or more genes. The speed of transitions between these states is controlled by a molecular factor (referred to as the “speed regulator,” shown in grey at the top of Figure 6A). At very low or zero concentrations of the speed regulator, the transitions slow down to the point where the states are indefinitely stabilized. When a group of cells experiences a gradient of the speed regulator (Figure 6A’; left panel: ‘for non-elongating tissues’), all cells transition through successive states, but at progressively slower speeds as they encounter lower concentrations of the gradient. This creates the illusion of cellular states propagating as waves from high to low gradient concentrations. These waves, known as “kinematic” or “pseudo-waves,” do not rely on diffusion or cell-to-cell communication (BECK and VÁRADI, 1972; Rohde et al., 2024). We refer to this version of the model as “gradient-based speed regulation,” which is particularly suited for patterning non-elongating tissues (Figure 6A’, left). The model can also pattern elongating tissues if the gradient retracts as a wavefront, a process we call “wavefront-based speed regulation” (Figure 6A’, right). The SR model can also generate periodic structures if the sequential gene activation process is driven by a biological clock instead of a sequential transition of states (Diaz-Cuadros et al., 2021). Notably, in the wavefront-based version of SR, if the wavefront is in the form of a tapered gradient (a superposition of the gradient-based and wavefront-based SR models), kinematic waves will propagate from high to low concentrations of the gradient, in the opposite direction of wavefront retraction, as the tissue elongates.

Figure 6

Figure 6. Temporal Patterning or the Speed Regulation Model. (A) Schematic of the temporal sequence of cellular states in temporal patterning. Cells progress through successive states (depicted in different colors), with transitions regulated by a speed regulator (gray). (A′) Gradient-based and wavefront-based speed regulation in tissue patterning. Left: In a non-elongating tissue, a morphogen gradient modulates the speed of transitions between cellular states, resulting in the induction of non-periodic waves that eventually stabilize into a non-periodic pattern. Right: In an elongating tissue, a retracting wavefront regulates the speed of state transitions, resulting in the formation of a non-periodic pattern. (B) Realization of the Speed Regulation model by modulating transcription and decay rates. The morphogen gradient jointly affects the overall transcription and gene product degradation rates to control the timing of gene expression. (B′) Simulation of the GRN realization in (B). This simulation corresponds to Computer Simulation 23 (Supplementary Text S2). (C) Enhancer switching model as a GRN realization of the Speed Regulation model. Each patterning gene is regulated by the additive sum of a dynamic GRN (driving sequential activities) and a static GRN (stabilizing expression) through two separate enhancers. The speed regulator modulates the balance between these GRNs. (C′) Simulation of the GRN realization in (C). This simulation corresponds to Computer Simulation 24 (Supplementary Text S2). Some elements of this figure was created by BioRender.com.

Like the French Flag (FF) model, the SR model is phenomenological and does not propose a molecular mechanism or GRN realization that can transform temporal sequences (or oscillations) into patterns. A straightforward realization of the model is for a morphogen gradient (the speed regulator) to modulate the overall transcription rate, including gene product decay rates (Boos et al., 2018) (Figure 6B) M2 [see Computer Simulation 23 (Supplementary Text S2)] and its simulation in Figure 6B’; Note that the circuit depicted in Figure 6B is for a 3 genes realization of this regulatory scheme, whereas Simulation 23 and Figure 5B’ are for a 5-genes realization).

For $j=1$:

$G_j=\frac{M}{1+M}\times \left({\displaystyle \prod _{i=2}^J}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}-\lambda G_j\right)(26.1)$

For $j\ne 1$:

$G_j=\frac{M}{1+M}\times \left(\frac{1}{1+{\left(\frac{G_{m-1}}{{\theta }_w}\right)}^n}{\displaystyle \prod _{i=1,i\ne j,j-1}^J}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}-\lambda G_j\right)(26.2)$

While this approach works well, it requires that the decay rates of all gene products involved in patterning be regulated by the morphogen gradient in proportion to the regulation of their transcription. Achieving this at the molecular level may be unfeasible.

An alternative mechanism, recently proposed to modulate the timing of GRNs, is the Enhancer Switching Model (Diaz-Cuadros et al., 2021; Zhu et al., 2017; Mau et al., 2023). This model posits that each patterning gene is simultaneously wired into two GRNs (Figure 6C): (i) a dynamic GRN that drives periodic or sequential gene activities, and (ii) a static GRN that stabilizes gene expression patterns. The concentration of the speed regulator (shown in grey in Figure 6C) activates the dynamic GRN while repressing the static GRN, thus balancing the contribution of each GRN to the overall dynamics and consequently regulating the speed of gene expression. At high concentrations of the speed regulator, the dynamic GRN is dominant, facilitating rapid oscillations or sequential gene activities. At low concentrations, the static GRN dominates, resulting in slower oscillations or sequential gene activities. Each gene is connected to these two GRNs through two enhancers: (i) a dynamic enhancer that encodes the wiring of the gene within the dynamic GRN, and (ii) a static enhancer that encodes the wiring of the gene within the static GRN (Figure 1D) (Equations 27.1–27.4; see Computer Simulation 24 (Supplementary Text S2) and its simulation in Figure 6C’; Note that the circuit depicted in Figure 6C is for a 3 genes realization of this regulatory scheme, whereas Simulation 24 and Figure 5C’’ are for a 5-genes realization).

$G_j=k_{td}{E}_{dj}+k_{ts}{E}_{sj}-\lambda G_j(27.1)$

For the dynamic enhancer, $E_{dj}$:

For $j=1$:

$E_{dj}=\frac{M}{1+M}\times {\displaystyle \prod _{\forall i}^{i\ne 1}}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}(27.2)$

For $j\ne 1$:

$E_{dj}=\frac{M}{1+M}\times \frac{1}{1+{\left(\frac{G_{j-1}}{{\theta }_w}\right)}^n}\times {\displaystyle \prod _{\forall i}^{i\ne j,j-1}}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}(27.3)$

For the static enhancer, $E_{sj}$:

$E_{sj}=\frac{1}{1+M}\times {\displaystyle \prod _{i=1,i\ne j}^J}\frac{1}{1+{\left(\frac{G_i}{{\theta }_s}\right)}^n}(27.4)$

The Enhancer Switching model has been proposed as the molecular mechanism underlying AP patterning in the intermediate-germ insect Tribolium castaneum, with partial experimental support (Boos et al., 2018; Mau et al., 2023). This model also aligns with recent findings on AP patterning in Drosophila (El-Sherif and Levine, 2016; Clark and Akam, 2016; Koromila et al., 2020; Soluri et al., 2020). However, further rigorous validation is needed to confirm its broader applicability.

It is worth noting here that timing control observed in the enhancer switching model arises from its influence at the GRN level rather than at the molecular DNA binding timescale. This occurs because the static module, a multi-stable network, introduces attractor points in the phase space of the dynamic module, which governs cyclical or sequential expression. The relative strength between the dynamic and static modules, determined by the speed regulator, dictates the extent to which cyclical or sequential attractors are influenced by fixed-point attractors, effectively slowing down their progression (Jutras-Dubé et al., 2020).

Many of the presented models exhibit robustness to variations in key parameters. For example, the Enhancer Switching model continues to produce largely accurate gene expression patterns despite one or more fold changes in the DNA binding dissociation constants (${\theta }_s$) of strong repressors within the Enhancer Switching GRN, modifications in gene product decay rates (λ), or variations in the relative strength of static and dynamic modules (k_ts/k_td) (Supplementary Figure S3).

3.3 Case study: impact of enhancer integration mechanisms on embryonic patterning

So far, we have discussed various embryonic pattern formation mechanisms using basic models of transcription. To that end, we employed simple toy models at each level of the modeling process. One of the main advantages of this approach is that it allows us to examine the effect of fundamental transcriptional and gene regulatory mechanisms on the overall performance of embryonic patterning systems. In this section, we will explore a specific example to demonstrate this. We will examine how the integration of multiple enhancer activities influences the performance of the Enhancer Switching model.

In our modeling of the enhancer switching mechanism (Figure 6C; Equations 27.1–27.4), we used a weighted sum of the two enhancer activities (Equation 27.1). The physical interpretation of this relationship is that each enhancer can transcribe the gene independently (Figure 7A). In other models of enhancer function, however, enhancers compete for binding to the promoter, and the enhancer that wins the competition sets the transcriptional state of the promoter (Equation 28; Figure 7B) (Bothma et al., 2015). In the first scenario (Figure 7A; Equation 27.1), each enhancer determines its own final transcriptional rate (denoted as $k_{td}$ and $k_{ts}$ in Equation 27.1), whereas in the second scenario, the normalized activities of competing enhancers set the transcriptional state of the promoter, but the maximum final transcriptional rate (denoted as $k_t$ in Equation 28) is set by the promoter itself (Figure 7B).

$G_j=k_t\left(\frac{k_{ed}{E}_{dj}+k_{es}{E}_{sj}}{1+k_{ed}{E}_{dj}+k_{es}{E}_{sj}}\right)-\lambda G_j(28)$

Figure 7

Figure 7. Case study: Impact of enhancer integration mechanisms on embryonic patterning. (A) Additive enhancer integration model. Multiple enhancers independently contribute to the transcription rate of a gene. The total transcription is the weighted sum of the individual enhancer activities. (A′) Simulation of gene expression with additive enhancer integration under normal conditions. The performance of the enhancer switching model is shown, demonstrating effective temporal patterning. This simulation corresponds to Computer Simulation 25 (Supplementary Text S2) with mut_flag = 0 (wild-type). (A″) Impact of an overactive enhancer in the additive model. The simulation illustrates how an overactive enhancer disrupts gene expression, leading to aberrant patterning and loss of GRN function. This simulation corresponds to Computer Simulation 25 (Supplementary Text S2) with mut_flag = 1 (overactive dynamic enhancer of Gene 3). (B) Enhancer competition integration model. Enhancers compete for binding to the promoter of a gene. The promoter’s transcriptional state is determined by the enhancer that successfully binds, normalizing the transcriptional output. (B′) Simulation of gene expression with enhancer competition integration under normal conditions. The model shows comparable performance to the additive case, maintaining proper temporal patterning. This simulation corresponds to Computer Simulation 26 (Supplementary Text S2) with mut_flag = 0 (wild-type). (B″) Mitigation of overactive enhancer effects in the competition model. The simulation demonstrates that the competitive mechanism normalizes transcription levels, preserving GRN function despite the presence of an overactive enhancer. This simulation corresponds to Computer Simulation 26 (Supplementary Text S2) with mut_flag = 1 (overactive dynamic enhancer of Gene 3).

The normalization effect of the enhancer competition scenario may be a crucial mechanism for integrating enhancer activities. The overall transcription level of a gene is a key parameter in any GRN, as the activation thresholds of genes within a GRN depend on the average transcription rates of the genes regulating them. Enhancers are rapidly evolving genetic elements and are the primary drivers of novelty in animal evolution. However, this also introduces a risk, as allowing a newly evolved enhancer to take control of a key transcriptional parameter, such as the overall average transcription rate, could disrupt normal function. Alternatively, in the enhancer competition model, enhancers first compete to activate the promoter, with the promoter acting as the final arbiter for the transcription rate.

To investigate this further, we simulated our enhancer switching model for both the additive case [Figures 7A, A’; see Computer Simulation 25 (Supplementary Text S2)] and the enhancer competition case [Figures 7B, B’; see Computer Simulation 26 (Supplementary Text S2)]. Both scenarios resulted in comparable performance (Figure 7A’, B’). We then examined the effect of an overactive enhancer (Figure 7A”, B”). In the additive case, the overactive enhancer disrupted the overall function of the GRN (Figure 7A”), whereas in the enhancer competition case, the normalization effect mitigated the impact of the overactive enhancer, preserving the overall performance of the GRN (Figure 7B”).

4 Discussion

In this paper, we developed a comprehensive framework to model the emergence of embryonic patterning, linking molecular gene regulation to tissue-level organization. We began by modeling transcription at the single-gene level using basic chemical reaction models and extended this to model GRNs that govern specific cellular functions. We then introduced phenomenological models of embryonic pattern formation, such as the French Flag model and Speed Regulation models, integrating these with molecular and GRN realizations. To facilitate understanding and application of our models, we accompanied our mathematical framework with computer simulations, providing intuitive and simple code for each model. Through our case study on enhancer integration, we demonstrated that a two-step gene regulation strategy—enhancer activation followed by competitive integration at the promoter—ensures robustness and evolvability in gene expression patterns, emphasizing the adaptability of eukaryotic transcriptional regulation.

Modeling gene regulation during embryonic development has traditionally followed two distinct strategies. The first involves using finely tuned GRN models to fit large datasets, which we refer to as the “simulation approach” (Jaeger et al., 2004; Verd et al., 2018; Manu et al., 2009). The second relies on simplified, abstract models that capture the overall behavior of the observed data, which we call the “toy model approach” (Alon, 2006; Zhu et al., 2017; Balaskas et al., 2012; François et al., 2007; François and Siggia, 2010). In our view, the simulation approach is only appropriate when detailed biochemical data are available, which is often not the case. By contrast, while the toy model approach is inherently oversimplified, it offers the advantage of conceptual clarity. However, both approaches share a significant limitation in modeling gene regulation during development: they rely on gene regulatory models that are largely inspired by bacterial systems. These models may not accurately reflect the complexity of eukaryotic transcriptional machinery, particularly in animals. In this paper, we adopted the toy model approach, acknowledging its inherent limitations. However, our aim was to explicitly articulate the various assumptions made when modeling gene regulation, from the molecular scale to the tissue level. This transparency allows for critical examination of these assumptions. The simplicity and low computational cost of the toy model approach—unlike the simulation approach, which requires extensive data collection and optimization—facilitates this exploration. Moreover, we provided accompanying computer simulations with intuitive and simple code, making our models accessible and facilitating their use as educational tools or starting points for further research. In fact, we present our model not as a definitive framework, but as a basis for questioning and refining the assumptions at each step. Our hope is that future work by our group and others will further challenge and improve upon these ideas.

One specific aspect of gene regulation we examine in this paper, which relates to the complex transcriptional machinery of animals, is the impact of multiple enhancer integration on the performance of GRNs during development—an aspect largely overlooked in previous modeling efforts. This represents a small departure from the bacteria-inspired models traditionally used in gene regulation studies, moving instead toward the more intricate transcriptional mechanisms found in animals, particularly during development. However, this is neither the only nor the most critical aspect to reconsider in efforts to model gene regulation more realistically. A widely used assumption, which we also adopted in our modeling, is the simplistic regulatory logic of individual genes within the GRN, typically modeled with basic OR and/or AND relationships. Under this assumption, the complexity of GRN performance stems from the network’s structure rather than from the computational capabilities of individual genes. Yet, it has been observed that single genes can exhibit far more complex regulatory logic (Yuh et al., 1998). For instance, a single TF may act as an activator or a repressor, depending on the concentration of co-regulating TFs (Hanna-Rose et al., 1997). These intricate regulatory interactions are likely evolutionarily conserved and could serve as modular building blocks, akin to Lego pieces, which can be used as is or slightly modified to construct more extensive GRNs.

Another key aspect of gene regulation that warrants further investigation using our modeling framework is how the timing of GRNs is modulated. Recent findings suggest that many phenomena in embryonic patterning arise from the modulation of GRN timing by morphogen gradients, supporting a temporal or speed regulation model. In this paper, we modeled this effect using two distinct hypotheses (shown in Figures 6B, C, respectively). However, given the largely unexplored computational complexity of transcriptional machinery of developmental genes, it is conceivable that additional mechanisms for timing regulation exist. Moreover, recent experimental findings suggest other assumptions in current models of gene regulation may need re-evaluation. For example, the effects of transcriptional bursting (Fukaya et al., 2016; Dar et al., 2012) and the formation of TF clusters or condensates (Cho et al., 2016; Sabari et al., 2018; Hamamoto and Fukaya, 2022) on GRN performance—particularly in embryonic patterning—remain unclear. These newly discovered phenomena could significantly influence how GRNs function.

With so many uncertainties surrounding current gene regulation and GRN modeling approaches, how can we approach the complex task of modeling embryonic patterning? While the answer is not straightforward, careful and minimal use of simple, transparent models—despite their inaccuracies—can yield partial insights. This is achievable as long as we remain aware of the models’ assumptions and are willing to critically evaluate them. In this paper, our goal has been to accomplish exactly that. To support further exploration, we offered a simple and transparent computational framework for modeling gene regulation during embryonic pattern formation, complete with accessible simulations and intuitive code. This enables researchers to test and expand upon our models, promoting a deeper understanding of the mechanisms underlying gene regulation in development.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

JG-G: Writing–original draft, Writing–review and editing. EE-S: Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that Generative AI was used in the creation of this manuscript. To fine-tune the writing.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2025.1522725/full#supplementary-material

SUPPLEMENTARY FIGURE S1 | Modeling the kinetics of simple chemical reactions. (A, A') Simulation of the kinetics of a simple chemical reaction where chemical species A and B react to form chemical species C. The reaction follows the equation, with a reaction rate proportional to the concentrations of A and B. The graph shows the concentrations of A (blue curve), B (red curve), and C (orange curve) over time. As the reaction progresses, the concentrations of reactants A and B decrease while the concentration of product C increases, illustrating the consumption of reactants and formation of product over time. This simulation corresponds to Computer Simulation 1 (Supplementary Text S1). (B, B') Simulation of a reversible chemical reaction where chemical species A and B react to form species C, which can also decompose back into A and B. The reaction follows the equation, incorporating both forward and reverse reactions with rate constants and, respectively. The graph displays the concentrations of A (blue curve), B (red curve), and C (orange curve) over time. The system approaches equilibrium where the rates of the forward and reverse reactions balance, resulting in constant concentrations of A, B, and C. This simulation demonstrates how reversible reactions reach equilibrium within a given timeframe and corresponds to Computer Simulation 2 (Supplementary Text S2).

SUPPLEMENTARY FIGURE S2 | Modeling gene activation dynamics and the effects of transcription factor binding kinetics. (A) Simulation of gene activation by a single activator TF with slow binding dynamics (low and). The graph depicts the concentrations of X (blue curve), D1 (orange curve), and Y (red curve) over time. Slower TF binding and unbinding result in delayed formation of the active DNA state and a gradual increase in mRNA levels. This demonstrates that slow TF-DNA binding kinetics can lead to a delayed gene expression response, affecting the timing of downstream cellular processes. This simulation corresponds to Computer Simulation 3 (Supplementary Text S2) with low and values. (B) Simulation of gene activation by a single activator transcription factor (TF) with slow binding dynamics (high and). The model illustrates the interaction between TF X, the active DNA state D1, and the mRNA transcript Y. The graph shows the concentrations of X (blue curve), D1 (orange curve), and Y (red curve) over time. Rapid TF binding and unbinding lead to quick fluctuations in the active DNA state and a prompt increase in mRNA production, resulting in a swift response in gene expression levels. This highlights how fast TF-DNA interactions can affect the timing of gene activation. This simulation corresponds to Computer Simulation 3 (Text S2) with high and values. (C) Simulation of gene activation using Michaelis-Menten kinetics. The model describes how the concentration of mRNA transcript Y changes over time in response to a constant concentration of activator TF X. The graph shows the concentration of Y (orange curve) increasing over time and reaching a steady state. The use of Michaelis-Menten kinetics captures the saturation effect where, at high TF concentrations, the rate of mRNA production approaches a maximum due to the limited number of available binding sites on the DNA. This results in a hyperbolic relationship between TF concentration and gene expression level, illustrating the principles of enzyme kinetics applied to transcriptional activation. This simulation corresponds to Computer Simulation 4 (Supplementary Text S2).

SUPPLEMENTARY FIGURE S3 | Sensitivity analysis of key parameters of the Enhancer Switching model. Study of 0.2x-5x fold changes in key parameters of the Enhancer Switching GRN model: (A) the binding dissociation constants ( ) of strong repressions within the GRN, (B) gene product decay rates (λ), and (C) relative strength of static and dynamic modules (k_ts/k_td).

References

Adelman, K., and Lis, J. T. (2012). Promoter-proximal pausing of RNA polymerase II: emerging roles in metazoans. Nat. Rev. Genet. 13 (10), 720–731. doi:10.1038/nrg3293

PubMed Abstract | CrossRef Full Text | Google Scholar

Akam, M. (1989). Drosophila development: making stripes inelegantly. Nature 341 (6240), 282–283. doi:10.1038/341282a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Alon, U. (2006). An introduction to systems biology: Design principles of biological circuits. illustrated ed. New York: CRC Press.

Google Scholar

Alon, U. (2007). Network motifs: theory and experimental approaches. Nat. Rev. Genet. 8 (6), 450–461. doi:10.1038/nrg2102

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, P. W. (1972). More is different. Science. 177 (4047), 393–396. doi:10.1126/science.177.4047.393

PubMed Abstract | CrossRef Full Text | Google Scholar

Averbukh, I., Lai, S.-L., Doe, C. Q., and Barkai, N. (2018). A repressor-decay timer for robust temporal patterning in embryonic Drosophila neuroblast lineages. eLife 7, e38631. doi:10.7554/eLife.38631

PubMed Abstract | CrossRef Full Text | Google Scholar

Balaskas, N., Ribeiro, A., Panovska, J., Dessaud, E., Sasai, N., Page, K. M., et al. (2012). Gene regulatory logic for reading the Sonic Hedgehog signaling gradient in the vertebrate neural tube. Cell 148 (1–2), 273–284. doi:10.1016/j.cell.2011.10.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Beck, M. T., and Váradi, Z. B. (1972). One, two and three-dimensional spatially periodic chemical reactions. Nature 235, 15–16. doi:10.1038/physci235015a0

CrossRef Full Text | Google Scholar

Bintu, L., Buchler, N. E., Garcia, H. G., Gerland, U., Hwa, T., Kondev, J., et al. (2005a). Transcriptional regulation by the numbers: models. Curr. Opin. Genet. Dev. 15 (2), 116–124. doi:10.1016/j.gde.2005.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bintu, L., Buchler, N. E., Garcia, H. G., Gerland, U., Hwa, T., Kondev, J., et al. (2005b). Transcriptional regulation by the numbers: applications. Curr. Opin. Genet. Dev. 15 (2), 125–135. doi:10.1016/j.gde.2005.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Bollenbach, T., Pantazis, P., Kicheva, A., Bökel, C., González-Gaitán, M., and Jülicher, F. (2008). Precision of the Dpp gradient. Development 135 (6), 1137–1146. doi:10.1242/dev.012062

PubMed Abstract | CrossRef Full Text | Google Scholar

Boos, A., Distler, J., Rudolf, H., Klingler, M., and El-Sherif, E. (2018). A re-inducible gap gene cascade patterns the anterior-posterior axis of insects in a threshold-free fashion. eLife 7, e41208. doi:10.7554/eLife.41208

PubMed Abstract | CrossRef Full Text | Google Scholar

Bothma, J. P., Garcia, H. G., Ng, S., Perry, M. W., Gregor, T., and Levine, M. (2015). Enhancer additivity and non-additivity are determined by enhancer strength in the Drosophila embryo. eLife 12, e07956. doi:10.7554/eLife.07956

PubMed Abstract | CrossRef Full Text | Google Scholar

Briscoe, J., and Small, S. (2015). Morphogen rules: design principles of gradient-mediated embryo patterning. Development 142 (23), 3996–4009. doi:10.1242/dev.129452

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Boyaci, H., and Campbell, E. A. (2021). Diverse and unified mechanisms of transcription initiation in bacteria. Nat. Rev. Microbiol. 19 (2), 95–109. doi:10.1038/s41579-020-00450-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, W.-K., Jayanth, N., English, B. P., Inoue, T., Andrews, J. O., Conway, W., et al. (2016). RNA Polymerase II cluster dynamics predict mRNA output in living cells. eLife 5, e13617. doi:10.7554/eLife.13617

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, E., and Akam, M. (2016). Odd-paired controls frequency doubling in Drosophila segmentation by altering the pair-rule gene regulatory network. eLife 5, e18215. doi:10.7554/eLife.18215

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, E., Peel, A. D., and Akam, M. (2019). Arthropod segmentation. Development 146 (18), dev170480. doi:10.1242/dev.170480

PubMed Abstract | CrossRef Full Text | Google Scholar

Clyde, D. E., Corado, M. S. G., Wu, X., Paré, A., Papatsenko, D., and Small, S. (2003). A self-organizing system of repressor gradients establishes segmental complexity in Drosophila. Nature 426 (6968), 849–853. doi:10.1038/nature02189

PubMed Abstract | CrossRef Full Text | Google Scholar

Corson, F., and Siggia, E. D. (2012). Geometry, epistasis, and developmental patterning. Proc. Natl. Acad. Sci. U. S. A. 109 (15), 5568–5575. doi:10.1073/pnas.1201505109

PubMed Abstract | CrossRef Full Text | Google Scholar

Cramer, P. (2019). Organization and regulation of gene transcription. Nature 573 (7772), 45–54. doi:10.1038/s41586-019-1517-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Dar, R. D., Razooky, B. S., Singh, A., Trimeloni, T. V., McCollum, J. M., Cox, C. D., et al. (2012). Transcriptional burst frequency and burst size are equally modulated across the human genome. Proc. Natl. Acad. Sci. U. S. A. 109 (43), 17454–17459. doi:10.1073/pnas.1213530109

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, E. H. (2010). Emerging properties of animal gene regulatory networks. Nature 468 (7326), 911–920. doi:10.1038/nature09645

PubMed Abstract | CrossRef Full Text | Google Scholar

Dessaud, E., Yang, L. L., Hill, K., Cox, B., Ulloa, F., Ribeiro, A., et al. (2007). Interpretation of the sonic hedgehog morphogen gradient by a temporal adaptation mechanism. Nature 450 (7170), 717–720. doi:10.1038/nature06347

PubMed Abstract | CrossRef Full Text | Google Scholar

Diaz-Cuadros, M., Pourquié, O., and El-Sherif, E. (2021). Patterning with clocks and genetic cascades: segmentation and regionalization of vertebrate versus insect body plans. PLoS Genet. 17 (10), e1009812. doi:10.1371/journal.pgen.1009812

PubMed Abstract | CrossRef Full Text | Google Scholar

Durston, A., Wacker, S., Bardine, N., and Jansen, H. (2012). Time space translation: a hox mechanism for vertebrate a-p patterning. Curr. Genomics 13 (4), 300–307. doi:10.2174/138920212800793375

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebisuya, M., and Briscoe, J. (2018). What does time mean in development? Development 145 (12), dev164368. doi:10.1242/dev.164368

PubMed Abstract | CrossRef Full Text | Google Scholar

Elowitz, M. B., and Leibler, S. (2000). A synthetic oscillatory network of transcriptional regulators. Nature 403 (6767), 335–338. doi:10.1038/35002125

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Sherif, E., Averof, M., and Brown, S. J. (2012). A segmentation clock operating in blastoderm and germband stages of Tribolium development. Development 139 (23), 4341–4346. doi:10.1242/dev.085126

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Sherif, E., and Levine, M. (2016). Shadow enhancers mediate dynamic shifts of gap gene expression in the drosophila embryo. Curr. Biol. 26 (9), 1164–1169. doi:10.1016/j.cub.2016.02.054

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Sherif, E., Zhu, X., Fu, J., and Brown, S. J. (2014). Caudal regulates the spatiotemporal dynamics of pair-rule waves in Tribolium. PLoS Genet. 10 (10), e1004677. doi:10.1371/journal.pgen.1004677

PubMed Abstract | CrossRef Full Text | Google Scholar

Filippopoulou, K., and Konstantinides, N. (2024). Evolution of patterning. FEBS J. 291 (4), 663–671. doi:10.1111/febs.16995

PubMed Abstract | CrossRef Full Text | Google Scholar

François, P., Hakim, V., and Siggia, E. D. (2007). Deriving structure from evolution: metazoan segmentation. Mol. Syst. Biol. 3, 154. doi:10.1038/msb4100192

PubMed Abstract | CrossRef Full Text | Google Scholar

François, P., and Siggia, E. D. (2010). Predicting embryonic patterning using mutual entropy fitness and in silico evolution. Development 137 (14), 2385–2395. doi:10.1242/dev.048033

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuda, N. J., Ardehali, M. B., and Lis, J. T. (2009). Defining mechanisms that regulate RNA polymerase II transcription in vivo. Nature 461 (7261), 186–192. doi:10.1038/nature08449

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukaya, T., Lim, B., and Levine, M. (2016). Enhancer control of transcriptional bursting. Cell 166 (2), 358–368. doi:10.1016/j.cell.2016.05.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Furlong, E. E. M., and Levine, M. (2018). Developmental enhancers and chromosome topology. Science 361 (6409), 1341–1345. doi:10.1126/science.aau0320

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaertner, B., and Zeitlinger, J. (2014). RNA polymerase II pausing during development. Development 141 (6), 1179–1183. doi:10.1242/dev.088492

PubMed Abstract | CrossRef Full Text | Google Scholar

Goto, T., Macdonald, P., and Maniatis, T. (1989). Early and late periodic patterns of even skipped expression are controlled by distinct regulatory elements that respond to different spatial cues. Cell 57 (3), 413–422. doi:10.1016/0092-8674(89)90916-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Green, J. B. A., and Sharpe, J. (2015). Positional information and reaction-diffusion: two big ideas in developmental biology combine. Development 142 (7), 1203–1211. doi:10.1242/dev.114991

PubMed Abstract | CrossRef Full Text | Google Scholar

Green, M. R. (2005). Eukaryotic transcription activation: right on target. Mol. Cell 18 (4), 399–402. doi:10.1016/j.molcel.2005.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamamoto, K., and Fukaya, T. (2022). Molecular architecture of enhancer-promoter interaction. Curr. Opin. Cell Biol. 74, 62–70. doi:10.1016/j.ceb.2022.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanna-Rose, W., Licht, J. D., and Hansen, U. (1997). Two evolutionarily conserved repression domains in the Drosophila Kruppel protein differ in activator specificity. Mol. Cell Biol. 17 (8), 4820–4829. doi:10.1128/mcb.17.8.4820

PubMed Abstract | CrossRef Full Text | Google Scholar

Ilsley, G. R., Fisher, J., Apweiler, R., De Pace, A. H., and Luscombe, N. M. (2013). Cellular resolution models for even skipped regulation in the entire Drosophila embryo. eLife 2, e00522. doi:10.7554/eLife.00522

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishihara, S., Fujimoto, K., and Shibata, T. (2005). Cross talking of network motifs in gene regulation that generates temporal pulses and spatial stripes. Genes cells. 10 (11), 1025–1038. doi:10.1111/j.1365-2443.2005.00897.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaeger, J., and Monk, N. (2014). Bioattractors: dynamical systems theory and the evolution of regulatory processes. J. Physiol. (Lond). 592 (11), 2267–2281. doi:10.1113/jphysiol.2014.272385

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaeger, J., Surkova, S., Blagov, M., Janssens, H., Kosman, D., Kozlov, K. N., et al. (2004). Dynamic control of positional information in the early Drosophila embryo. Nature 430 (6997), 368–371. doi:10.1038/nature02678

PubMed Abstract | CrossRef Full Text | Google Scholar

Jutras-Dubé, L., El-Sherif, E., and François, P. (2020). Geometric models for robust encoding of dynamical information into embryonic patterns. eLife 9, e55778. doi:10.7554/eLife.55778

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalir, S., Mangan, S., and Alon, U. (2005). A coherent feed-forward loop with a SUM input function prolongs flagella expression in Escherichia coli. Mol. Syst. Biol. 1, 2005.0006. doi:10.1038/msb4100010

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, Y. J., Rhee, K., Liu, J., Jeammet, S., Turner, M. A., Small, S. J., et al. (2022). Predictive modeling reveals that higher-order cooperativity drives transcriptional repression in a synthetic developmental enhancer. eLife 11, 11. doi:10.7554/elife.73395

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohwi, M., and Doe, C. Q. (2013). Temporal fate specification and neural progenitor competence during development. Nat. Rev. Neurosci. 14 (12), 823–838. doi:10.1038/nrn3618

PubMed Abstract | CrossRef Full Text | Google Scholar

Koromila, T., Gao, F., Iwasaki, Y., He, P., Pachter, L., Gergen, J. P., et al. (2020). Odd-paired is a pioneer-like factor that coordinates with Zelda to control gene expression in embryos. eLife 9, e59610. doi:10.7554/eLife.59610

PubMed Abstract | CrossRef Full Text | Google Scholar

Lammers, N. C., Kim, Y. J., Zhao, J., and Garcia, H. G. (2020). A matter of time: using dynamics and theory to uncover mechanisms of transcriptional bursting. Curr. Opin. Cell Biol. 67, 147–157. doi:10.1016/j.ceb.2020.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Landge, A. N., Jordan, B. M., Diego, X., and Müller, P. (2020). Pattern formation mechanisms of self-organizing reaction-diffusion systems. Dev. Biol. 460 (1), 2–11. doi:10.1016/j.ydbio.2019.10.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Levine, M., and Davidson, E. H. (2005). Gene regulatory networks for development. Proc. Natl. Acad. Sci. U. S. A. 102 (14), 4936–4942. doi:10.1073/pnas.0408031102

PubMed Abstract | CrossRef Full Text | Google Scholar

Longabaugh, W. J. R., Davidson, E. H., and Bolouri, H. (2005). Computational representation of developmental genetic regulatory networks. Dev. Biol. 283 (1), 1–16. doi:10.1016/j.ydbio.2005.04.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Lord, N. D., Carte, A. N., Abitua, P. B., and Schier, A. F. (2021). The pattern of nodal morphogen signaling is shaped by co-receptor expression. eLife 10, 10. doi:10.7554/elife.54894

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, J. A., El-Sherif, E., and Brown, S. J. (2012). Comparisons of the embryonic development of Drosophila, Nasonia, and Tribolium. Wiley Interdiscip. Rev. Dev. Biol. 1 (1), 16–39. doi:10.1002/wdev.3

PubMed Abstract | CrossRef Full Text | Google Scholar

Mangan, S., and Alon, U. (2003). Structure and function of the feed-forward loop network motif. Proc. Natl. Acad. Sci. U. S. A. 100 (21), 11980–11985. doi:10.1073/pnas.2133841100

PubMed Abstract | CrossRef Full Text | Google Scholar

Manu, , Surkova, S., Spirov, A. V., Gursky, V. V., Janssens, H., Kim, A.-R., et al. (2009). Canalization of gene expression and domain shifts in the Drosophila blastoderm by dynamical attractors. PLoS Comput. Biol. 5 (3), e1000303. doi:10.1371/journal.pcbi.1000303

PubMed Abstract | CrossRef Full Text | Google Scholar

Mau, C., Rudolf, H., Strobl, F., Schmid, B., Regensburger, T., Palmisano, R., et al. (2023). How enhancers regulate wavelike gene expression patterns. eLife 12, 12. doi:10.7554/elife.84969

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayo, A. E., Setty, Y., Shavit, S., Zaslaver, A., and Alon, U. (2006). Plasticity of the cis-regulatory input function of a gene. PLoS Biol. 4 (4), e45. doi:10.1371/journal.pbio.0040045

PubMed Abstract | CrossRef Full Text | Google Scholar

McQueen, C., and Towers, M. (2020). Establishing the pattern of the vertebrate limb. Development 147 (17), dev177956. doi:10.1242/dev.177956

PubMed Abstract | CrossRef Full Text | Google Scholar

Meinhardt, H. (1991). “Models of biological pattern formation and their application to the early development of drosophila,” in Complexity, chaos, and biological evolution. Editors E. Mosekilde, and L. Mosekilde (New York, NY: Springer US), 303–322.

CrossRef Full Text | Google Scholar

Moussian, B., and Roth, S. (2005). Dorsoventral axis formation in the Drosophila embryo--shaping and transducing a morphogen gradient. Curr. Biol. 15 (21), R887–R899. doi:10.1016/j.cub.2005.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, P., and El-Sherif, E. (2020). A systems-level view of pattern formation mechanisms in development. Dev. Biol. 460 (1), 1. doi:10.1016/j.ydbio.2019.10.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Naidu, V. G., Zhang, Y., Lowe, S., Ray, A., Zhu, H., and Li, X. (2020). Temporal progression of Drosophila medulla neuroblasts generates the transcription factor combination to control T1 neuron morphogenesis. Dev. Biol. 464, 35–44. doi:10.1016/j.ydbio.2020.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Nechaev, S., and Adelman, K. (2011). Pol II waiting in the starting gates: regulating the transition from transcription initiation into productive elongation. Biochim. Biophys. Acta 1809 (1), 34–45. doi:10.1016/j.bbagrm.2010.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Negrete, J., and Oates, A. C. (2021). Towards a physical understanding of developmental patterning. Nat. Rev. Genet. 22 (8), 518–531. doi:10.1038/s41576-021-00355-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Oates, A. C., Morelli, L. G., and Ares, S. (2012). Patterning embryos with oscillations: structure, function and dynamics of the vertebrate segmentation clock. Development 139 (4), 625–639. doi:10.1242/dev.063735

PubMed Abstract | CrossRef Full Text | Google Scholar

Palmeirim, I., Henrique, D., Ish-Horowicz, D., and Pourquié, O. (1997). Avian hairy gene expression identifies a molecular clock linked to vertebrate segmentation and somitogenesis. Cell 91 (5), 639–648. doi:10.1016/s0092-8674(00)80451-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Panovska-Griffiths, J., Page, K. M., and Briscoe, J. (2013). A gene regulatory motif that generates oscillatory or multiway switch outputs. J. R. Soc. Interface 10 (79), 20120826. doi:10.1098/rsif.2012.0826

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, J., Estrada, J., Johnson, G., Vincent, B. J., Ricci-Tam, C., Bragdon, M. D., et al. (2019). Dissecting the sharp response of a canonical developmental enhancer reveals multiple sources of cooperativity. eLife 8, e41266. doi:10.7554/eLife.41266

PubMed Abstract | CrossRef Full Text | Google Scholar

Perez-Carrasco, R., Barnes, C. P., Schaerli, Y., Isalan, M., Briscoe, J., and Page, K. M. (2018). Combining a toggle switch and a repressilator within the AC-DC circuit generates distinct dynamical behaviors. Cell Syst. 6 (4), 521–530. doi:10.1016/j.cels.2018.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Pourquié, O. (2003). The segmentation clock: converting embryonic time into spatial pattern. Science 301 (5631), 328–330. doi:10.1126/science.1085887

PubMed Abstract | CrossRef Full Text | Google Scholar

Roensch, K., Tazaki, A., Chara, O., and Tanaka, E. M. (2013). Progressive specification rather than intercalation of segments during limb regeneration. Science 342 (6164), 1375–1379. doi:10.1126/science.1241796

PubMed Abstract | CrossRef Full Text | Google Scholar

Rogers, K. W., and Schier, A. F. (2011). Morphogen gradients: from generation to interpretation. Annu. Rev. Cell Dev. Biol. 27, 377–407. doi:10.1146/annurev-cellbio-092910-154148

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohde, L. A., Bercowsky-Rama, A., Valentin, G., Naganathan, S. R., Desai, R. A., Strnad, P., et al. (2024). Cell-autonomous timing drives the vertebrate segmentation clock’s wave pattern. eLife 13. doi:10.7554/eLife.93764

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, S. (2011). Mathematics and biology: a Kantian view on the history of pattern formation theory. Dev. Genes Evol. 221 (5–6), 255–279. doi:10.1007/s00427-011-0378-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudolf, H., Zellner, C., and El-Sherif, E. (2020). Speeding up anterior-posterior patterning of insects by differential initialization of the gap gene cascade. Dev. Biol. 460 (1), 20–31. doi:10.1016/j.ydbio.2019.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sabari, B. R., Dall’Agnese, A., Boija, A., Klein, I. A., Coffey, E. L., Shrinivas, K., et al. (2018). Coactivator condensation at super-enhancers links phase separation and gene control. Science 361 (6400), eaar3958. doi:10.1126/science.aar3958

PubMed Abstract | CrossRef Full Text | Google Scholar

Samee, M. A. H., Lydiard-Martin, T., Biette, K. M., Vincent, B. J., Bragdon, M. D., Eckenrode, K. B., et al. (2017). Quantitative Measurement and Thermodynamic modeling of Fused enhancers support a two-Tiered mechanism for interpreting regulatory DNA. Cell Rep. 21 (1), 236–245. doi:10.1016/j.celrep.2017.09.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarrazin, A. F., Peel, A. D., and Averof, M. (2012). A segmentation clock with two-segment periodicity in insects. Science 336 (6079), 338–341. doi:10.1126/science.1218256

PubMed Abstract | CrossRef Full Text | Google Scholar

Scholes, C., DePace, A. H., and Sánchez, Á. (2017). Combinatorial gene regulation through kinetic control of the transcription cycle. Cell Syst. 4 (1), 97–108. doi:10.1016/j.cels.2016.11.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Schroeder, M. D., Greer, C., and Gaul, U. (2011). How to make stripes: deciphering the transition from non-periodic to periodic patterns in Drosophila segmentation. Development 138 (14), 3067–3078. doi:10.1242/dev.062141

PubMed Abstract | CrossRef Full Text | Google Scholar

Setty, Y., Mayo, A. E., Surette, M. G., and Alon, U. (2003). Detailed map of a cis-regulatory input function. Proc. Natl. Acad. Sci. U. S. A. 100 (13), 7702–7707. doi:10.1073/pnas.1230759100

PubMed Abstract | CrossRef Full Text | Google Scholar

Shea, M. A., and Ackers, G. K. (1985). The OR control system of bacteriophage lambda. A physical-chemical model for gene regulation. J. Mol. Biol. 181 (2), 211–230. doi:10.1016/0022-2836(85)90086-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, M. S., and Cohen, B. A. (2012). Thermodynamic state ensemble models of cis-regulation. PLoS Comput. Biol. 8 (3), e1002407. doi:10.1371/journal.pcbi.1002407

PubMed Abstract | CrossRef Full Text | Google Scholar

Simsek, M. F., and Özbudak, E. M. (2022). Patterning principles of morphogen gradients. Open Biol. 12 (10), 220224. doi:10.1098/rsob.220224

PubMed Abstract | CrossRef Full Text | Google Scholar

Soluri, I. V., Zumerling, L. M., Payan Parra, O. A., Clark, E. G., and Blythe, S. A. (2020). Zygotic pioneer factor activity of Odd-paired/Zic is necessary for late function of the Drosophila segmentation network. eLife 9, e53916. doi:10.7554/eLife.53916

PubMed Abstract | CrossRef Full Text | Google Scholar

Spitz, F., and Furlong, E. E. M. (2012). Transcription factors: from enhancer binding to developmental control. Nat. Rev. Genet. 13 (9), 613–626. doi:10.1038/nrg3207

PubMed Abstract | CrossRef Full Text | Google Scholar

Stapornwongkul, K. S., and Vincent, J.-P. (2021). Generation of extracellular morphogen gradients: the case for diffusion. Nat. Rev. Genet. 22 (6), 393–411. doi:10.1038/s41576-021-00342-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tufcea, D. E., and François, P. (2015). Critical timing without a timer for embryonic development. Biophys. J. 109 (8), 1724–1734. doi:10.1016/j.bpj.2015.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Verd, B., Clark, E., Wotton, K. R., Janssens, H., Jiménez-Guri, E., Crombach, A., et al. (2018). A damped oscillator imposes temporal order on posterior gap gene expression in Drosophila. PLoS Biol. 16 (2), e2003174. doi:10.1371/journal.pbio.2003174

PubMed Abstract | CrossRef Full Text | Google Scholar

Vincent, B. J., Estrada, J., and DePace, A. H. (2016). The appeasement of Doug: a synthetic approach to enhancer biology. Integr. Biol. (Camb). 8 (4), 475–484. doi:10.1039/c5ib00321k

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolpert, L. (1969). Positional information and the spatial pattern of cellular differentiation. J. Theor. Biol. 25 (1), 1–47. doi:10.1016/s0022-5193(69)80016-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolpert, L. (1989). Positional information revisited. Development 107 (Suppl. l), 3–12. doi:10.1242/dev.107.Supplement.3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yokoshi, M., and Fukaya, T. (2019). Dynamics of transcriptional enhancers and chromosome topology in gene regulation. Dev. Growth Differ. 61 (5), 343–352. doi:10.1111/dgd.12597

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuh, C. H., Bolouri, H., and Davidson, E. H. (1998). Genomic cis-regulatory logic: experimental and computational analysis of a sea urchin gene. Science 279 (5358), 1896–1902. doi:10.1126/science.279.5358.1896

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Rudolf, H., Healey, L., François, P., Brown, S. J., Klingler, M., et al. (2017). Speed regulation of genetic cascades allows for evolvability in the body plan specification of insects. Proc. Natl. Acad. Sci. U. S. A. 114 (41), E8646-E8655. doi:10.1073/pnas.1702478114

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Y., Qiu, Y., Chen, W., Nie, Q., and Lander, A. D. (2020). Scaling a Dpp morphogen gradient through feedback control of Receptors and Co-receptors. Dev. Cell 53 (6), 724–739. doi:10.1016/j.devcel.2020.05.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuo, P., Stanojević, D., Colgan, J., Han, K., Levine, M., and Manley, J. L. (1991). Activation and repression of transcription by the gap proteins hunchback and Krüppel in cultured Drosophila cells. Genes Dev. 5 (2), 254–264. doi:10.1101/gad.5.2.254

PubMed Abstract | CrossRef Full Text | Google Scholar