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

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.


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.
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). 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 Figure 1 continued on next page 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 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. 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  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, B; Jaeger, 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).
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 nonoverlapping 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 multistability; 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).
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   (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). Timevariance 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 nucleusalmost 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 ( 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 3figure 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  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 1-3), 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-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 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. (A-F) 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. DOI: https://doi.org/10.7554/eLife.42832.008 ( 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 selforganising 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. 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., 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.
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 . 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  (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). DOI: https://doi.org/10.7554/eLife.42832.009 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. g a i ð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: dg a i ðtÞ dt ¼ R a fðu a Þ À l a g a i ðtÞ; where R a is the rate of protein production, and l a the rate of protein decay. f is a sigmoid regulation-expression function which is used to represent the coarse-grained saturating kinetics of transcriptional regulation. It is defined as follows: where G ¼ fhb; Kr; kni; gtg denotes the set of trunk gap genes, and M ¼ fBcd; Cadg 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 g m i . Interconnectivity matrices W and E, with elements w ba and e ma , 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. h a 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. R a 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
Values for parameters R a , l a , W, E, and h a were obtained by fitting the model to quantitative spatiotemporal 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
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: ðHb j Ti ; Kr j Ti ; Kni j Ti ; Gt j Ti represent simulated concentrations of gap genes at time class T i 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
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
Since they implement time-variable maternal inputs, AC/DC subcircuits are time-variant or nonautonomous 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. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. 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Þ ¼ g x g z g y ð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 ðl þ dÞ 3 is a real negative number, therefore it has one real negative, and two complex roots: l 1 , l 2 and l 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. 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:

Author contributions
This system has a unique steady state. Defining g i 1 d i G i , for i ¼ fX; Y; Zg, the steady states of the system are solutions of the following system of equations: By solving for the state variables in the system, and substituting, we obtain: X ¼ g x g z g y ð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: and evaluate its eigenvalues (Strogatz, 2015;Hirsch et al., 2012). This linearisation at steady state (X Ã , Y Ã , Z Ã ) is written in matrix form as: where The g i are all positive due to the fact that G x , G y and G z are monotonically decreasing continuous functions; they measure the sensitivity of the regulation functions at the steady systems, the !-limit set of any bounded orbit can be embedded in R 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 G<2d¼)R e ðl 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ðl 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: Regulation-expression functions G x and G z are assumed to be bounded, monotonically decreasing, and continuous. G y ð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: Y Ã depends on the strength of the backward reaction from Z on Y, as given by G z ð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: Àg z Àd z 2 6 4 3 7 5 where and g yx ¼ À dG y dX steadystate ; g yz ¼ À dG y dZ steadystate >0: g yx measures the sensitivity of the repressive interaction of X on Y, and g yz , to the repressive interaction of Z on Y. Again, for the sake of simplicity, let us assume that d x ¼ d y ¼ d z ¼ d and calculate the associated characteristic equation: Àd À l 0 Àg x Àg yx Àd À l Àgyz 0 Àg z Àd À l ¼ 0 From this, it follows that ðÀd À lÞ 3 þ ðÀg x ÞðÀg yx ÞðÀg z Þ À ðÀd À lÞðÀg yz ÞðÀg z Þ ¼ 0; (24) ðd þ lÞ 3 þ g x g yx g z À ðd þ lÞg yz g z ¼ 0: If we let d þ l ¼ h andG 3 ¼ g x g yx g z >0, then Equation (25) becomes h 3 À g yz g z h þG 3 ¼ 0: Now, let x ¼ g yz g z >0, and Equation (26) becomes h 3 À xh þG 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 andG are evaluated at the steady state. They remain positive since x ¼ g yz j ss g zj ss >0 (28 G 3 ¼ g xj ss g yx j ss g zj 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 andG 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 h, 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 h i ; i 2 f1; 2; 3g, where h i ¼ l i þ d i and l i ; i 2 f1; 2; 3g, are the eigenvalues associated with the corresponding steady state of the AC/ DC circuit. If we consider d i ¼ d, 8i, 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), saddles 1,2 (Appendix 1-figure 2L), or saddles 2,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 4-6 show phase portraits that have been calculated using a connectionist version of the AC/DC model (see Equation 33 below), where production rates R x ¼ R y ¼ R z ¼ 1. Diffusion paramters d, expression thresholds h, and repressive interaction strengths A; B; C; D are shown in the corresponding Appendix 1-tables 1-3.
. 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).
In what follows, we will assume that R x ¼ R y ¼ R z ¼ R, d x ¼ d y ¼ d z ¼ d, and R d ¼ 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: Again, FðY is monotonically increasing (F 0 ðYÞ ! 0; 8Y), 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 Àd x 0 Àg x Àg yx Àd y Àg yz 0 Àg z Àd z 2 6 4 3 7 5 where q½R x gðÀCZ þ hÞ À d If we now assume that d x ¼ d y ¼ d z ¼ d, we have the same characteristic equation for the connectionist formulation as in the general case (see Equation 23):