1. Computational and Systems Biology
  2. Developmental Biology
Download icon

Modularity, criticality, and evolvability of a developmental gene regulatory network

  1. Berta Verd  Is a corresponding author
  2. Nicholas AM Monk
  3. Johannes Jaeger  Is a corresponding author
  1. The Barcelona Institute of Science and Technology, Spain
  2. Universitat Pompeu Fabra (UPF), Spain
  3. Konrad Lorenz Institute for Evolution and Cognition Research (KLI), Austria
  4. University of Cambridge, United Kingdom
  5. University of Sheffield, United States
  6. Wissenschaftskolleg zu Berlin, Germany
  7. Center for Systems Biology Dresden (CSBD), Germany
  8. Complexity Science Hub (CSH), Austria
  9. Centre de Recherches Interdisciplinaires (CRI), France
Research Article
  • Cited 0
  • Views 2,146
  • Annotations
Cite this article as: eLife 2019;8:e42832 doi: 10.7554/eLife.42832

Abstract

The existence of discrete phenotypic traits suggests that the complex regulatory processes which produce them are functionally modular. These processes are usually represented by networks. Only modular networks can be partitioned into intelligible subcircuits able to evolve relatively independently. Traditionally, functional modularity is approximated by detection of modularity in network structure. However, the correlation between structure and function is loose. Many regulatory networks exhibit modular behaviour without structural modularity. Here we partition an experimentally tractable regulatory network—the gap gene system of dipteran insects—using an alternative approach. We show that this system, although not structurally modular, is composed of dynamical modules driving different aspects of whole-network behaviour. All these subcircuits share the same regulatory structure, but differ in components and sensitivity to regulatory interactions. Some subcircuits are in a state of criticality, while others are not, which explains the observed differential evolvability of the various expression features in the system.

https://doi.org/10.7554/eLife.42832.001

Introduction

Systems biology aims to understand the function and evolution of complex regulatory networks. This requires some sort of hierarchical decomposition of these networks into manageable and intelligible subsystems, whose properties and behaviour can be analysed and understood in relative isolation (Simon, 1962; Riedl, 1975; Lewontin, 1978; Bonner, 1988; Raff, 1996; West-Eberhard, 2003; Schlosser and Wagner, 2004; Callebaut et al., 2005). If each subsystem possesses a clearly delimited and discernible function, the network can be subdivided into functional modules (Raff, 1996; von Dassow and Munro, 1999; Hartwell et al., 1999; Wagner et al., 2007; Mireles and Conrad, 2018). In the Introduction of our paper, we provide a careful argument showing that the most common approach to identify functional modules has severe limitations, and propose an alternative method, which we then use to dissect and analyse a specific pattern-forming network, the gap gene system of the vinegar fly, Drosophila melanogaster.

The most common strategy to identify functional modules is to partition the graph representing a network into simple motifs (Shen-Orr et al., 2002; Alon, 2007) or subcircuits (also called subnetworks or communities; (Girvan and Newman, 2002; Oliveri and Davidson, 2004; Babu et al., 2004; Levine and Davidson, 2005; Newman, 2006; Davidson and Erwin, 2006; Oliveri and Davidson, 2007; Erwin and Davidson, 2009; Davidson, 2010). Network motifs are small subgraphs that are identified through their statistical enrichment (Alon, 2007; Alon, 2006), while subcircuits are characterised by a high connection density among their component nodes contrasting with sparse connections to the outside (Girvan and Newman, 2002; Radicchi et al., 2004; Newman, 2006; Wagner et al., 2007; Fortunato, 2010). In both cases, subsystems are defined in terms of the regulatory structure or network topology: they are structural modules. This approach presupposes a strong connection between functional and structural modularity (see, for example, Lim et al., 2013).

Strictly interpreted, structural modules are mutually exclusive: they are disjoint subgraphs of a complex regulatory network that do not share nodes between each other (Girvan and Newman, 2002; Radicchi et al., 2004; Palla et al., 2005). And yet, such modules can never be fully isolated: their context within the larger network influences behaviour and function. The structural approach therefore relies on the assumption that context-dependence is weak, and structural modularity is generally pronounced enough, to preserve the salient properties and behaviour of a motif or subcircuit in its native network context.

Structural modularity is widely regarded as a necessary condition for the evolvability of complex networks. ‘Evolvability,’ in the general sense of the term, is defined as the ability to evolve (Dawkins, 1989; Wagner and Altenberg, 1996; Hendrikse et al., 2007; Pigliucci, 2008). More specifically, evolvability refers to the capacity of an evolving system to generate or facilitate adaptive change (Wagner and Altenberg, 1996; Pigliucci, 2008). Structural modularity can boost this capacity in several ways. Entire modules can be co-opted into new pathways during evolution, generating innovative change (Raff, 1996; von Dassow and Munro, 1999; True and Carroll, 2002; Davidson and Erwin, 2006; Erwin and Davidson, 2009; Monteiro and Podlaha, 2009; Wagner, 2011). Furthermore, each module can vary relatively independently, and it has been argued that this accounts for the individuality, origin, and homology of morphological characters as well as their trait-specific variational properties (Wagner and Altenberg, 1996; Wagner et al., 2007; Wagner, 2014). Finally, structural modularity allows for a fine-tuned response to specific selective pressures by minimizing off-target pleiotropic effects (Wagner and Altenberg, 1996; Pavlicev et al., 2008; Wagner and Zhang, 2011).

The identification and analysis of structural modules has been very successful in many cases. For example, it has been used to understand the regulatory principles of segment determination in Drosophila (von Dassow et al., 2000; Ingolia, 2004), the origin and evolution of butterfly wing spots (Carroll et al., 1994; Brakefield et al., 1996; Keys et al., 1999; Beldade et al., 2002; Monteiro et al., 2003; Monteiro et al., 2006) and beetle horns (Moczek, 2006), and the mechanism and evolution of larval skeleton formation in sea urchins and sea stars (Hinman et al., 2003; Hinman and Davidson, 2007; Oliveri et al., 2008; Gao and Davidson, 2008). Other examples abound in the literature (see Raff, 1996; Schlosser and Wagner, 2004; Callebaut et al., 2005; Peter and Davidson, 2015 for comprehensive reviews).

In spite of its usefulness, structural modularity has a number of serious limitations. Some modelling studies suggest that it is not necessary for evolvability (see, for example, Crombach and Hogeweg, 2008). Furthermore, it is notoriously difficult to identify structural modules and delimit their boundaries with any precision. One reason for this may be that the definition of (sub)system boundaries is fundamentally context- and problem-dependent (see, for example, Chu et al., 2003; Chu, 2011). More to the point, even the simplest subcircuits tend to exhibit a rich dynamic repertoire comprising a range of different behaviours depending on context (boundary conditions), quantitative strength of parameter values (determining genetic interactions as well as production and decay rates), and the specific form of the regulation-expression functions used to integrate multiple regulatory inputs (Mangan and Alon, 2003; Wall et al., 2005; Ingram et al., 2006; Siegal et al., 2007; Payne and Wagner, 2015; Ahnert and Fink, 2016; Perez-Carrasco et al., 2018; Page and Perez-Carrasco, 2018). Because of this, it is usually not possible to single out subgraphs exhibiting specific behaviours and functions that are robustly independent of their native network context. In cases like these, looking for structural modules is not the most fruitful approach to subdivide a complex regulatory network.

A recent simulation-based screen of multifunctional gene regulatory networks nicely illustrates this important point (Jiménez et al., 2017). The authors performed a systematic computational search for network structures able to implement two qualitatively different dynamical behaviours. They then identified the particular subcircuits that were responsible for either of the two behaviours (or functions). What they found is an entire spectrum of structural overlap among functional modules. At one end of the spectrum, ‘hybrid’ networks represent the sum of two completely disjoint structural modules. At the other end, ‘emergent’ networks use exactly the same nodes and connections to implement both dynamical behaviours (note that Salazar-Ciudad et al., 2000; Salazar-Ciudad et al., 2001a; Salazar-Ciudad et al., 2001b) introduce a different definition of ‘emergent’ network structure, indicating a flat hierarchy of regulatory interactions rich in regulatory feedback). Most networks identified by the screen fell somewhere in between these two extremes, that is, they show partial structural overlap between functional modules. This suggests that most functionally modular networks are not modular in the strict structural sense.

The limitations of structural modularity can be further illustrated using the real-world example of the gap gene system. This gene regulatory network is involved in pattern formation and segment determination during the blastoderm stage of early embryogenesis in dipteran insects such as D. melanogaster (see Jaeger, 2011 for review). Its regulatory structure is summarized in Figure 1A. Nodes in this network represent genes encoding transcription factors that activate or repress each other as indicated by their connections. The gap gene network reads and interprets morphogen gradients formed along the major or antero-posterior (A–P) axis of the embryo by the protein products of the maternal co-ordinate genes bicoid (bcd), caudal (cad), and hunchback (hb). This results in broad, overlapping expression domains for the trunk gap genes hb, Krüppel (Kr), knirps (kni), and giant (gt) (Figure 1B, C). Extensive gap-gene cross-regulation, especially during the last cleavage cycle (cycle 14A, C14A) of the blastoderm stage, is essential for the correct dynamic positioning of expression domain boundaries, in particular the dynamic kinematic shifts of posterior gap domains towards the anterior of the embryo (Jaeger et al., 2004b; Jaeger et al., 2004a; Surkova et al., 2008; Manu et al., 2009a; Crombach et al., 2012; Verd et al., 2018). Maternal co-ordinate and gap genes together then regulate pair-rule and segment-polarity genes, which form a molecular pre-pattern at the time of gastrulation that precedes and determines the formation of morphological segments during later stages of development.

The gap gene system: structure, regulation, and expression.

(A) The regulatory structure of the gap gene network. Nodes represent transcription factors encoded by maternal coordinate genes bicoid (bcd, purple) and caudal (cad, cyan), as well as trunk gap genes hunchback (hb, yellow), Krüppel (Kr, green), knirps (kni, red) and giant (gt, blue). Activating interactions are indicated by arrows, repressive interactions by T-bars. Circular arrows indicate auto-activation. Dashed lines represent weak, solid lines stronger regulatory interactions. (B) Spatial representation of gap gene regulation: boxes indicate relative positions of gap domains along the antero-posterior (A–P) axis of the embryo. Anterior is to left, posterior to the right. Background colour indicates predominant activating inputs by maternal gradients. T-bars indicate gap gene cross-repression. Circular arrows represent gap gene auto-activation. (C) Gene expression dynamics of maternal coordinate genes bcd (purple) and cad (cyan), as well as trunk gap genes hb (yellow), Kr (green), kni (red) and gt (blue). Quantified expression patterns are shown as lines (maternal coordinate) and coloured areas (gap genes). Output of the fitted full model is shown as dots. Y-axes represent relative protein concentration in arbitrary units (a.u.); X-axes represent A–P position in %, where 0% is the anterior pole. C12 and C14 refer to cleavage cycle 12 and 14A, respectively. C14A is subdivided into eight time classes T1–8 of equal length. Only the trunk region of the embryo from 35 to 75% is shown. (D) Dynamical regimes driven by the gap gene system. Multistable (switch-like) behaviour in the anterior is indicated by a phase space with trajectories converging to multiple point attractors (represented by circles). Damped oscillations leading to dynamic anterior shifts of gap domains in the posterior are indicated by a colour wheel with trajectories cycling through successions of different gap genes. The bifurcation boundary between the two regimes at 52% A–P position is indicated by a dashed line. See text for details. Figure adapted from Verd et al. (2017); Verd et al. (2018).

https://doi.org/10.7554/eLife.42832.002

Due to its small size and high connection density (see Figure 1A), it is not possible to identify structural clusters in the gap gene system. Previous simulation-based analyses identified a number of mechanisms driving gap gene expression. One of these consists of two double-negative (positive) feedback loops between hb/kni and Kr/gt that are responsible for the basic staggered arrangement of gap domains along the A–P axis. Another one comprises the asymmetric repressive interactions between overlapping gap genes (e. g. kni on Kr and Kr on hb) that are driving dynamic anterior shifts of gap domains over time (Figure 1A, B) (Jaeger et al., 2004a; Jaeger et al., 2004b; Jaeger, 2011; Crombach et al., 2012). Note that these mechanisms are not structurally modular, as they involve different interactions among the same set of genes. Moreover, attempts to explain gap gene expression dynamics with network motifs have consistently failed to correctly account for the integrated behaviour of the whole system (Zinzen and Papatsenko, 2007; Ishihara and Shibata, 2008; Papatsenko, 2009).

All these difficulties arise because the gap gene network is a strongly emergent network (sensu Jiménez et al., 2017). If there are functional modules in this network, they show a high degree of structural overlap. For such emergent networks, analyses based on structural modularity fail, because this approach is only valid for systems of markedly hybrid character. As was argued in Jiménez et al. (2017), such systems only constitute a tiny subset of all possible naturally occurring regulatory networks.

An alternative way to identify network modules that does not rely directly on their structure is to define them in terms of their activity (Kauffman, 1993; Hartwell et al., 1999; von Dassow and Munro, 1999; Newman and Bhat, 2008; Newman and Bhat, 2009; Alexander et al., 2009; Benítez and Alvarez-Buylla, 2010). Such dynamical modules or dynamical subsystems consist of a group of connected network nodes that implement a particular behaviour (Irons and Monk, 2007; Benítez and Alvarez-Buylla, 2010). Their spatio-temporal pattern of activity specifies a certain type of dynamics—such as bistable or oscillatory—termed a ‘dynamical regime’ (Verd et al., 2017; Verd et al., 2018). Different dynamical regimes are distinguished by the composition of the underlying phase portraits (Strogatz, 2015; Perez-Carrasco et al., 2018). For instance, they can be generated by different subsets of a system’s attractors (with associated basins) (Irons and Monk, 2007). With this type of classification, nodes may be shared between overlapping modules, simultaneously driving different dynamics in the context of different subsystems. No disjoint modular network structure is required or even expected (Jiménez et al., 2017). Still, each module’s behaviour is relatively independent and clearly distinguishable from that of other dynamical subsystems (Irons and Monk, 2007).

Note that dynamical modules are not the same as co-expression modules (sometimes also called regulatory modules; see, for example, Eisen et al., 1998; Bar-Joseph et al., 2003; Segal et al., 2003; Stuart et al., 2003; Kim et al., 2009; Koch et al., 2017). The latter are defined by the correlated or anti-correlated expression of their components. In contrast, the mutual dependence of component expression patterns in dynamical modules is causal, rather than correlative. It can be very complex and obscured by the non-linearity of regulatory interactions. Dynamical modules are not identified by individual expression patterns, but rather by the coherent collective activity of the module as a whole (Kauffman, 1993; Irons and Monk, 2007; Alexander et al., 2009).

Just like their structural counterparts, dynamical modules also influence evolvability. However, the way they achieve this is fundamentally different. Structural subcircuits can vary independently since they are only loosely interconnected. Dynamical modules, in contrast, drive distinct behaviours that exhibit different levels of sensitivity to changes in system parameters (Alexander et al., 2009). In particular, some dynamical modules may be in a state of criticality, close to a bifurcation point beyond which their dynamics may change drastically and abruptly; others may be structurally stable, that is, not critical and far from any bifurcation, and therefore insensitive or robust to changes in parameters (Thom, 1976; Kauffman, 1993; Jaeger et al., 2012; Jaeger and Monk, 2014; Jaeger and Sharpe, 2014; Strogatz, 2015). Recall that specific network components may be involved in governing more than one dynamic behaviour. Mutations affecting these components will tend to have a strong effect on those behaviours that are critical, leaving others unaltered because they are structurally stable. This type of dynamic modularity implies that the network is more likely to evolve in certain directions than others. While mutations may be random, their effects on network dynamics are certainly not. Behaviours driven by modules that are critical (close to bifurcation) will be more labile than those driven by modules that are structurally stable. This amounts to functional modularity, even though the overall structure of the network is not modular.

Dynamical modularity transcends the inherent limitations of structural approaches and allows us to gain insights into heavily emergent regulatory networks. However, it is far from trivial to apply this framework to the experimental study of specific evolving developmental processes. Irons and Monk have developed an algorithmic method to identify dynamical modules in Boolean network models (Irons and Monk, 2007). Unfortunately, this method is difficult to generalize and adapt to continuous mathematical frameworks (for instance, models formulated as systems of ordinary differential equations) used to study morphogen-driven pattern-forming networks such as the gap gene system.

For this reason, we adopt a pragmatic approach to identify dynamical modules in the gap gene network of D. melanogaster. Our approach is based on the observation that only subsets of gap genes are expressed and exert their regulatory influence in any one region of the embryo. This allowed us to identify four localized subsystems, each containing three trunk gap genes that are active in different but overlapping regions of the embryo. Surprisingly, all four subsystems share the same regulatory structure, which identifies them as AC/DC subcircuits (Panovska-Griffiths et al., 2013; Perez-Carrasco et al., 2018). AC/DC subcircuits represent one of the simplest known genetic systems able to produce both switch-like (multistable) and oscillatory behaviour. We show that these AC/DC modules drive distinct dynamical regimes: static domain boundaries in the anterior, versus anteriorly shifting gap domains in the posterior of the embryo (Figure 1C, D; Verd et al., 2017; Verd et al., 2018). The boundary between these two dynamical regimes is positioned by an AC/DC circuit in a state of criticality. This makes the position of this boundary especially sensitive to changes in the strength of regulatory interactions, and we argue that this has shaped the evolvability of the gap gene system within the Dipteran order (flies, midges, and mosquitoes).

Results and discussion

Modularity

The gap gene system is composed of three dynamical modules

We have taken a pragmatic approach to identifying dynamical modules in the gap gene network (Figure 1A, BJaeger, 2011). The analysis asks which of the four trunk gap genes (hb, Kr, kni, and gt) are required—or, more accurately, which ones are not—to drive correct expression dynamics in nuclei at different positions along the antero-posterior (A–P) axis during cleavage cycle 14A (C14A). It is based on a detailed dynamical model of the gap gene network (the ‘full model’), which has been fit to quantitative spatio-temporal data of trunk gap gene expression (see Figure 1C) (Verd et al., 2017; Verd et al., 2018). This model implements accurate dynamic mechanisms of gap gene regulation that have been extensively validated against empirical evidence (Jaeger et al., 2004b; Jaeger et al., 2004a; Manu et al., 2009b; Manu et al., 2009a; Ashyraliyev et al., 2009; Jaeger, 2011; Crombach et al., 2012; Verd et al., 2017; Verd et al., 2018). We consider the sensitivity of the network to a gap gene negligible, if the node representing that gene in the network can be removed from the model at the onset of C14A without significant consequences to the resulting expression dynamics (see ‘Node sensitivity analysis’ in ‘Materials and methods’). If sensitivity is negligible in a given nucleus, we conclude that the gene is not required to drive gap gene expression in that nucleus during C14A.

The results of our analysis are shown in Figure 2A. They reveal three regions that are insensitive to specific gap genes: (1 in the region between 35 and 47% A–P position, developmental trajectories are insensitive to kni (red background); (2) between 49 and 59% A–P position, they are insensitive to gt (blue), and (3) between 61 and 75% A–P position, they are insensitive to hb (yellow). Therefore, the gap gene regulatory network can be reduced from four to three trunk gap genes in each of these three regions (Figure 2A and B). Regional boundaries reflect the position of expression boundaries, but differ from those in that they remain constant, while expression patterns change over time during C14A (Surkova et al., 2008; Jaeger, 2011).

Identification of dynamical modules in the gap gene system.

(A) Node sensitivity analysis. The plot shows sensitivity of model output to the removal of hb (yellow), Kr (green), kni (red) and gt (blue). Y-axis represents node sensitivity as defined in ‘Materials and methods,’ Equation (4); X-axis represents A–P position in %, where 0% is the anterior pole. Regions insensitive to the absence of specific gap genes are indicated by background colour. (B) Dissecting the gap gene network into dynamical modules. Network schemata as in Figure 1B. Subcircuits active in each region identified in (A) are highlighted. (C) AC/DC subcircuits. All subcircuits identified in (A) and (B) share the same regulatory structure, indicated by T-bar connectors. Numbers indicate strength of interactions (in arbitrary units). Maternal inputs and auto-activation are omitted for clarify. Note that there is a fourth AC/DC subcircuit posterior of the region included in the present analysis (not shown, see ‘Materials and methods’). See text for details.

https://doi.org/10.7554/eLife.42832.003

Next, we single out a minimum set of gap gene interactions that are still able to drive correct expression dynamics in each of these three regions. Surprisingly, the structure of all three resulting subcircuits is qualitatively the same, even though each involves a different (overlapping set of gap genes (Figure 2B), and therefore distinct strengths of regulatory interactions (Figure 2C). This structure combines a negative feedback loop between all three genes of a subcircuit with a double-negative (positive) feedback loop between two of them (Figure 2B, C). Such a combination between positive and negative feedback loops is called the ‘AC/DC circuit,’ first described in the context of dorso-ventral patterning in the vertebrate neural tube (Balaskas et al., 2012; Panovska-Griffiths et al., 2013; Perez-Carrasco et al., 2018). Positive feedback loops occur between non-overlapping gap genes Kr and Gt, as well as Hb and Kni; the interactions involved are much stronger than the negative repressive interactions between gap genes with overlapping expression domains (Figure 2C). Previous work has shown that strong positive feedback is required for the basic staggered arrangement of gap domains (‘alternating cushions’), while weaker (and hence slower) negative feedback drives anterior shifts in gap domain position over time in the posterior trunk region of the embryo (Jaeger et al., 2004b; Jaeger et al., 2004a; Perkins et al., 2006; Ashyraliyev et al., 2009; Jaeger, 2011; Crombach et al., 2012; Verd et al., 2018).

The dynamics of AC/DC subcircuits faithfully reproduce the dynamics from the full model

The next step is to establish whether AC/DC subcircuits are sufficient for patterning in each of the three embryonic regions identified in Figure 2A. In order to qualify as a true dynamical module sensu Irons and Monk (2007), each AC/DC subcircuit must recover the expression dynamics as well as the underlying dynamical regime of the full model in the region where it is active. Anterior to the bifurcation boundary at 52% A–P position, this means static gap domain borders governed by multi-stability; posterior to 52% A–P position, this means kinematically shifting domain boundaries governed by a monostable dynamical regime that drives a stereotypical temporal succession of gap gene expression in each nucleus (Figure 1B, C and Figure 3A) (Verd et al., 2017; Verd et al., 2018).

Figure 3 with 2 supplements see all
AC/DC subcircuits faithfully reproduce the dynamics of the full model.

(A) Spatio-temporal dynamics of gap gene expression in the trunk region of the embryo. Coloured areas show the position of gap domains for hb (yellow), Kr (green), kni (red), and gt (blue). Y-axis represents time (flowing downwards) during cleavage cycle C14A. Time classes T1–eight as defined in Figure 1. X-axis represents % A–P position, where 0% is the anterior pole. Regions of influence for each AC/DC subcircuit is indicated by grey lines. Black dotted line indicates the bifucation at 52% position separating static from shifting gap domain boundaries. (B) Comparative dynamical analysis of AC/DC subcircuits and the full model. Instantaneous phase portraits of AC/DC1 (nucleus at 37%), AC/DC2 (nuclei at 49 and 59%), and AC/DC3 (nucleus at 63% A–P position) are shown. Point attractors are shown as spheres, spiral sinks as cylinders, saddles as cubes. Colour code indicates time class, from T1 (dark red) to T8 (yellow). Trajectories from simulations of AC/DC subcircuits are shown in turquoise, trajectories from simulations of the full model in black. Axes represent concentrations of gap proteins (in arbitrary units) as indicated. See ‘Materials and methods’ for model definition and details on phase space analysis.

https://doi.org/10.7554/eLife.42832.004

All three AC/DC subcircuits—AC/DC1 in the anterior, AC/DC2 in the middle, and AC/DC3 in the posterior—reproduce the expression dynamics of the full model with reasonable accuracy in their respective regions of influence throughout cleavage cycle 14A, which we subdivide into eight time classes T1–8 of equal length for the purpose of our analysis (Figure 3—figure supplement 1) (Jaeger et al., 2004b; Jaeger et al., 2004a; Surkova et al., 2008). The only major expression defect involves the anterior boundary of the posterior gt domain, which fails to shift in the AC/DC3 model, starting to deviate from the full model around time class T3 (Figure 3—figure supplement 1). Other defects affect aspects of domain shape or levels of expression (Figure 3—figure supplement 1), but not the dynamic positioning of boundary interfaces, on which our analysis is focused.

The range of dynamical regimes that can be implemented by a circuit is constrained by its regulatory structure (see ‘Introduction’. Which particular dynamical regime is realized depends on the strength of regulatory interactions as well as initial and boundary conditions provided by maternal factors. Since all three AC/DC subcircuits exhibit different interaction strengths (Figure 2C) and receive different maternal inputs, they can drive different expression dynamics (see Appendix 1 for a systematic mathematical analysis).

We determine the dynamical regime of each AC/DC subcircuit by calculating and analysing their associated phase portraits. To achieve this, we simulate each AC/DC subcircuit in isolation, including realistic time-variable inputs from maternal gradients (Verd et al., 2017) and auto-activation. It has been shown previously that auto-activation does not significantly contribute to patterning in the gap gene system, but is required to regulate accurate levels of gap gene expression (Perkins et al., 2006). The time-dependent boundary conditions render these systems time-variant or non-autonomous, and cause their phase portraits to change their geometry over time (Verd et al., 2017). Time-variance or non-autonomy means that not only the variables, but also the parameters or boundary conditions change over time, and therefore implicitly the rules governing the regulatory system themselves. Previously, we have developed a successful method for the classification of transient trajectories in time-variant, non-autonomous dynamical systems (Verd et al., 2014), and have used it to study the dynamics of gap gene expression in the full model (Verd et al., 2017; Verd et al., 2018). Here, we redeploy this technique to study each AC/DC subcircuit independently (see ‘Materials and methods’ for details).

AC/DC1 consists of nodes representing hb, Kr, and gt (Figure 2C, left panel. Its region of influence lies between 35 and 47% A–P position (Figure 2A). Phase portraits in this region exhibit multistability with two or more point attractors (shown for the nucleus at 37% A–P position in Figure 3B, left; additional nuclei are shown in Figure 3—figure supplement 2A). These trajectories are very similar to those simulated with the full model, suggesting that flow direction and magnitude is conserved between both models (Figure 3B and Figure 3—figure supplements 2A, A’). In each case, trajectories are shaped by the pursuit of a moving point attractor (Verd et al., 2014), located at equivalent positions in phase space in AC/DC1 and the full model. The only difference between the two models is an additional attracting steady state in the full system, situated at high Kr levels. This attractor is positioned too far from the trajectories to influence their shape.

The nodes in AC/DC2 represent hb, Kr, and kni (Figure 2C, central panel), and it is active between 49 and 59% A–P position (Figure 2A). This region straddles the bifurcation that occurs at 52% A–P position in the full model, designating the division between static anterior and shifting posterior gap domains (Manu et al., 2009a; Gursky et al., 2011; Verd et al., 2017). To be considered a dynamical module throughout this region, the AC/DC2 subcircuit must recover the two distinct dynamical regimes on either side of the bifurcation boundary. The phase portraits of a nucleus anterior and another one posterior to the bifurcation boundary are shown in Figure 3B (centre). Additional nuclei are shown in Figure 3—figure supplement 2B.

The nucleus at 49% is located just anterior to the bifurcation. The phase portrait of AC/DC2 in this nucleus is bistable for all time points except the last one (at T8, when it becomes monostable and only an attractor close to the origin persists (Figure 3B). This is extremely similar to what happens in the full model, where by T8 the steady states at high Kr values have disappeared, and only steady states close to the origin remain. This transition to monostability occurs too late to significantly alter the course of the trajectory. Instead, the geometry of the trajectory in this nucleus—almost identical in AC/DC2 and the full model—first converges towards a saddle during early C14A, and later directly towards a moving attractor (Figure 3B). In the full model, a geometric capture occurs at stages prior to C14A (Verd et al., 2014; Verd et al., 2017), outside the time range of the AC/DC models (Verd et al., 2017).

Steady states in the phase portrait of AC/DC2 at 49% A–P position are restricted to the Hb-Kr plane. They can be mapped onto a subset of steady states present in the full model (Figure 3—figure supplements 2B, B’). Additional steady states in the full model are located on the Hb-Gt plane, at positions which are very similar to those present in AC/DC1 in more anterior nuclei (Figure 3—figure supplement 2A). This suggests that phase portraits of the full model are composites of AC/DC1 and AC/DC2 in the region anterior to the bifurcation boundary at 52% A–P position.

Posterior to the bifurcation at 52% A–P position, the full model exhibits trajectories that curve towards a spiral sink attractor in a monostable phase portrait (Verd et al., 2018). The trajectories of AC/DC2 in this region are very similar (Figure 3B for the nucleus at 59% A–P position, see Figure 3—figure supplements 2B, B’ for trajectories in other nuclei). Both exhibit a spiralling geometry which is confined to the Kr-Kni plane, with some minor deviations in peak concentration at late stages.

Interestingly, however, the underlying topology of phase space is not equivalent in the AC/DC2 subcircuit and the full model. While both models exhibit monostability in this region, and trajectories in both are shaped by the pursuit of a moving attractor, there are no spiral sinks in AC/DC2 (Figure 3B, and Figure 3—figure supplement 2B). Spiralling trajectories in this subcircuit are shaped mainly by the movement of a conventional point attractor, with spiral sinks appearing only late, at time points T7 and T8. Attractor movement is much more pronounced than in the full model (Figure 3B, and Figure 3—figure supplements 2B, B’). Furthermore, the steady states of AC/DC2 are located along the Kni axis with early attractors at high concentration values of Kni, while the full model shows steady states arranged along the Gt axis with early attractors at high concentration values of Gt (Figure 3—figure supplement 2B’). This is similar to steady states in AC/DC3 in the nucleus at 63% A–P positions (Figure 3B, right panel). At later stages, when the influence of attractor position on the geometry of the trajectory becomes more pronounced, the positions of steady states converge between the two models: in each case, they are located at low Kni values, along the Kni axis. This illustrates that identical transient behaviour can be caused by different types of moving attractors in (biological) time-variant or non-autonomous dynamical systems.

The nodes in AC/DC3 represent Kr, kni, and gt (Figure 2C, right panel). It is active between 61 and 75% A–P position (Figure 2A). In the nuclei at 61 and 63% A–P position, AC/DC3 faithfully reproduces the dynamics of gap gene expression obtained from the full model (Figure 3B, right, and Figure 3—figure supplement 2C). As in the case of AC/DC2, the system is monostable, and trajectories are shaped by the pursuit of a moving attractor. This attractor is located at high Gt concentrations early on, moving closer to the origin over time with only residual levels of Kni left at late stages (Figure 3B, and Figure 3—figure supplement 2C). Similar to AC/DC2, the type of attractor differs between the full model, where it is a spiral sink at all time points, and AC/DC2, where it is a conventional point attractor during time classes T1–6, only turning into a spiral sink at time points T7 and T8 (Figure 3B, and Figure 3—figure supplements 2C, C’).

Posterior to 63% A–P position, AC/DC3 no longer recapitulates gap gene expression dynamics accurately, because its trajectories fail to switch from the Kr-Kni to the Kni-Gt plane as they do in the full model (Figure 3—figure supplement 1C). It is possible that this switch requires overlap with a fourth AC/DC subcircuit in the posterior sub-terminal region of the embryo, which could not be characterised further since most of its region of influence lies outside the spatial range we can analyse in the full model (see ‘Materials and methods’).

In summary, phase space analysis of AC/DC subcircuits establishes that they are true dynamical modules of the gap gene network in the region between 35% and 63% A–P position (Irons and Monk, 2007). They faithfully recover the geometry of trajectories in the full model, whose phase portrait can be seen as an overlapping composite of those of the subcircuits in this region. This simple picture is complicated by the fact that spiralling trajectories occur posterior to 52% A–P position despite the absence of spiral sinks in the phase portraits of AC/DC models. The lack of spiral sinks is compensated by more pronounced movements of conventional point attractors, suggesting that distinct time-variant or non-autonomous phase space topologies can generate equivalent transient dynamics. We discuss this somewhat degenerate relationship between phase space topology and trajectory shapes further in the ‘Conclusions.’

Criticality and evolvability

AC/DC1 in the anterior and AC/DC3 in the posterior each have a consistent dynamical regime across their respective regions of influence (Figures 13), indicating that these two dynamical modules are structurally stable with respect to inputs from maternal gradients. In contrast, AC/DC2 correctly reproduces the bifurcation observed in the full model at 52% A–P position, which separates static anterior patterning from kinematically shifting domains in the posterior (Figure 3) (Manu et al., 2009a; Gursky et al., 2011; Verd et al., 2017). The presence of a bifurcation implies that this dynamical module is in a state of criticality with respect to the inputs from maternal gradients. In other words, the parameters of the AC/DC2 circuit place it very close to a bifurcation boundary, which it will cross due to the different maternal inputs it receives along the A–P axis of the embryo within its region of influence. Our analysis therefore suggests that the gap gene network of D. melanogaster is critical to maternal inputs only in the middle of the embryo, where AC/DC2 is active between 49 and 59% A–P position, while it is structurally stable outside of this region. This contrasts with an earlier proposition—based on the quantification of cross-correlations between expression patterns and a set of purely theoretical models of gap gene regulation—which suggested that the system shows signs of criticality along the entire antero-posterior axis (Krotov et al., 2014).

Understanding criticality in complex regulatory networks is far from trivial, since bifurcations may depend on more than just one of the system’s parameters (see, for example, Thom, 1976; Scheffer, 2009; Kuznetsov, 2004). This can lead to counter-intuitive effects. For instance, gap gene patterning in the central region of the D. melanogaster blastoderm is robust towards variation in the levels of maternal gradients (Figure 4) (Houchmandzadeh et al., 2002; Gregor et al., 2007a; Manu et al., 2009b; Gursky et al., 2011). This is surprising in light of the fact that AC/DC2 is in a critical state. In contrast, the position of the bifurcation boundary between static and shifting gap domains is labile between different dipteran species (Jaeger et al., 2004b; Surkova et al., 2008; Manu et al., 2009a; García-Solache et al., 2010; Jaeger, 2011; Crombach et al., 2014; Wotton et al., 2015). In the scuttle fly Megaselia abdita (Phoridae), it is located more anteriorly compared to D. melanogaster: the region where gap domain shifts occur includes the Hb-Kr interface at around 40% A–P position (see Figure 5 below) (Wotton et al., 2015). The moth midge Clogmia albipunctata (Psychodidae) shows even more extended and pronounced gap domain shifts (García-Solache et al., 2010; Crombach et al., 2014). In the following subsections, we will provide an analysis that resolves this apparent paradox.

Intra-species robustness of gap gene patterning to perturbations in levels of maternal gradients.

(A) The position of gap gene expression features is relatively robust towards changes in maternal gradient concentrations. This sketch—modified from Figure 1 of Manu et al. (2009b)—shows concentration variation in the maternal Bcd gradient (above) and the zygotic expression domain of Kr (below). Solid line and coloured area indicate averaged expression patterns, dashed lines indicate expression variation. Red arrows show the difference between the range of positional variation in Bcd (red background) and Kr. (B) Phase diagram for AC/DC2 in response to variation in Bcd concentration. Pink dots and error bars show average Bcd concentration with standard deviation between 47 and 61% A–P position. Dashed lines show the maximum range of Bcd profiles in the data. AC/DC2 subcircuits for nuclei indicated by the X-axis, where simulated with Cad concentration fixed to its value at T1, and Bcd concentration fixed to the values given by the Y-axis. Background colour indicates the resulting dynamical regime: the multistable anterior regime is shown in black, the monostable posterior regime in grey. (C) Phase diagram for AC/DC2 in response to variation in Cad concentration. White dots indicate Cad concentrations between 47 and 61% A–P position. AC/DC2 subcircuits for nuclei indicated by the X-axis, where simulated with Bcd concentration fixed to its value at T1, and Cad concentration fixed to the values given by the Y-axis. Background indicates dynamical regime as in (B). See text for details.

https://doi.org/10.7554/eLife.42832.007
Inter-species lability of the bifurcation boundary depends on altered gap-gap interactions.

(A) The position of the bifurcation boundary separating static and shifting gap gene expression domains differs between D. melanogaster (upper panel, 52%) and M. abdita (lower panel, 40% A–P position). The difference between species is highlighted by a grey arrow. Phylogenetic distance between the two species is indicated to the left in ‘million years ago’ (mya). (B) Phase diagram for AC/DC2 in response to altering the strength of hb repression by Kr (plotted against A–P position). (C, D) Phase diagram for AC/DC2 in response to altering both hb repression by Kr (X-axis) and Kr repression by Kni (Y-axis) shown for subcircuits in nuclei at 49 (C) and 59% (D) A–P position. (E) Phase diagram for AC/DC1 in response to altering the strength of hb repression by Kr (plotted against A–P position). (F) Phase diagram for AC/DC1 in response to altering both hb repression by Kr (X-axis) and gt repression by Hb (Y-axis) shown for the subcircuit in the nucleus at 37% A–P position. (AF) Background indicates dynamical regimes as in Figure 4: the multistable anterior regime is shown in black, the monostable posterior regime in grey. In addition, there is a narrow strip of a multistable oscillatory regime in (F) (shown in red). All Phase diagrams were calculated with maternal gradient concentrations fixed to their values at T1. See text for details.

https://doi.org/10.7554/eLife.42832.008

Intra-specific robustness to maternal gradient concentration

We used the AC/DC2 model to assess the effect of changing the maternal Bcd and Cad gradients on the position of the bifurcation boundary (cf. Figure 4). In the region between 47 and 61% A–P position, Bcd spans a range of concentration from 24 to 7 arbitrary units in our data (Figure 4B). For each nucleus in this region, we fixed Cad levels to their values during early C14A (time point T1), while varying Bcd levels in steps of 0.01 from to 30 arbitrary units. We then calculated the phase portraits of AC/DC2 for each combination of maternal input concentrations (see ‘Materials and methods’).

The resulting dynamical regimes are shown in the phase diagram displayed in Figure 4B. There is a threshold response at a Bcd concentration of approximately 15 arbitrary units: above this threshold, the system is bistable; below, it is monostable. Bcd concentrations in nuclei anterior of 52% fall into the multistable regime, nuclei posterior of this position fall into the monostable regime (see dots in Figure 4B), recovering the same bifurcation position as in the full model (Manu et al., 2009a; Gursky et al., 2011; Verd et al., 2017).

Quantitative measurements of the Bcd gradient reveal a standard deviation of about 15% relative concentration (Gregor et al., 2007b; Liu et al., 2013). Plotting this as error bars in Figure 4B, we find that they cross the boundary between dynamical regimes only in nuclei between 51 and 53% A–P position. This very narrow spatial range is in agreement with observed levels of variability in the data and in previous analyses of gap gene circuits (Gregor et al., 2007b; Manu et al., 2009b; Gursky et al., 2011). If we alter the concentration of Bcd by a larger amount, as happens in mutants with varying numbers of bcd copies (Driever and Nüsslein-Volhard, 1988; Houchmandzadeh et al., 2002; Liu et al., 2013), we see a more pronounced displacement of the boundary as indicated by the intersection of the dashed gradient profiles with the bifurcation threshold in Figure 4B.

The situation is very different for the maternal Cad gradient. In the region between 47 and 61% A–P position, Cad spans a range of concentrations from 70 to 110 arbitrary units (Figure 4C). For each nucleus in this region, we fixed Bcd levels to their values during early C14A (time point T1), while varying Cad levels in steps of 0.1 from 60 to 120 arbitrary units. We then calculated the phase portraits of AC/DC2 for each combination of maternal input concentrations (see ‘Materials and methods’).

The resulting phase diagram shows that, in contrast to Bcd, the position of the bifurcation boundary is entirely insensitive to Cad concentration. There is an abrupt vertical divide between dynamical regimes around 52% A–P position (Figure 4C). This implies that varying Cad concentration has no effect on the position of the bifurcation, and therefore the extent of the posterior region that exhibits dynamic gap domain shifts. This complements earlier studies which suggested that Cad serves as a general activator in the posterior, but is not directly involved in setting the position of gap gene expression features (Struhl et al., 1992; Manu et al., 2009a; Verd et al., 2017), or in controlling the rate of gap domain shifts in this region of the embryo (Verd et al., 2018).

Inter-specific lability of the bifurcation boundary

The analysis described in the previous section predicts that increased levels of Bcd should lead to a posterior displacement of the bifurcation boundary between the static anterior and the shifting posterior patterning regime (Figure 4B). Surprisingly, exactly the opposite is observed in M. abdita. In this species, the anterior localization domain of bcd mRNA is expanded compared to D. melanogaster, which presumably leads to an expanded Bcd protein gradient (Stauber et al., 1999; Stauber et al., 2000; Wotton et al., 2015; Crombach et al., 2016). However, despite more Bcd being present, the bifurcation boundary is located further anterior than in D. melanogaster (Figure 5A) (Wotton et al., 2015; Crombach et al., 2016). This is impossible to explain if we only take maternal inputs to the gap gene system into account.

Previous studies have shown that the strength of gap-gap cross-regulatory interactions differs between M. abdita and D. melanogaster (Wotton et al., 2015; Crombach et al., 2016). Based on this, we sought to identify the specific gap-gap interactions that could account for the altered position of the bifurcation in our models. The dissection of the gap gene network into dynamical modules narrows this search to interactions within the AC/DC1 and AC/DC2 subcircuits (Figure 2B, C), since their regions of influence cover the relevant region of the embryo. First, we focus on AC/DC2, because the bifurcation boundary in D. melanogaster lies within its region of influence. We know from previous analyses, that the strong positive feedback between hb and kni is heavily conserved between the two species (Wotton et al., 2015; Crombach et al., 2016). We also know that asymmetric interactions between overlapping gap genes (kni on Kr, and Kr on hb) are involved in regulating domain shifts in both species (Jaeger et al., 2004b; Jaeger et al., 2004a; Ashyraliyev et al., 2009; Crombach et al., 2012; Crombach et al., 2016). Therefore, we concentrate our analysis on these two interactions.

In our D. melanogaster model, repression of hb by Kr is negligible (Figure 2C). In contrast, both genetic evidence (Wotton et al., 2015) and models (Crombach et al., 2016) for M. abdita indicate that there is a significant net repressive effect of Kr on hb in this species. In Figure 5B, we plot the dynamical regimes of AC/DC2 in nuclei between 47 and 61% A–P position, while varying this regulatory parameter by decreasing steps of 0.0005 across a narrow range of values from to −0.0035. If repression remains minimal (closer to zero than −0.0005), we recover the bifurcation at 52% A–P position, as nuclei anterior to this position are multistable, while more posterior nuclei are monostable (Figure 5B). If the strength of repression is further increased, the bifurcation boundary moves anteriorly and vanishes altogether around a repression strength of −0.0025. In contrast, the position of the bifurcation boundary is largely insensitive to the interaction between kni and Kr (Figure 5, D). Taken together, the analysis predicts that AC/DC2 in M. abdita should be structurally stable, as net repression of hb by Kr is increased in this species (Wotton et al., 2015; Crombach et al., 2016), eliminating the bifurcation boundary present in D. melanogaster.

However, the bifurcation boundary in M. abdita is located even further anterior, as shifting boundaries extend to around 40% A–P position (Wotton et al., 2015; Crombach et al., 2016), far into the region between 35 and 45% A–P position covered by AC/DC1. Since the repression of hb by Kr is shared between AC/DC1 and AC/DC2, we asked if increasing its strength would induce a bifurcation in AC/DC1, abolishing its structural stability and rendering it critical. Interestingly, this is not the case, as AC/DC1 remains multistable regardless of repression strength due to the very strong positive feedback between Kr and gt (Figure 5). In this subcircuit, further alterations to regulatory interactions are required to render the circuit monostable. This can be achieved, for example, by simultaneously changing both repressive interactions of Kr on hb and Hb on gt (Figure 5). In this case, the latter interaction is the critical one. Therefore, at least two repressive interactions must be stronger in M. abdita than in D. melanogaster to render AC/DC1 critical in this species.

In summary, our analysis suggests that in D. melanogaster only AC/DC2 is critical, while in M. abdita the bifurcation boundary between static and dynamic patterning regimes falls into the region of influence of AC/DC1. This change in the stability of dynamical modules is caused by changes in the strength of particular gap-gap cross-repressive interactions. Dynamical modules and their criticality—how close they are to a bifurcation boundary—have important consequences for the evolvability of the gap gene network. Our work predicts that small changes in the strength of gap-gap cross-regulatory interactions specifically affect the extent of static versus shifting patterning regimes. In contrast, other aspects of gap gene patterning, such as the alternating-cushions mechanism of positive feedback between Kr and gt as well as hb and kni, are extremely robust towards changes in interaction strengths.

Conclusions

Stuart Kauffman’s notion of adaptive systems at the ‘edge of chaos’ first encapsulated the idea that evolving regulatory networks exhibit modular dynamics and are in a state of criticality (Kauffman, 1993). (Irons and Monk, 2007) later made the idea of dynamical modules precise, and formulated an algorithm to detect them in Boolean network models. Here, in the first part of our analysis, we generalize the notion of a dynamical module to continuous systems and use a pragmatic approach based on sensitivity analysis to identify dynamical modules in the empirically tractable gap gene system of D. melanogaster. The three modules described include distinct sets of regulators, but share a common regulatory network structure. They all correspond to AC/DC subcircuits, able to drive multistable (switch-like) and oscillatory dynamics depending on parameter values and boundary conditions, amongst others (Panovska-Griffiths et al., 2013). In this paper, we show that each circuit is active in a particular region of influence along the antero-posterior axis, where it is able to reproduce the geometry of transient gap gene expression trajectories, and hence the overall expression dynamics, of a full gap gene circuit (Figures 2 and 3).

(Benítez and Alvarez-Buylla, 2010; Benítez and Alvarez-Buylla, 2010) have used a similar approach to study the robustness of pattern formation in the root and leaf epidermis of the mustard cress Arabidopsis thaliana (see also Benítez et al., 2008; Benítez et al., 2011). These authors define ‘dynamic modules’ in a manner quite similar to Irons and Monk (2007). Their modules are composed of overlapping sets of components and interactions that generate a specific dynamic behaviour, or set of attractors, in the context of a given tissue (leaf or root). In contrast to our analysis, all of these modules are able to faithfully reproduce the behaviour of the full model. They are therefore heavily redundant and their linkage is responsible for making the patterning outcome of the whole network robust, rather than generating a diversified pattern in different regions of a tissue. This indicates that dynamical modularity can not only be used to analyse the multifunctionality of gene regulatory networks, but also the robustness or stability of their behaviour.

Dynamical modules as defined and used here should not be confused with dynamic patterning modules (Newman and Bhat, 2008; Newman and Bhat, 2009; Hernández-Hernández et al., 2012), which require an interaction of conserved molecular regulatory networks with generic biophysical properties of aggregated cells. Dynamic patterning modules can be considered a specific sub-type of dynamical modules, which exhibit characteristic and robust behaviour based on the self-organising reciprocal interaction of genetic and tissue-level regulation. While dynamic patterning modules are predominantly used to explain the emergence of biological form during the early evolution of multicellular organisms, dynamical modules are intended for the decomposition and functional characterization of general cellular or developmental processes and their evolutionary potential (see below).

Our dynamical modules should also not be confused with the recent postulation of ‘static’ and ‘dynamic’ modules in the segment determination network of the flour beetle Tribolium castaneum (Zhu et al., 2017). ‘Static’ and ‘dynamic’ in that case refer to slow versus rapid time scales of expression dynamics driven by distinct classes of subcircuits, which are all defined in terms of their disjoint structural components. In contrast, our results indicate that structural modularity is not essential for the evolution of insect segmentation. The search for structural modules may be in vain, as it is in the case of the gap gene network of D. melanogaster where only dynamical modules are present in the system. In this sense, dynamical modules provide a powerful complementary alternative to identifying structural modules.

Concerning the dynamical analysis of the AC/DC subcircuits, it is interesting to note that the geometry of transient trajectories generated by different time-variant or non-autonomous models can be equivalent despite underlying discrepancies in features of phase space (see ‘Results and discussion’. In particular, posterior subcircuits AC/DC2 and AC/DC3 exhibit spiralling trajectories in the absence of spiral sink attractors (Figure 3, and Figure 3B, B’, C, C’) (Verd et al., 2017; Verd et al., 2018). In these models, the spiral geometry of the trajectories is generated by a corresponding movement of a point attractor in space, while in the full model it is the consequence of (much less pronounced) attractor movement in addition to the complex eigenvalues of the spiral sink. In other words, there is a certain degeneracy or disconnect between the attractors present in phase space and the resulting transient dynamics driven by the system.

This in turn has important implications for the way in which we analyse dynamical systems models in biology. To understand the dynamics of a regulatory network, it is not generally sufficient to perform a steady-state analysis of a time-invariant or autonomous version of the system. Transient dynamics, and the explicit dependence of regulatory structure on time implied by time-variance or non-autonomy, should be considered the default assumption. Steady-state dynamics and time-invariant autonomy must first be established before conclusions from classical attractor analysis can be considered valid and applicable.

In the second part of our analysis, we focus on the consequences of dynamical modularity and criticality on the evolvability of the system. Our results shed light on the apparent paradox that gap gene expression dynamics are surprisingly insensitive to maternal gradient concentrations during development within a species, but quite labile when comparing different species on an evolutionary time scale (Figures 4 and 5) (Briscoe and Small, 2015). We show that the paradox is resolved if we consider the fact that this lability depends on changes in downstream regulatory interactions between gap genes rather than evolutionary changes to maternal gradients.

Bifurcation analysis reveals that AC/DC1 and AC/DC3 are structurally stable in D. melanogaster, and therefore insensitive to changes in the parameter values that affect the strength of interactions. In contrast, AC/DC2 is in a state of criticality, close to a bifurcation boundary, with regard to variation in the strength of the repressive effect of Kr on hb. Increasing repression between these two genes leads to a bifurcation event, which changes the dynamical regime of AC/DC2 from multistable, switch-like, behaviour (generating stable gap domain boundaries to spiralling trajectories (generating shifting transient boundaries). Evolvability of the gap gene system therefore depends on the fact that some dynamical modules are critical, while others are structurally stable. Evolution induces changes in the stability of specific modules. This explains why the extent and dynamics of gap domain shifts appear to be highly evolvable, while the basic staggered arrangement of domain boundaries remains remarkably stable, at least during the evolution of cyclorrhaphan (or ‘higher’) flies (Bonneton et al., 1997; Stauber et al., 2000; Shaw et al., 2002; Lemke et al., 2008; Lemke and Schmidt-Ott, 2009; Lemke et al., 2010; Wotton et al., 2015; Crombach et al., 2016, see Jaeger, 2011 for review).

Our analysis of D. melanogaster allows us to infer a number of characteristics of gap gene regulation in other dipteran species, such as M. abdita, even though the quality of the models we have for that species unfortunately does not allow a direct comparison of phase spaces (Crombach et al., 2012; Crombach et al., 2016). Dynamic boundary shifts extend much further anterior in this species than in D. melanogaster (Figure 5A) (Wotton et al., 2015). The anteriorly displaced position of the bifurcation suggests that AC/DC1, not AC/DC2, must be critical in M. abdita. Our analysis reveals that changes in several gap-gap interactions are required to render AC/DC1 structurally unstable (Figure 5). This provides a plausible lineage explanation (Calcott, 2009) (or evolutionary trajectory) for changes in dynamical modules that accurately fit what we know about changes in gap gene expression and regulation between M. abdita and D. melanogaster (Wotton et al., 2015; Crombach et al., 2016).

Outside the cyclorrhaphan lineage, the arrangement of gap domains and the patterning output of the system changes. In the moth midge C. albipunctata, for example, there is no posterior domain of Gt protein expression, and the posterior domain of hb mRNA only forms after gastrulation (Rohr et al., 1999; García-Solache et al., 2010). Gap domain shifts are much more pronounced in this species compared to D. melanogaster or M. abdita (García-Solache et al., 2010; Crombach et al., 2014). The beetle T. castaneum shows an even more dynamic mode of segment determination: only the anterior-most domains form simultaneously, while more posterior domains are generated sequentially by sustained oscillations in segmentation gene expression (Sarrazin et al., 2012; El-Sherif et al., 2012).

An extended analysis of our AC/DC circuits (see Appendix 1) reveals that they can be induced to drive sustained limit-cycle oscillations with relatively small additional changes in the values of parameters that determine cross-regulatory interactions (summarized in Figure 6) (Panovska-Griffiths et al., 2013; Perez-Carrasco et al., 2018; Page and Perez-Carrasco, 2018). Although gap genes do not seem to be directly involved in the process of segment determination in T. castaneum, they do show repeated waves of kinematically shifting gene expression in the blastoderm and germband of the embryo (Zhu et al., 2017). So do pair-rule genes (Sarrazin et al., 2012; El-Sherif et al., 2012; El-Sherif et al., 2014), which are known to be essential for segment determination in all arthropods. Pair-rule genes in T. castaneum may also be regulated by AC/DC-like regulatory subcircuits (Choe et al., 2006). These surprising resemblances suggest that ancestral gap and pair-rule expression may have relied on similar regulatory principles, and that the regulatory changes required to turn sequential (short-germband) into simultaneous (long-germband) segmentation may be much more subtle than commonly thought (see Tautz, 2004; Clark, 2017). Improved empirical and modelling evidence from many additional species will be required to rigorously test this prediction.

AC/DC subcircuits and their possible role in the evolution of long-germband from short-germband segment determination.

(A) Phylogenetic relationships between the intermediate-germband insect T. castaneum and the two long-germband dipteran species analysed in this paper. Coloured bars indicate which dynamical regimes are active in which relative region along the A–P axis of the embryo (see key). (B) Visualisation of the dynamical regimes of an AC/DC circuit. We combine equation parameters into two composite control parameters which correspond to the x and Γ~3 parameters in Equations (28) and (29) of Appendix 1 and represent the strength of the positive feedback between two of the genes in the circuit (X-axis) and the strength of the negative feedback involving all three genes (Y-axis) respectively. We used simplified connectionist models (see ‘Materials and methods’)—with all degradation rates set to equal values, time-constant basal activation terms, and no auto-activation—to evaluate the dynamical regimes for varying values of these two control parameters. The results are colour-coded as indicated in the key. Possible dynamical regimes are: monostability and multistability without (grey), and with damped oscillations (light/dark red), as well as sustained limit-cycle oscillations (bright red) (see Appendix 1 for details).

https://doi.org/10.7554/eLife.42832.009

Taken together, our results demonstrate the potential of analysing the developmental and evolutionary dynamics of regulatory networks using an approach based on dynamical rather than structural modularity. The main limitation of this approach remains the small number of systems in which this type of analysis is currently possible. While structural analysis only requires the mechanistic decomposition of a system into its components and their interactions, our analysis relies on the ‘recomposition’ or reconstitution of the system with a dynamical model of some sort (Bechtel and Abrahamsen, 2005; Bechtel and Abrahamsen, 2010; Bechtel, 2011). This model can be discrete and qualitative as in the Boolean analysis carried out by Irons and Monk (2007), or it can be continuous and quantitative as in our example of the gap gene system. In the former case, we require good qualitative evidence on the expression of the relevant factors and their interactions; in the latter, we need quantitative sets of expression data suitable for model fitting (Jaeger and Crombach, 2012). Significant advances in ‘omics’ and ‘big-data’ techniques, genetic perturbation by gene knock-down, precise genome editing, as well as the methodologies of quantitative microscopy and image bioinformatics should render the acquisition of such evidence feasible in a wide range of cellular and developmental systems. The most important consideration for these efforts is that empirical and modelling efforts are tightly integrated and adjusted to each other for an accurate representation and analysis of systems dynamics.

One peculiarity of our gap gene example is the fact that we are dealing with a spatio-temporal patterning system. Accordingly, the pragmatic approach for identifying dynamical modules we present here relies on the fact that different subcircuits are relevant in different regions of an embryo or tissue (see Figure 2). Clearly, this approach is not appropriate for regulatory systems without spatially differentiated dynamics, and in these cases a different strategy is needed. There are two possibilities worth pursuing in the regard. Analogous to the spatial case, node sensitivity analysis can be performed to identify network nodes that are not required during certain periods of time. This subdivides the system into modules active during different periods of influence, rather than regions of influence. However, this approach is limited since it will miss modules that simultaneously drive different dynamical regimes in a multistable network. To overcome this limitation in Boolean systems (whether spatial or not), the approach by Irons and Monk (2007) can be straightforwardly applied. It is based on mapping system components to simple and fundamental dynamical patterns (called subsystems) that combinatorially compose the attractors of a system (see also Jaeger and Monk, 2019). This general strategy could also be used as a guide to identifying dynamical modules in continuous systems: what are the qualitative distinguishable behaviours of the system? What are the simplest dynamical patterns they share? Which network components contribute to which of these patterns? A general strategy of this kind fits well with the type of approach we employed here in the spatial context of the gap gene system.

It is important to note that, in continuous systems, dynamical modularity will always be a matter of degree, just as in the case of structural and functional modules (Wagner and Altenberg, 1996). Subsystems interact and dynamical regimes depend on each other to varying degrees (see Koseska and Bastiaens, 2017), for an interesting perspective on the importance of this phenomenon). This is no major impediment to our approach. As long as the behaviour of a subsystem can be distinguished from other types of dynamics, we can map it to those network components and their interactions that are mainly responsible for driving it. Our analysis of the gap gene network nicely illustrates this point, as the individual AC/DC subcircuits reproduce the main features of the full system’s dynamics, while subtle differences, caused by interactions between subsystems, remain detectable (see Results and Discussion).

Finally, we must ask ourselves whether our approach is scalabale to larger networks and, if so, to what extent. Scalability in this case crucially depends on the availablity of a dynamical model for the full system. If there is a Boolean model, the approach by Irons and Monk (2007) can be used without specific size limitations (see, for example, the networks described in Kauffman, 1993). For larger-scale continuous models, for instance genome-scale dynamic models of metabolism (Smallbone and Mendes, 2013; Srinivasan et al., 2015), what is needed is the ability to identify all (or at least the most prominent or biologically relevant) attractors and dynamical regimes. The numerical analysis of such models can be challenging but, given enough computing power, there is no reason to assume that it cannot be scaled to such larger network models. Indeed, since our approach depends on attractors and dynamical regimes, it is important to note that these increase in number only slowly with an increase in the number of nodes in an interaction network (Kauffman, 1993).

Luckily, all of the limitations discussed above are of a purely practical nature and the number of potential systems amenable to dynamical analysis is increasing. AC/DC circuits involved in neural tube patterning in vertebrates can be considered dynamical modules (Balaskas et al., 2012; Panovska-Griffiths et al., 2013; Perez-Carrasco et al., 2018). So can the growing number of experimentally verified Turing-type pattern generators, for example, those driving digit patterning in the vertebrate limb from sharks to mammals (Raspopovic et al., 2014; Onimaru et al., 2016). We provide a generalised account of dynamical modularity with additional examples in Jaeger and Monk (2019). Apart from these examples, there are good theoretical reasons to believe that dynamical modularity is a very widespread phenomenon. The in silico screen performed by Jiménez et al. (2017) suggests that many regulatory systems are capable of multistable behaviours, while different dynamical regimes rarely map to cleanly separable clusters in the structure of the network. In other words, most regulatory networks are of emergent rather than hybrid type, exhibiting dynamical but not necessarily structural modularity. Based on this, we conclude that the gap gene system is probably not an isolated example of a system where our approach is useful. As a matter of fact, we see an urgent and growing need to identify and characterise dynamical modules in development and evolution (Jaeger and Monk, 2019). The theory and the approach we present here greatly extend the reach of traditional methods by capturing modularity in systems that show no overt clusters in regulatory structure or co-expression patterns. They bring us closer to the aim of identifying true functional modules in evolving developmental systems (Wagner, 2014), as dynamical behaviour is much more tightly integrated with biological function than the regulatory structure of a network.

Materials and methods

The full model: a diffusion-less gap gene circuit

What we refer to as the ‘full model’ in this paper corresponds to a diffusion-less gap gene circuit, which is formulated in the connectionist modelling framework first proposed by Mjolsness et al. (1991). It is derived from gap gene circuits with diffusion (Ashyraliyev et al., 2009), and was previously published and described (Verd et al., 2017; Verd et al., 2018). Here, we only provide a brief description of the model. See these previous publications for details.

Gap gene circuits consist of a one-dimensional row of nuclei, arranged along the antero-posterior (A–P axis of the embryo. They are hybrid models that implement continuous dynamics during interphase and mitosis, with discrete instantaneous nuclear divisions occurring at the end of each mitosis. The spatial domain of the model used here extends over a range of 35 to 75% A–P position (where 0% is the anterior pole), covering the trunk region of the embryo. The full model includes cleavage cycles C13 and C14A of the blastoderm stage during early development of D. melanogaster (Foe and Alberts, 1983). C14A is further subdivided into eight time classes of equal duration (T1–T8) (Surkova et al., 2008). Division takes place at the end of C13.

The state variables of the model represent the concentrations of transcription-factor proteins encoded by trunk gap genes hb, Kr, kni, and gt. gia(t) represents the concentration of protein a in nucleus i at time t. The rate of change in gap protein concentration over time is given by the following system of ordinary differential equations:

(1) dgia(t)dt=Raϕ(ua)-λagia(t),

where Ra is the rate of protein production, and λa the rate of protein decay. ϕ is a sigmoid regulation-expression function which is used to represent the coarse-grained saturating kinetics of transcriptional regulation. It is defined as follows:

(2) ϕ(ua)=12(ua(ua)2+1+1),

where

(3) ua=bGWbagib(t)+mMEmagim(t)+ha.

G={ℎ𝑏,𝐾𝑟,𝑘𝑛𝑖,𝑔𝑡} denotes the set of trunk gap genes, and M={Bcd,Cad} the set of external inputs from protein gradients encoded by maternal coordinate genes (which are not themselves regulated by gap genes. We linearly interpolate quantified spatio-temporal protein expression data (Surkova et al., 2008; Pisarev et al., 2009; Ashyraliyev et al., 2009) to obtain the concentration profiles of the maternal regulators gim.

Interconnectivity matrices W and E, with elements wba and ema, represent regulatory weights of interactions between gap genes, and external inputs from maternal gradients, respectively. The effect of regulator b or m on gap gene target a is activating, if the corresponding weight is positive, repressive, if the weight is negative; there is no interaction if the weight is (near zero. ha is a threshold parameter encoding the basal activity of gap gene a in the absence of any spatially or temporally specific regulators. The system of Equations (1) governs regulatory dynamics during interphase. Ra is set to zero during mitosis.

Earlier studies have established that diffusion of gap proteins is not essential for pattern formation (Jaeger et al., 2004b; Manu et al., 2009a; Verd et al., 2017). Therefore, it is not included in this version of the model. Omitting diffusion renders each nucleus independent of the others, and reduces the dimensionality of system from 164 (4 gene products in 41 nuclei) to 41 independent systems with 4 dimensions each. This makes the system amenable to phase (or state) space analysis.

Fitting of the full model

Request a detailed protocol

Values for parameters Ra, λa, W, E, and ha were obtained by fitting the model to quantitative spatio-temporal gene expression data as previously described (Reinitz and Sharp, 1995; Jaeger et al., 2004b; Ashyraliyev et al., 2009; Verd et al., 2017). Briefly: we solve model Equations (1) numerically, and compare the resulting model output to data by calculating a weighted root-mean-square (RMS) score (Ashyraliyev et al., 2009). This is repeated while changing parameter values until the fit is no longer improving. The difference between model output and data is minimized using a global optimization algorithm called parallel Lam Simulated Annealing (pLSA) (Chu et al., 1999). Model fitting was carried out on the Mare Nostrum cluster at the Barcelona Supercomputing Centre (http://www.bsc.es). The circuit used here is identical to that published and described in Verd et al. (2017); Verd et al. (2018). It accurately reproduces the spatio-temporal expression dynamics of gap gene expression (see Figure 1B), and is fully consistent with the available experimental evidence on gap gene regulation (see Jaeger et al., 2004b; Jaeger et al., 2004a; Manu et al., 2009b; Manu et al., 2009a; Ashyraliyev et al., 2009; Jaeger, 2011; Crombach et al., 2012; Verd et al., 2017; Verd et al., 2018).

Identifying dynamical modules

The full gap gene circuit described above drives two distinct dynamical regimes anterior and posterior of a bifurcation, which occurs at 52% A–P position (Verd et al., 2017; Verd et al., 2018). Anterior to the bifurcation, the system is multistable, positioning static gap domain boundaries through switch-like dynamic behaviour; posterior to the bifurcation, the system is monostable, and the only attractor present is a spiral sink implementing a damped oscillator mechanism driving the observed dynamic anterior shifts of posterior gap domains (Verd et al., 2017; Verd et al., 2018). In the absence of any evident structural modules in the network (see the ‘Introduction’), we ask whether there are dynamical modules or subcircuits that suffice to reproduce the observed dynamical regimes in different regions of the embryo. We define dynamical modules as subcircuits embedded in the gap gene network that are capable of recovering the dynamics of gap gene expression at a given A–P position (Irons and Monk, 2007). Dynamical modules may show overlap in their components, interactions, and regions of influence.

Node sensitivity analysis

Request a detailed protocol

We observe that any specific nucleus along the A–P axis of the embryo only ever expresses two gap genes at the same time, and never more than three different gap genes over the whole duration of cleavage cycle C14A. This implies that only a subset of the four trunk gap genes are required for patterning at any given position. To identify these subsets and their regions of influence, we assessed the sensitivity of simulated developmental trajectories to the removal of gap genes in all nuclei between 35 and 75% A–P position. Early gap gene patterning is largely governed by regulatory inputs from maternal gradients, while gap-gap interactions become increasingly predominant during late C13 and C14A. Since our analysis focuses on gap-gap cross-regulation, we limit our analysis to C14A.

Using maternal gradient and simulated gap gene concentrations at the onset of C14A (time class T1 as initial conditions, we numerically solve the full model, and compare it to simulations where one of the gap genes and all its regulatory interactions have been erased from the model. The ‘node sensitivity’ of the system in a given nucleus with respect to a particular deleted gap gene (or node of the network) is then given by the distance between the trajectory simulated using the full model (curve a), and the trajectory obtained from a simulation lacking the node representing that gene (curve b). Differences are summed over all time classes (T1–T8) in C14A, which results in the following distance metric:

(4) d=i=18(HbTia-HbTib)2+(KrTia-KrTib)2+(KniTia-KniTib)2+(GtTia-GtTib)28.

(HbTij,KrTij,KniTij,GtTij) represent simulated concentrations of gap genes at time class Ti indicates whether the concentration is derived from the full model (a) or simulations with specific nodes removed (b). When a particular node is removed the concentration of that gene in the trajectory simulated from the resulting AC/DC is . The smaller the distance between trajectories, the smaller the contribution of the removed gap gene to the overall dynamics of gene expression. Sensitivities close to zero indicate that the removed gene is not required for patterning in a given nucleus during C14A. The resulting regions of insensitivity are shown in Figure 2A. The corresponding AC/DC subcircuits associated which each of these regions are shown in Figure 2B and C.

AC/DC subcircuits

Formulation and simulation

Request a detailed protocol

We implemented models of the three subcircuits consisting of the components and interactions shown in Figure 2C using the gene circuit modelling formalism described above. Parameter values, including time-variable maternal inputs and autoregulatory weights are taken from the full model. AC/DC1 includes hb, Kr, and gt, spanning the region of 35 to 47% A–P position; AC/DC2 includes hb, Kr, and kni, spanning the region of 49 to 59% A–P position; AC/DC3 includes Kr, kni, and gt, spanning the region of 61 to 75% A–P position. There is a fourth AC/DC subcircuit, consisting of hb, kni, and gt (AC/DC4, whose region of influence lies further posterior, outside the spatial domain of the full model. The AC/DC4 subcircuit was not analysed further in this study.

We use maternal gradients and simulated gap protein concentrations from the full model at C14-T1 as initial conditions. Simulations of AC/DC subcircuits include C14A only (no mitosis or division. During C14A, zygotic gap gene cross-regulation has become essential for the positioning of domain boundaries, largely supplanting initial positional cues by gradients of maternal activators (Jaeger et al., 2007; Jaeger, 2011).

Phase space analysis

Request a detailed protocol

Since they implement time-variable maternal inputs, AC/DC subcircuits are time-variant or non-autonomous dynamical models (Strogatz, 2015). We previously developed a methodology to characterise and classify simulated transient trajectories in such models (Verd et al., 2014), which we have used to analyse the full gap gene circuit model (Verd et al., 2017; Verd et al., 2018). Here, we redeploy this method for the comparative analysis of AC/DC subcircuits and the full model. Briefly, we calculate instantaneous phase portraits for different time points by ‘freezing’ parameter values (maternal inputs) to the value at each given time. Steady states were calculated using a Newton-Raphson algorithm as in Manu et al. (2009a). For each nucleus within an AC/DC subcircuit’s region of influence, we compare trajectories and attractor positions between the subcircuit and the full model. The results are shown in Figure 3, and Figure 3—figure supplement 2, and are discussed in further detail the ‘Results and discussion’ section of our paper. A more systematic analysis of dynamical regimes that can be driven by the AC/DC subcircuit is provided in Appendix 1.

Appendix 1

Mathematical analysis of the AC/DC circuit

The aim of this analysis is to understand the overall dynamic repertoire of the AC/DC circuit. More specifically, we want to understand the dynamical regimes (as defined in the main text that this circuit can implement with different sets of parameters determining the strength of its repressive regulatory interactions. We start from the well-studied behaviour of the three-gene repressilator circuit (Elowitz and Leibler, 2000; Müller et al., 2006; Garcia-Ojalvo et al., 2004; Strelkowa and Barahona, 2010), and then extend the analysis to the AC/DC network (Panovska-Griffiths et al., 2013).

For both repressilator and AC/DC circuit, we characterize the number and type of steady states present for any given set of parameter values. From this, we establish the range of different dynamical regimes that each of these circuits can implement. This comparative approach illustrates how the addition of the backward repressive interaction that distinguishes the AC/DC from the repressilator circuit affects the dynamical repertoire of the network. This analysis reveals that the AC/DC circuit (unlike the repressilator can implement all the dynamical regimes of gene expression that have been observed during segment determination in different groups of insects.

Analysis of the repressilator circuit

Let us first consider the three-gene repressilator circuit shown in Appendix 1—figure 1A. Transcription factor X represses gene y, whose product, transcription factor Y, represses gene z, whose product, transcription factor Z, in turn represses gene x. We represent repressive regulatory interactions by bounded, monotonically decreasing, and continuous functions G. For example, any sigmoid-shaped regulation-expression function conforms to these general conditions.

Appendix 1—figure 1
Analysis of the repressilator circuit.

(A) Structure of the repressilator circuit: the nodes of the network represent genes and their products X, Y, and Z (circles). Repressive interactions among these nodes are indicated by T-bars. We assume constitutive activation of all three nodes. (B) Uniqueness of solutions: solutions of Equation (7) are defined by intersections of the curves given by H(X)=X and H(X)=gxgzgy(X). In this particular case, the shape of both curves results in them intersecting at only one value of X, X*, and therefore there is only one solution. (C) Complex roots of (λ+d)3 is a real negative number, therefore it has one real negative, and two complex roots: λ1, λ2 and λ3. (D, E) Two different sets of possible eigenvalues: (D) Eigenvalues of a saddle point. One of the eigenvalues associated with the steady state has a negative real part, while the other two have positive real parts. In this case, the steady state is a saddle point. The uniqueness of the steady state, and the nature of the regulation-expression functions we’ve chosen, allows us to infer the presence of a limit cycle driving stable oscillations. (E) Eigenvalues of a stable spiral sink. All eigenvalues associated with the steady state have negative real parts and two of them are complex. In this case, the steady state is a stable spiral sink, which drives damped oscillations.

https://doi.org/10.7554/eLife.42832.012

We use the following abstract formulation of the repressilator circuit, where G stands for a regulation-expression function conforming to the conditions prescribed above, and d are gene product degradation rates:

(5) X˙=Gx(Z)dxXY˙=Gy(X)dyYZ˙=Gz(Y)dzZ.

This system has a unique steady state. Defining gi1diGi, for i={X,Y,Z}, the steady states of the system are solutions of the following system of equations:

(6) 0=gx(Z)X0=gy(X)Y0=gz(Y)Z.

By solving for the state variables in the system, and substituting, we obtain:

(7) X=gxgzgy(X).

Due to the fact that the three regulation-expression functions are bounded, monotonically decreasing, and continuous, so is their composition. Therefore, there is a unique positive solution X* to Equation (7) (Appendix 1—figure 1B). By substituting backwards, following analogous reasoning, we can establish that there is a unique positive steady state given by (X*,Y*,Z*).

Types of steady states: eigenvalue analysis

In order to evaluate the nature of the steady state solution of the repressilator model (Equation 7), we linearise the system around its steady state: X=X*+x, Y=Y*+y, Z=Z*+z, and evaluate its eigenvalues (Strogatz, 2015; Hirsch et al., 2012). This linearisation at steady state (X*, Y*, Z*) is written in matrix form as:

(8) [x˙y˙z˙]=[-dx0-γx-γy-dy00-γz-dz][xyz],

where

(9) γx=dGxdZ|Z,γy=dGydX|X,γz=dGzdY|Y.

The γi are all positive due to the fact that Gx, Gy and Gz are monotonically decreasing continuous functions; they measure the sensitivity of the regulation functions at the steady state (so, for example γx measures the sensitivity of the rate of production of the product of X to changes in its regulator Z.

The local stability of the steady state (X*,Y*,Z*) is determined by the nature of the eigenvalues λ associated with the linearised system at the steady state. These are the roots of the characteristic polynomial given by

(10) |-dx-λ0-γx-γy-dy-λ00-γz-dz-λ|=0
(11) (dxλ)(dyλ)(dzλ)γxγyγz=0
(12) (λ+dx)(λ+dy)(λ+dz)=γxγyγzΓ3<0,

where we define Γ to be a negative real number.

For simplicity, we consider the special case where all degradation rates are equal, dx=dy=dz=d:

(13) (λ+d)3=-Γ3.

The three distinct complex roots of (λ+d)3=-Γ3 are given by (see also Appendix 1—figure 1C):

(14) λ+d=(Γ3)13={ΓΓe±iπ3

It follows that the three eigenvalues associated with the linearised system are

(15) λ1=Γd,λ2=Γeiπ3d,λ3=Γeiπ3d,

λ1 is always negative and real. However, the real parts of λ2 and λ3 can take either sign, since

(16) Re(λ2,3)=Γ2-d.

Proof:

(17) λ2=Γeiπ3d=Γ(cos(π3)+isin(π3))dRe(λ2)=Γcos(π3)d=Γ2d,λ3=Γeiπ3d=Γ(cos(π3)isin(π3))dRe(λ3)=Γcos(π3)d=Γ2d.

This leaves us with two alternative possibilities. From the equations in Equation (17) we know that Re(λ2,3)>0 if Γ2-d>0. From this, it follows that Re(λ2,3)>0 if Γ>2d (this case is illustrated in Appendix 1—figure 1D). A steady state with one negative eigenvalue and two eigenvalues with positive real parts is an unstable saddle point. Since the repressilator only has one steady state, this saddle point is unique.

By applying the Poincaré-Bendixson theorem, we can establish that the system must have a limit cycle (Strogatz, 2015; Hirsch et al., 2012). The Poincaré-Bendixson theorem in two dimensions predicts the presence of limit cycles. It tells us that if we can find a trapping region—a closed and bounded region of the phase plane that contains a trajectory that remains within this region from a certain time onwards—and this trapping region contains no equilibrium points, then there must be at least one limit cycle within this region. The Poincaré-Bendixson theorem does not normally hold for systems of dimension larger than two. However, (Mallet-Paret and Smith, 1990) have shown that, in monotone cyclic feedback systems, the ω-limit set of any bounded orbit can be embedded in 2. Under these conditions, which are met in our case, and in this two-dimensional subspace, the Poincaré-Bendixson theorem applies, and we can infer that there must be a limit cycle present. In this case, therefore, the dynamical regime of the repressilator is to drive stable limit-cycle oscillations.

There is a second possibility. If Γ<2dRe(λ2,3) (Appendix 1—figure 1E). The fact that all three eigenvalues have negative real part, with one real and one complex conjugate pair of eigenvalues, implies that the unique steady state is a stable spiral sink. Spiral sinks are the hallmark of damped oscillations. In this case, the dynamical regime of the repressilator is to drive damped oscillations.

A third, very unlikely, possibility occurs when Re(λ2,3)=0. In this case, linear stability analysis suggests that the steady state has a centre in a two dimensional subspace. However, further analysis of the full nonlinear system is required to determine the true nature of the steady state in this case.

In summary, our analysis shows that the repressilator circuit only has one unique steady state. Depending on the parameter values of the model, this steady state is either a saddle point or a spiral sink. Therefore, the dynamical repertoire of this repressilator consists of two distinct dynamical regimes: stable limit-cycle oscillations, or damped oscillations (in agreement with Page and Perez-Carrasco, 2018).

Analysis of the AC/DC circuit

We obtain an AC/DC circuit (Panovska-Griffiths et al., 2013) by adding an additional backward repressive interaction from Z to Y to the repressilator (compare Appendix 1—figure 1A with Appendix 1—figure 2A). The mathematical formulation of the AC/DC circuit is then as follows:

(18) X˙=Gx(Z)dxXY˙=Gy(X,Z)dyYZ˙=Gz(Y)dzZ.
Appendix 1—figure 2
Analysis of the AC/DC circuit.

(A) Structure of the AC/DC circuit: the nodes of the network represent genes and their products X, Y, and Z (circles). T-bars represent the repressive interactions between the genes, where A,B,C and D denote the strengths of those respressions. Repressive interactions among these nodes are indicated by T-bars. We assume constitutive activation of all three nodes. (B) General shape of the function given by F(Y). (C–E) Effect of increasing the back reaction on the number of steady states. Functions as defined in the text. (C) No back reaction: there is only one intersection point. (D) Weak back reaction: still, there is only one intersection point. (E) Strong back reaction: now there are three intersection points. (F) The simplified characteristic Equation (27) defines a depressed cubic. Its dependence on parameter values is shown. (G,H) Roots of the depressed cubic given by Equation (27). (G) Option 1: one real negative root, and two complex roots with positive real part. (H) Option 2: real roots, where one has negative real part and two have positive real parts. (I–M) Combinations of eigenvalues for steady states in the AC/DC circuit. (I) Unstable spiral sink: a steady state with one real negative eigenvalue and two complex eigenvalues with positive real parts. (J) Stable spiral sink: a steady state with one real negative eigenvalue and two complex eigenvalues with negative real parts. (K) Point attractor: a steady state with real negative eigenvalues. (L) Saddle1,2: a steady state with real eigenvalues where one is positive and two are negative. (M) Saddle2,1: a steady state with real eigenvalues where two are positive and one is negative.

https://doi.org/10.7554/eLife.42832.013

Regulation-expression functions Gx and Gz are assumed to be bounded, monotonically decreasing, and continuous. Gy(X,Z) will typically be taken as the multiplication of the regulation-expression function of X on Y, and that of Z on Y, and its form is explained in more detail below.

Let us begin by finding the steady states of the system:

(19)  0=Gx(Z)dxX0=Gy(X,Z)dyY0=Gz(Y)dzZ}{X=1dxGx(Z)=1dxGx(1dzGz(Y))=1dxF(Y)Y=1dyGy(X,Z)=1dyGy(1dxF(Y),1dzGz(Y))Z=1dzGz(Y)

Y* depends on the strength of the backward reaction from Z on Y, as given by Gz(Y*), and on the combination of the reactions of Z on X and X on Y represented by F(Y) (Equation (19)) is monotonically increasing, as it results from the composition of two monotonically decreasing functions, with a general shape as shown in Appendix 1—figure 2B. As the backward reaction on Y is increased, the three scenarios shown in Appendix 1—figure 2C–E become possible. Let us now look at the Jacobian matrix associated with the system:

(20) [-dx0-γx-γyx-dy-γyz0-γz-dz]

where

(21) γx=dGxdZ|steady state,γz=dGzdY|steady state>0,

and

(22) γyx=dGydX|steady state,γyz=dGydZ|steady state>0.

γyx measures the sensitivity of the repressive interaction of X on Y, and γyz, to the repressive interaction of Z on Y. Again, for the sake of simplicity, let us assume that dx=dy=dz=d and calculate the associated characteristic equation:

(23) |-d-λ0-γx-γyx-d-λ-γyz0-γz-d-λ|=0

From this, it follows that

(24) (dλ)3+(γx)(γyx)(γz)(dλ)(γyz)(γz)=0,
(25) (d+λ)3+γxγyxγz(d+λ)γyzγz=0.

If we let d+λ=η and Γ~3=γxγyxγz>0, then Equation (25) becomes

(26) η3-γyzγzη+Γ~3=0.

Now, let x=γyzγz>0, and Equation (26) becomes

(27) η3-xη+Γ~3=0.

This implies that the AC/DC circuit has either one or three steady states.

Types of steady states: eigenvalue analysis

The eigenvalues of all steady states are roots of the polynomial in Equation (27). In particular, x and Γ~ are evaluated at the steady state. They remain positive since

(28 x=γyz|ssγz|ss>0
(29) Γ~3=γx|ssγyx|ssγz|ss>0.

The simplified characteristic equation shown in Equation (27) defines a depressed cubic (Appendix 1—figure 2F). Its exact shape depends on the values of x and Γ~ when evaluated at each steady state. Given that these parameters always remain positive, we can look at the shape of Equation (27), and from it infer the possible roots, as well as the type and sign of every eigenvalue.

From Appendix 1—figure 2F we can see that, when solving for η, Equation (27) will always have a negative real root and, depending on the magnitude of x, either two positive real roots, or two conjugate complex roots, both with positive real parts (Appendix 1—figure 2G, H). The roots of the depressed cubic are given by ηi,i{1,2,3}, where ηi=λi+di and λi,i{1,2,3}, are the eigenvalues associated with the corresponding steady state of the AC/DC circuit. If we consider di=d, i, we get five possible combinations of eigenvalues, indicating that steady states found in the AC/DC circuit can be unstable spiral sinks (Appendix 1—figure 2I), stable spiral sinks (Appendix 1—figure 2J), point attractors (Appendix 1—figure 2K), saddles1,2 (Appendix 1—figure 2L), or saddles2,1 (Appendix 1—figure 2M).

In summary, the AC/DC circuit can exhibit dynamical regimes that consist of one or three steady states from among the five types described above. This indicates a rich dynamical repertoire, which includes the regimes shown below, known to be involved in segment determination in different groups of insects.

Appendix 1—figures 46 show phase portraits that have been calculated using a connectionist version of the AC/DC model (see Equation 33 below), where production rates Rx=Ry=Rz=1. Diffusion paramters d, expression thresholds h, and repressive interaction strengths A,B,C,D are shown in the corresponding Appendix 1—tables 13.

  • Sustained oscillations are observed if the underlying phase portrait has a unique unstable spiral, and therefore a limit cycle (Appendix 1—figure 3, parameter values used for simulation are shown in Appendix 1—table 1).

Appendix 1—figure 3
Sustained limit-cycle oscillations: simulation showing a limit cycle (blue curve) driving sustained oscillations in the phase portrait of an AC/DC circuit.

This circuit was simulated using the connectionist formulation given by Equation 33, and the parameter values shown in Appendix 1—table 1.

https://doi.org/10.7554/eLife.42832.014
Appendix 1—table 1
Parameter values used to simulate the limit cycle shown in Appendix 1—figure 3.

See Equation 33 for parameter definitions.

https://doi.org/10.7554/eLife.42832.015
Dynamical regimeParameters
dxdydzHABCD
Sustained Oscillations0.090.090.091.5-0.01-0.01-0.01
Appendix 1—figure 4
Damped oscillations: simulation of the damped oscillatory dynamical regime.

Two spiralling trajectories (blue and green curves) are shown converging towards a stable spiral sink. This circuit was simulated using the connectionist formulation given by Equation 33, and the parameter values shown in Appendix 1—table 2.

https://doi.org/10.7554/eLife.42832.016
Appendix 1—table 2
Parameter values used to simulate the damped oscillations shown in Appendix 1—figure 4.

See Equation 33 for parameter definitions.

https://doi.org/10.7554/eLife.42832.017
Dynamical regimeParameters
dxdydzHABCD
Damped Oscillations0.090.090.091.5-0.025-0.025-0.025
  • In addition, the AC/DC circuit can also exhibit a bistable dynamical regime, where the underlying phase portrait will have two point attractors and a saddle1,2 (Appendix 1—figure 5, parameter values used for simulation are shown in Appendix 1—table 3).

Appendix 1—figure 5
Bistability: simulation of the bistable dynamical regime.

Trajectories with varying initial conditions (coloured curves) are shown to converge towards the two point attractors. This circuit was simulated using the connectionist formulation given by Equation 33, and the parameter values shown in Appendix 1—table 3.

https://doi.org/10.7554/eLife.42832.018
Appendix 1—table 3
Parameter values used to simulate the bistable regime shown in Appendix 1—figure 4.

See Equation 33 for parameter definitions.

https://doi.org/10.7554/eLife.42832.019
Dynamical regimeParameters
dxdydzHABCD
Bi-stability0.090.090.091.5-0.01-0.09-0.01-0.09

We were able to find additional dynamical regimes by sampling the parameter space numerically. These are listed below. Note that a numerical analysis cannot guarantee an exhaustive list of all possible dynamical regimes. The biological relevance of these regimes has not yet been explored.

  • An unstable spiral can co-exist with a saddle and a point attractor.

  • A stable spiral sink can co-exist with a saddle2,1 and a point attractor in a bistable regime. The saddle probably needs to have a 2-dimensional unstable manifold, for there to be a basin of attraction with a stable spiral.

  • A bistable regime with two stable spiral sinks.

We conclude that the dynamical repertoire of the AC/DC circuit contains at least six different monostable and multistable regimes, three of which correspond to behaviour observed in insect segmentation. Similar dynamical regimes have been observed when the ACDC circuit was explored using different regulatory functions (Perez-Carrasco et al., 2018; Page and Perez-Carrasco, 2018).

Validation of the analysis when using a connectionist model formulation

The analysis above uses an abstract general formulation of the AC/DC circuit to characterize its dynamical repertoire. We validate this analysis for the connectionist formalism as described in ‘Materials and methods’, which was used for all the simulations in the main paper.

Let us consider the AC/DC circuit shown in Appendix 1—figure 6A, and this time formulate it mathematically as follows:

(30) X˙=Rxg(u)dxXY˙=Ryg(u)dyYZ˙=Rzg(u)dzZ

where:

(31) g(u)=0.5(uu2+1+1)
Appendix 1—figure 6
Connectionist AC/DC circuit.

(A) Regulatory structure and parameters of the simulated AD/DC circuit. Nodes represent genes and their transcription-factor products. T-bars represent repressive regulatory interactions. All genes are constitutively activated. (B) Sigmoidal regulation-expression function used for the simulations (Equation 31). (C) The number of steady states in the connectionist AD/DC circuit, can be one (black curves), or three (coloured curves) depending on the strength of the back reaction.

https://doi.org/10.7554/eLife.42832.020

The sigmoidal function Equation (31) is shown in Appendix 1—figure 6B. It is bounded, monotonically increasing, and continuous. u represents the sum of all regulatory contributions to a given target gene. For example, the contribution of gene product X on gene y is evaluated as the product between the strength of the interaction of X on y, the sign of that interaction (negative for repression, positive for activation, denoted by δ+-), and the concentration of X. R represents the maximum rate of activation for every gene. More generally, we write u as:

(32) u=i(δ+-)αiXi,

based on which Equation (30) becomes:

(33) X˙=Rxg(CZ+h)dxXY˙=Ryg(AXDZ+h)dyYZ˙=Ryg(BY+h)dzZ

In what follows, we will assume that Rx=Ry=Rz=R, dx=dy=dz=d, and Rd=R¯ . As in the general case explored above, it is possible to have more than one steady state with the connectionist formulation of an AC/DC circuit:

(34)  0=Rxg(CZ+h)dxX0=Ryg(AXDZ+h)dyY0=Rzg(BY+h)dzZ}{X=R¯g(CR¯g(BY+h)+h))=R¯F(Y)Y=R¯g(AR¯g(CR¯g(BY+h)+h)DR¯g(BY+h)+h)=R¯G(Y)Z=R¯g(BY+h)}

Again, F(Y) is monotonically increasing (F(Y)0,Y), since it is the composition of two monotonically decreasing functions (g is increasing with respect to u, but decreasing with respect to X, Y and Z). Depending on the strength of the back reaction—which is given by the magnitude of parameter D in the connectionist model formulation (Appendix 1—figure 6A)—we will have one or three steady states (Appendix 1—figure 6C). This is in agreement with our general analysis above.

We will now show that the characteristic equation is the same in both cases. Once this is established, it follows that the same types of steady states occur in both model formulations. Therefore, the dynamical regimes revealed by our general analysis also apply in the connectionist case. Let us look at the associated Jacobian

(35) [-dx0-γx-γyx-dy-γyz0-γz-dz]

where

(36) [Rxg(CZ+h)dxX]Z=dRxg(u)dududZ=dRxg(u)dud(CZ+h)dZ=γx¯(C)=γx<0
(37) [Rzg(BY+h)dzZ]Y=dRzg(u)dududY=dRzg(u)dud(BY+h)dY=γz¯(B)=γz<0
(38) [Ryg(AXDZ+h)dyY]X=dRyg(u)duuX=dRyg(u)du(AXDZ+h)X=γyx¯(A)=γyx<0
(39) [Ryg(AXDZ+h)dyY]Z=dRyg(u)duuZ=dRyg(u)du(AXDZ+h)Z=γyz¯(D)=γyz<0.

If we now assume that dx=dy=dz=d, we have the same characteristic equation for the connectionist formulation as in the general case (see Equation 23):

(40) |-d-λ0-γx-γyx-d-λ-γyz0-γz-d-λ|=0

This demonstrates that the results of the analysis of a general AC/DC circuit are still valid when the AC/DC circuit is formulated as a connectionist model. Since the characteristic equation is the same in both cases, the dynamical repertoire will be equivalent as well.

References

  1. 1
  2. 2
  3. 3
    An Introduction to Systems Biology: Design Principles of Biological Circuits
    1. U Alon
    (2006)
    Boca Raton: Chapman and Hall/CRC.
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
    Explanation: a mechanist alternative
    1. W Bechtel
    2. A Abrahamsen
    (2005)
    Studies in History and Philosophy of Science Part C: Studies in History and Philosophy of Biological and Biomedical Sciences 36:421–441.
    https://doi.org/10.1016/j.shpsc.2005.03.010
  11. 11
  12. 12
  13. 13
  14. 14
    Epidermal patterning in arabidopsis: models make a difference
    1. M Benítez
    2. NAM Monk
    3. ER Alvarez-Buylla
    (2011)
    Journal of Experimental Zoology Part B: Molecular and Developmental Evolution 316B:241–253.
    https://doi.org/10.1002/jez.b.21398
  15. 15
  16. 16
    The Evolution of Complexity by Means of Natural Selection
    1. JT Bonner
    (1988)
    Princeton: Princeton University Press.
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
    The evolution of evolvability
    1. R Dawkins
    (1989)
    In: C Langton, editors. Artificial Life: The Proceedings of an Interdisciplinary Workshop on the Synthesis and Simulation of Living Systems.  Redwood City: Addison-Wesley. pp. 201–220.
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
    Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesis
    1. VE Foe
    2. BM Alberts
    (1983)
    Journal of Cell Science 61:31–70.
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
    Differential Equations, Dynamical Systems, and an Introduction to Chaos (Third Edition)
    1. M Hirsch
    2. S Smale
    3. R Devaney
    (2012)
    Amsterdam, Netherlands: Elsevier.
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
    The gap gene network
    1. J Jaeger
    (2011)
    Cellular and Molecular Life Sciences 68:243–274.
    https://doi.org/10.1007/s00018-010-0536-y
  65. 65
    The inheritance of process: a dynamical systems approach
    1. J Jaeger
    2. D Irons
    3. N Monk
    (2012)
    Journal of Experimental Zoology Part B: Molecular and Developmental Evolution 318:591–612.
    https://doi.org/10.1002/jez.b.22468
  66. 66
    Life’s Attractors
    1. J Jaeger
    2. A Crombach
    (2012)
    In: O. S Soyer, editors. Evolutionary Systems Biology, 751. New York: Springer. pp. 93–119.
    https://doi.org/10.1007/978-1-4614-3567-9_5
  67. 67
  68. 68
    Evolutionary Systems Biology 2.0
    1. J Jaeger
    2. NA Monk
    (2019)
    Dynamical modularity of the genotype-phenotype map, Evolutionary Systems Biology 2.0, Berlin, Springer.
  69. 69
    On the concept of mechanism in development
    1. J Jaeger
    2. J Sharpe
    (2014)
    In: A Minelli, T Pradeu, editors. Towards a Theory of Development. Oxford University Press. pp. 56–78.
    https://doi.org/10.1093/acprof:oso/9780199671427.001.0001
  70. 70
  71. 71
    The Origins of Order: Self-Organization and Selection in Evolution
    1. SA Kauffman
    (1993)
    Oxford: Oxford University Press.
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
  86. 86
  87. 87
  88. 88
  89. 89
  90. 90
  91. 91
  92. 92
  93. 93
  94. 94
  95. 95
  96. 96
  97. 97
  98. 98
  99. 99
  100. 100
  101. 101
  102. 102
  103. 103
  104. 104
  105. 105
  106. 106
  107. 107
  108. 108
  109. 109
  110. 110
  111. 111
    Genomic Control Process: Development and Evolution
    1. I Peter
    2. EH Davidson
    (2015)
    Amsterdam, NL: Academic Press.
  112. 112
  113. 113
  114. 114
  115. 115
    The Shape of Life: Genes, Development, and the Evolution of Animal Form
    1. RA Raff
    (1996)
    Chicago, IL, United States: University of Chicago Press.
  116. 116
  117. 117
  118. 118
    Die Ordnung des Lebendigen: Systembedingungen der Evolution
    1. R Riedl
    (1975)
    Hamburg, DE: P. Parey.
  119. 119
  120. 120
  121. 121
  122. 122
  123. 123
  124. 124
    Critical Transitions in Nature and Society
    1. M Scheffer
    (2009)
    Princeton: Princeton University Press.
  125. 125
    Modularity in Development and Evolution
    1. G Schlosser
    2. GP Wagner
    (2004)
    Chicago: University of Chicago Press.
  126. 126
  127. 127
  128. 128
  129. 129
  130. 130
    The architecture of complexity
    1. HA Simon
    (1962)
    Proceedings of the American Philosophical Society 106:467–482.
  131. 131
  132. 132
  133. 133
  134. 134
  135. 135
  136. 136
  137. 137
  138. 138
  139. 139
  140. 140
  141. 141
    Structural Stability and Morphogenesis
    1. R Thom
    (1976)
    Reading: Westview Press.
  142. 142
  143. 143
  144. 144
  145. 145
  146. 146
  147. 147
  148. 148
  149. 149
  150. 150
    Homology, Genes, and Evolutionary Innovation
    1. GP Wagner
    (2014)
    Princeton: Princeton University Press.
  151. 151
  152. 152
  153. 153
  154. 154
    Developmental Plasticity and Evolution
    1. MJ West-Eberhard
    (2003)
    Oxford: Oxford University Press.
  155. 155
  156. 156
  157. 157

Decision letter

  1. Lee Altenberg
    Reviewing Editor; The KLI Institute, United States
  2. Patricia J Wittkopp
    Senior Editor; University of Michigan, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Modularity, criticality, and evolvability of a developmental gene regulatory network" for consideration by eLife. Your article has been reviewed by three peer reviewers, who are generally positive about the paper, but raise some important issues that need to be addressed, and the evaluation has been overseen by a Reviewing Editor and Patricia Wittkopp as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

The three reviewers found much to like in this work and I thought the authors would appreciate hearing this in their own words, which I've combined together below. Following these summary comments is a single list of required revisions:

The manuscript "Modularity, criticality, and evolvability of a developmental gene regulatory network" is a theoretical work, based on previously published data and models. It is first concerned with the "modularity" keyword of the title: it presents a method to detect functional modules of interaction networks. The method defines dynamical modules as those where one or more of the nodes of the network are largely negligible and the whole network behaviour can be reasonable well captured by a reduced module without the irrelevant nodes.

Then it applies this method to the gap gene regulatory network of the vinegar fly Drosophila melanogaster, finding it can be simplified into simpler modules along the antero-posterior axis of the fly embryo and using these modules to analyse the system. This analysis draws conclusions on the other two keywords, "criticality and evolvability". It has long been discussed how different kind of networks, including gene regulatory ones, are "critical" or in a state of "criticality". This means very different things in different contexts, form properties of the topology of the networks to dynamical properties. Here a clear definition of criticality is made: criticality is to be "close to a bifurcation point beyond which" the "dynamics may change drastically and abruptly". The study shows that one of the modules of the gap gene network is critical, actually, at both sides of a bifurcation depending on the position along the antero-posterior axis of the fly embryo. This has implications on evolvability: small changes in the critical module might be the source of phenotypic change. To discuss this, D. melanogaster is compared to Megaselia abdita, showing that indeed the developmental differences between these flies can be explained in terms of the critical properties (bifurcations) of their regulatory networks.

The "modularity" part, with the definition of "dynamical modules" of a network as opposed to the "structural modules" found in the literature, is a clever trick that works very well in this systems, but I am afraid it might not be easy to generalise. The model used here is actually 41 independent networks of 4 nodes each. In this setting, it is relevant to find in which of those 41 networks we may leave out one or more nodes. However, in the more frequent modelling situation where we have a single network, this strategy would only single out nodes that are irrelevant and should have not been in the model to start with. It might be useful, however, in non-autonomous systems, to define different dynamical modules as time progresses.

However, I find the "criticality and evolvability" conclusions highly relevant. It does not depend on the modular decomposition of the network: the full model already shows the bifurcation discussed in the AC/DC2 module. However, pointing at the importance of dynamical properties as bifurcations with regard to the evolvability of developmental systems is a strong message, and it is cleanly shown in this example. I believe this way of seeing the underlying mechanics of evo-devo showcases a very generic result, relevant for the whole field, and as such worth of attention and the exposure that a publication in eLife would bring. For this reason I recommend publication.

This manuscript details a very important finding: subnetworks with largely overlapping genes can display functional modularity and evolvability. This conclusion complements the current views in the literature that focus on structural modularity, in which modular networks don't share genes. The results helps explain a common situation, in which genes are implicated in multiple related networks that act at different times, and hence would appear to have many pleiotropic effects. The results illustrate a compelling example of functional modularity, and I found the explanations of rather complex modeling quite clear such that a general readership could follow the logic and main points.

The paper by Verd and collaborators also offers a new angle to the study of modularity in molecular networks. Specifically, it contributes to developing the concept of dynamic modularity and the relation between phenotypes and dynamic modules. The authors were able to formally identify the repetition of a dynamic subcircuit in the gap gene network, and to study the dynamic behavior of each subcircuit in a non-autonomous context. I found particularly relevant their results and discussion on the potential mechanisms for the evolution and evolvability of the gap gene network, for which a comparative AND dynamic approach was followed. The careful analysis of dynamic systems that are partially shared by different lineages is a much needed perspective in evo-devo and is thus an important contribution of the paper.

I would like to discuss the suitability of these "technical" works for multidisciplinary journals. Nobody would comment on how "technical" an experimental effort would be, no matter how obscure its methods and how subtle and prone to misconception the results might be for a reader not absolutely familiar with the experimental technique used. However, theory papers based on mathematical ideas have to fight again and again against the prejudice of being "technical", no matter how clearly the premises and results are discussed and made amenable for a general audience. This has to come to an end if we expect the life sciences to one day be a truly multidisciplinary effort. Only the importance of the question made and the relevance of the conclusions, together with the clarity of presentation, should matter; the method used (experiment or theory) should be irrelevant as long as it can be assessed to be sound and correct. With respect to this, this manuscript is well written and makes its points clearly, and the authors have rightly considered that eLife was a reasonable venue to present their ideas and results.

Essential revisions:

1) My main suggestion is that the authors provide a little more guidance to back up their impact statement, which suggests that the paper provides a road map for others to apply this approach. Typical readers would not necessarily know whether they themselves can apply this approach in other systems. The authors acknowledge that few systems have enough empirical knowledge to apply these models. Is this paper meant to be a cautionary example that separate processes with overlapping genes could still act as functional modules? Or is this meant to be a roadmap to researchers? If the authors intend a roadmap, I would recommend adding to the conclusions a more explicit statement of what empirical data researchers would need in order to apply this approach. In either case, I suggest more extensive discussion of what this means for a reader who doesn't study gap genes in insects. What does non-modularity look like? Summarizing all the parameter sensitivity analyses, is it likely than any overlapping AC/DC circuits would exhibit functional modularity, or is this likely a special case? Is functional modularity a yes-no categorization, or are there degrees of functional modularity (e.g. do the departures from the full model imply some lack of functional modularity, and does that matter?). Adding a bit more to guide the reader here will clarify the broader implications of the work.

2) Another point that could benefit from more basic discussion is the link between criticality and evolvability. Criticality is likely a hard topic for the readers to picture. Perhaps one or more phase portraits could be added to a figure to help illustrate how features of the dynamical system promote evolvability? The Discussion is clear in general, but more high level summary statements could help some readers understand what dynamical systems would reduce evolvability and what would promote evolvability.

3) Another potentially challenging point is the non-autonomous aspect of the models. Perhaps a simpler statement when this first appears can help the reader understand autonomy more quickly and intuitively.

4) As the authors mention, their study model constitutes a very small network. Could they comment or speculate on how robust and/or dependent on the network scale and type of data is the dynamic vs. structural identification of modules?

5) How likely is it to follow a similar, pragmatic, approach for the identification of modules in other systems (it would be useful to point to some examples)? Is this approach not possible only because of the structural or dynamic peculiarities of the gap gene system (one of the two or three best characterized gene networks)? If such pragmatic approach were not widely applicable to other systems, it would be valuable to elaborate on what this case study can contribute to the more general dynamical/structural module study and on potential ways to test these ideas in other systems. For example, how could one test that identical transient behavior can be caused by different types of moving attractors in biological non-autonomous dynamical systems in other networks that are not as well studied.

https://doi.org/10.7554/eLife.42832.022

Author response

Essential revisions:

1) My main suggestion is that the authors provide a little more guidance to back up their impact statement, which suggests that the paper provides a road map for others to apply this approach. Typical readers would not necessarily know whether they themselves can apply this approach in other systems. The authors acknowledge that few systems have enough empirical knowledge to apply these models. Is this paper meant to be a cautionary example that separate processes with overlapping genes could still act as functional modules? Or is this meant to be a roadmap to researchers? If the authors intend a roadmap, I would recommend adding to the conclusions a more explicit statement of what empirical data researchers would need in order to apply this approach. In either case, I suggest more extensive discussion of what this means for a reader who doesn't study gap genes in insects. What does non-modularity look like? Summarizing all the parameter sensitivity analyses, is it likely than any overlapping AC/DC circuits would exhibit functional modularity, or is this likely a special case? Is functional modularity a yes-no categorization, or are there degrees of functional modularity (e.g. do the departures from the full model imply some lack of functional modularity, and does that matter?). Adding a bit more to guide the reader here will clarify the broader implications of the work.

We definitely intended our paper to be a roadmap, the beginning of a new research programme, and not a cautionary tale. We have added four new paragraphs to our “Conclusions” section of the manuscript, which discuss what type of empirical evidence is required to apply our approach, how it can be applied to systems that are not spatially differentiated, and how it scales to the analysis of larger networks (see also points 4 and 5 below).

We emphasise that non-modularity manifests itself in the absence of structural clusters in a multifunctional network. This is a very general criterion to identify systems with dynamical, but not structural, modularity.

We have also added a paragraph which discusses the fact that modularity (whether functional, structural, or dynamical) is always a matter of degree, at least in continuous dynamical systems.

We hope that these discussions clarify the broader implications and applicability of our work.

2) Another point that could benefit from more basic discussion is the link between criticality and evolvability. Criticality is likely a hard topic for the readers to picture. Perhaps one or more phase portraits could be added to a figure to help illustrate how features of the dynamical system promote evolvability? The Discussion is clear in general, but more high level summary statements could help some readers understand what dynamical systems would reduce evolvability and what would promote evolvability.

We have added several clarifying statements (Results and Discussion, and in the Conclusions) about the role of criticality and structural stability in determining the evolvability of a regulatory systems. The discussion, in particular in “Conclusions”, is intended to make the different sensitivities of subsystems understandable to the reader who is not familiar with dynamical systems theory.

3) Another potentially challenging point is the non-autonomous aspect of the models. Perhaps a simpler statement when this first appears can help the reader understand autonomy more quickly and intuitively.

We now introduce the concept of time-variance or non-autonomy more gently. Since “non-autonomy” is a potentially confusing term, we now use “time-variant” along with it throughout the manuscript. This both clarifies the concept and maintains continuity with our earlier publications on non-autonomous systems.

4) As the authors mention, their study model constitutes a very small network. Could they comment or speculate on how robust and/or dependent on the network scale and type of data is the dynamic vs. structural identification of modules?

We have added a paragraph to the “Conclusions” section discussing how our approach scales to larger systems.

5) How likely is it to follow a similar, pragmatic, approach for the identification of modules in other systems (it would be useful to point to some examples)? Is this approach not possible only because of the structural or dynamic peculiarities of the gap gene system (one of the two or three best characterized gene networks)? If such pragmatic approach were not widely applicable to other systems, it would be valuable to elaborate on what this case study can contribute to the more general dynamical/structural module study and on potential ways to test these ideas in other systems. For example, how could one test that identical transient behavior can be caused by different types of moving attractors in biological non-autonomous dynamical systems in other networks that are not as well studied.

As mentioned for point 1 above, we have added four paragraphs to the “Conclusions” section, that discuss how our approach can be applied to other systems, even if they are quite different in nature/scale to the example of the gap gene network which we use in this paper.

https://doi.org/10.7554/eLife.42832.023

Article and author information

Author details

  1. Berta Verd

    1. EMBL/CRG Systems Biology Research Unit, Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Barcelona, Spain
    2. Universitat Pompeu Fabra (UPF), Barcelona, Spain
    3. Konrad Lorenz Institute for Evolution and Cognition Research (KLI), Klosterneuburg, Austria
    4. Department of Genetics, University of Cambridge, Cambridge, United Kingdom
    Present address
    Department of Genetics, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Conceptualization, Resources, Software, Formal analysis, Funding acquisition, Validation, Investigation, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    bertaverd@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9835-009X
  2. Nicholas AM Monk

    School of Mathematics and Statistics, University of Sheffield, Sheffield, United States
    Contribution
    Formal analysis, analysis of the model and revising the article critically
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5465-4857
  3. Johannes Jaeger

    1. EMBL/CRG Systems Biology Research Unit, Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Barcelona, Spain
    2. Universitat Pompeu Fabra (UPF), Barcelona, Spain
    3. Konrad Lorenz Institute for Evolution and Cognition Research (KLI), Klosterneuburg, Austria
    4. School of Mathematics and Statistics, University of Sheffield, Sheffield, United States
    5. Wissenschaftskolleg zu Berlin, Berlin, Germany
    6. Center for Systems Biology Dresden (CSBD), Dresden, Germany
    7. Complexity Science Hub (CSH), Vienna, Austria
    8. Centre de Recherches Interdisciplinaires (CRI), Paris, France
    Present address
    Department of Molecular Evolution and Development, University of Vienna, Vienna, Austria
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    yoginho@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2568-2103

Funding

MINECO (BFU2009-10184/BFU2012-33775/SEV-2012-0208)

  • Johannes Jaeger

European Commission (FP7-KBBE-2011-5/289434 (BioPreDyn))

  • Johannes Jaeger

MEC-EMBL

  • Johannes Jaeger

La Caixa Savings Bank

  • Berta Verd

KLI Klosterneuburg

  • Berta Verd

Wissenschaftskolleg zu Berlin

  • Johannes Jaeger

Max-Planck-Gesellschaft

  • Johannes Jaeger

Herchel Smith Fund

  • Berta Verd

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Anton Crombach for insightful discussions, comments on the manuscript, and for providing the fitted full model used in this analysis. Other members of the Jaeger Lab in Barcelona, as well as fellows and visitors at the KLI Klosterneuburg, provided inspiring and motivating feedback on the project. We thank Erik Clark and Ruben Pérez-Carrasco for countless discussions and for their thoughtful comments on the manuscript. We thank Ovidiu Radulescu for initial discussions about the dynamical behaviour of the full circuit. The authors thankfully acknowledge the computer resources, technical expertise and assistance provided by the Barcelona Supercomputing Center—Centro Nacional de Supercomputación.

Senior Editor

  1. Patricia J Wittkopp, University of Michigan, United States

Reviewing Editor

  1. Lee Altenberg, The KLI Institute, United States

Publication history

  1. Received: January 14, 2019
  2. Accepted: June 5, 2019
  3. Accepted Manuscript published: June 6, 2019 (version 1)
  4. Version of Record published: July 22, 2019 (version 2)

Copyright

© 2019, Verd et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,146
    Page views
  • 445
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Physics of Living Systems
    Andreas Hofmann et al.
    Research Article Updated
    1. Computational and Systems Biology
    2. Neuroscience
    Tung D Phan et al.
    Research Article Updated