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 wholenetwork 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.001Introduction
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; WestEberhard, 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 patternforming 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 (ShenOrr 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 contextdependence 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 coopted 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 traitspecific variational properties (Wagner and Altenberg, 1996; Wagner et al., 2007; Wagner, 2014). Finally, structural modularity allows for a finetuned response to specific selective pressures by minimizing offtarget 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 problemdependent (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 regulationexpression 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; PerezCarrasco et al., 2018; Page and PerezCarrasco, 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 simulationbased 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 SalazarCiudad et al., 2000; SalazarCiudad et al., 2001a; SalazarCiudad 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 realworld 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 anteroposterior (A–P) axis of the embryo by the protein products of the maternal coordinate 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 gapgene crossregulation, 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 coordinate and gap genes together then regulate pairrule and segmentpolarity genes, which form a molecular prepattern 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 simulationbased analyses identified a number of mechanisms driving gap gene expression. One of these consists of two doublenegative (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 AlvarezBuylla, 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 AlvarezBuylla, 2010). Their spatiotemporal 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; PerezCarrasco 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 coexpression modules (sometimes also called regulatory modules; see, for example, Eisen et al., 1998; BarJoseph 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 anticorrelated 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 nonlinearity 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 morphogendriven patternforming 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 (PanovskaGriffiths et al., 2013; PerezCarrasco et al., 2018). AC/DC subcircuits represent one of the simplest known genetic systems able to produce both switchlike (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 anteroposterior (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 spatiotemporal 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 doublenegative (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 dorsoventral patterning in the vertebrate neural tube (Balaskas et al., 2012; PanovskaGriffiths et al., 2013; PerezCarrasco 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 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 timevariable inputs from maternal gradients (Verd et al., 2017) and autoactivation. It has been shown previously that autoactivation 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 timedependent boundary conditions render these systems timevariant or nonautonomous, and cause their phase portraits to change their geometry over time (Verd et al., 2017). Timevariance or nonautonomy 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 timevariant, nonautonomous 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 HbKr 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 HbGt 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 KrKni 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) timevariant or nonautonomous 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 KrKni to the KniGt 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 subterminal 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 timevariant or nonautonomous 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 crosscorrelations 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 anteroposterior 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 counterintuitive 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íaSolache 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 HbKr 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íaSolache et al., 2010; Crombach et al., 2014). In the following subsections, we will provide an analysis that resolves this apparent paradox.
Intraspecific 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üssleinVolhard, 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).
Interspecific 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 gapgap crossregulatory 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 gapgap 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 gapgap crossrepressive 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 gapgap crossregulatory interactions specifically affect the extent of static versus shifting patterning regimes. In contrast, other aspects of gap gene patterning, such as the alternatingcushions 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 (switchlike) and oscillatory dynamics depending on parameter values and boundary conditions, amongst others (PanovskaGriffiths et al., 2013). In this paper, we show that each circuit is active in a particular region of influence along the anteroposterior 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 AlvarezBuylla, 2010; Benítez and AlvarezBuylla, 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ándezHerná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 subtype of dynamical modules, which exhibit characteristic and robust behaviour based on the selforganising reciprocal interaction of genetic and tissuelevel 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 timevariant or nonautonomous 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 steadystate analysis of a timeinvariant or autonomous version of the system. Transient dynamics, and the explicit dependence of regulatory structure on time implied by timevariance or nonautonomy, should be considered the default assumption. Steadystate dynamics and timeinvariant 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, switchlike, 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 SchmidtOtt, 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 gapgap 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íaSolache et al., 2010). Gap domain shifts are much more pronounced in this species compared to D. melanogaster or M. abdita (GarcíaSolache et al., 2010; Crombach et al., 2014). The beetle T. castaneum shows an even more dynamic mode of segment determination: only the anteriormost domains form simultaneously, while more posterior domains are generated sequentially by sustained oscillations in segmentation gene expression (Sarrazin et al., 2012; ElSherif et al., 2012).
An extended analysis of our AC/DC circuits (see Appendix 1) reveals that they can be induced to drive sustained limitcycle oscillations with relatively small additional changes in the values of parameters that determine crossregulatory interactions (summarized in Figure 6) (PanovskaGriffiths et al., 2013; PerezCarrasco et al., 2018; Page and PerezCarrasco, 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 pairrule genes (Sarrazin et al., 2012; ElSherif et al., 2012; ElSherif et al., 2014), which are known to be essential for segment determination in all arthropods. Pairrule genes in T. castaneum may also be regulated by AC/DClike regulatory subcircuits (Choe et al., 2006). These surprising resemblances suggest that ancestral gap and pairrule expression may have relied on similar regulatory principles, and that the regulatory changes required to turn sequential (shortgermband) into simultaneous (longgermband) 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 (Jaeger and Crombach, 2012). Significant advances in ‘omics’ and ‘bigdata’ techniques, genetic perturbation by gene knockdown, 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 spatiotemporal 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 largerscale continuous models, for instance genomescale 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; PanovskaGriffiths et al., 2013; PerezCarrasco et al., 2018). So can the growing number of experimentally verified Turingtype 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 coexpression 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 diffusionless gap gene circuit
What we refer to as the ‘full model’ in this paper corresponds to a diffusionless 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 onedimensional row of nuclei, arranged along the anteroposterior (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 transcriptionfactor proteins encoded by trunk gap genes hb, Kr, kni, and gt. ${g}_{i}^{a}(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:
where ${R}^{a}$ is the rate of protein production, and ${\lambda}^{a}$ the rate of protein decay. $\varphi $ is a sigmoid regulationexpression function which is used to represent the coarsegrained saturating kinetics of transcriptional regulation. It is defined as follows:
where
$G=\{\mathrm{\u210e\mathit{b}},\mathrm{\mathit{K}\mathit{r}},\mathrm{\mathit{k}\mathit{n}\mathit{i}},\mathrm{\mathit{g}\mathit{t}}\}$ denotes the set of trunk gap genes, and $M=\{\mathrm{Bcd},\mathrm{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 spatiotemporal 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}_{i}^{m}$.
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
Request a detailed protocolValues for parameters ${R}^{a}$, ${\lambda}^{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 rootmeansquare (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 spatiotemporal 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 switchlike 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 protocolWe 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 gapgap interactions become increasingly predominant during late C13 and C14A. Since our analysis focuses on gapgap crossregulation, 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:
$({\mathrm{Hb}}_{{T}_{i}}^{j},{\mathrm{Kr}}_{{T}_{i}}^{j},{\mathrm{Kni}}_{{T}_{i}}^{j},{\mathrm{Gt}}_{{T}_{i}}^{j})$ 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
Request a detailed protocolWe 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 timevariable 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 C14T1 as initial conditions. Simulations of AC/DC subcircuits include C14A only (no mitosis or division. During C14A, zygotic gap gene crossregulation 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 protocolSince they implement timevariable maternal inputs, AC/DC subcircuits are timevariant 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 NewtonRaphson 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 wellstudied behaviour of the threegene repressilator circuit (Elowitz and Leibler, 2000; Müller et al., 2006; GarciaOjalvo et al., 2004; Strelkowa and Barahona, 2010), and then extend the analysis to the AC/DC network (PanovskaGriffiths 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 threegene 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 sigmoidshaped regulationexpression function conforms to these general conditions.
We use the following abstract formulation of the repressilator circuit, where G stands for a regulationexpression function conforming to the conditions prescribed above, and d are gene product degradation rates:
This system has a unique steady state. Defining ${g}_{i}\equiv {\displaystyle \frac{1}{{d}_{i}}}{G}_{i}$, for $i=\{X,Y,Z\}$, 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:
Due to the fact that the three regulationexpression 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:
where
The ${\gamma}_{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 state (so, for example ${\gamma}_{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 $\lambda $ associated with the linearised system at the steady state. These are the roots of the characteristic polynomial given by
where we define $\mathrm{\Gamma}$ to be a negative real number.
For simplicity, we consider the special case where all degradation rates are equal, ${d}_{x}={d}_{y}={d}_{z}=d$:
The three distinct complex roots of ${(\lambda +d)}^{3}={\mathrm{\Gamma}}^{3}$ are given by (see also Appendix 1—figure 1C):
It follows that the three eigenvalues associated with the linearised system are
${\lambda}_{1}$ is always negative and real. However, the real parts of ${\lambda}_{2}$ and ${\lambda}_{3}$ can take either sign, since
Proof:
This leaves us with two alternative possibilities. From the equations in Equation (17) we know that $Re({\lambda}_{2,3})>0$ if $\frac{\mathrm{\Gamma}}{2}d>0$. From this, it follows that $Re({\lambda}_{2,3})>0$ if $\mathrm{\Gamma}>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, (MalletParet and Smith, 1990) have shown that, in monotone cyclic feedback systems, the $\omega $limit set of any bounded orbit can be embedded in ${\mathbb{R}}^{2}$. Under these conditions, which are met in our case, and in this twodimensional 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 limitcycle oscillations.
There is a second possibility. If $\mathrm{\Gamma}<2d\u27f9{R}_{e}({\lambda}_{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({\lambda}_{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 limitcycle oscillations, or damped oscillations (in agreement with Page and PerezCarrasco, 2018).
Analysis of the AC/DC circuit
We obtain an AC/DC circuit (PanovskaGriffiths 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:
Regulationexpression 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 regulationexpression 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:
where
and
${\gamma}_{yx}$ measures the sensitivity of the repressive interaction of X on Y, and ${\gamma}_{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:
From this, it follows that
If we let $d+\lambda =\eta $ and ${\stackrel{~}{\mathrm{\Gamma}}}^{3}={\gamma}_{x}{\gamma}_{yx}{\gamma}_{z}>0$, then Equation (25) becomes
Now, let $x={\gamma}_{yz}{\gamma}_{z}>0$, and Equation (26) becomes
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 $\stackrel{~}{\mathrm{\Gamma}}$ are evaluated at the steady state. They remain positive since
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 $\stackrel{~}{\mathrm{\Gamma}}$ 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 $\eta $, 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 ${\eta}_{i},i\in \{1,2,3\}$, where ${\eta}_{i}={\lambda}_{i}+{d}_{i}$ and ${\lambda}_{i},i\in \{1,2,3\}$, are the eigenvalues associated with the corresponding steady state of the AC/DC circuit. If we consider ${d}_{i}=d$, $\forall 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), 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).
Damped oscillations are observed if the underlying phase portrait has a unique stable spiral sink (Appendix 1—figure 4, parameter values used for simulation are shown in Appendix 1—table 2).
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 saddle_{1,2} (Appendix 1—figure 5, parameter values used for simulation are shown in Appendix 1—table 3).
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 coexist with a saddle and a point attractor.
A stable spiral sink can coexist with a saddle_{2,1} and a point attractor in a bistable regime. The saddle probably needs to have a 2dimensional 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 (PerezCarrasco et al., 2018; Page and PerezCarrasco, 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:
where:
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 ${\delta}_{+}$), and the concentration of X. $R$ represents the maximum rate of activation for every gene. More generally, we write $u$ as:
based on which Equation (30) becomes:
In what follows, we will assume that ${R}_{x}={R}_{y}={R}_{z}=R$, ${d}_{x}={d}_{y}={d}_{z}=d$, and $\frac{R}{d}=\overline{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}^{\prime}(Y)\ge 0,\forall 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
where
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):
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
Form and function in gene regulatory networks: the structure of network motifs determines fundamental properties of their dynamical state spaceJournal of the Royal Society Interface 13:20160179.https://doi.org/10.1098/rsif.2016.0179
 2

3
An Introduction to Systems Biology: Design Principles of Biological CircuitsBoca Raton: Chapman and Hall/CRC.

4
Network motifs: theory and experimental approachesNature Reviews Genetics 8:450–461.https://doi.org/10.1038/nrg2102

5
Gene circuit analysis of the terminal gap gene huckebeinPLOS Computational Biology 5:e1000548.https://doi.org/10.1371/journal.pcbi.1000548

6
Structure and evolution of transcriptional regulatory networksCurrent Opinion in Structural Biology 14:283–291.https://doi.org/10.1016/j.sbi.2004.05.004
 7

8
Computational discovery of gene modules and regulatory networksNature Biotechnology 21:1337–1342.https://doi.org/10.1038/nbt890

9
Mechanism and biological explanation*Philosophy of Science 78:533–557.https://doi.org/10.1086/661513

10
Explanation: a mechanist alternativeStudies 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
Dynamic mechanistic explanation: computational modeling of circadian rhythms as an exemplar for cognitive scienceStudies in History and Philosophy of Science Part A 41:321–333.https://doi.org/10.1016/j.shpsa.2010.07.003
 12
 13

14
Epidermal patterning in arabidopsis: models make a differenceJournal of Experimental Zoology Part B: Molecular and Developmental Evolution 316B:241–253.https://doi.org/10.1002/jez.b.21398
 15

16
The Evolution of Complexity by Means of Natural SelectionPrinceton: Princeton University Press.
 17
 18
 19

20
Lineage explanations: explaining how biological mechanisms changeThe British Journal for the Philosophy of Science 60:51–78.https://doi.org/10.1093/bjps/axn047
 21
 22
 23

24
Parallel simulated annealing by mixing of statesJournal of Computational Physics 148:646–662.https://doi.org/10.1006/jcph.1998.6134
 25

26
Complexity: against systemsTheory in Biosciences 130:229–245.https://doi.org/10.1007/s1206401101214
 27

28
Efficient reverseengineering of a developmental gene regulatory networkPLOS Computational Biology 8:e1002589.https://doi.org/10.1371/journal.pcbi.1002589
 29

30
Gap gene regulatory dynamics evolve along a genotype networkMolecular Biology and Evolution 33:1293–1307.https://doi.org/10.1093/molbev/msw013

31
Evolution of evolvability in gene regulatory networksPLOS Computational Biology 4:e1000112.https://doi.org/10.1371/journal.pcbi.1000112
 32
 33

34
The evolution of evolvabilityIn: C Langton, editors. Artificial Life: The Proceedings of an Interdisciplinary Workshop on the Synthesis and Simulation of Living Systems. Redwood City: AddisonWesley. pp. 201–220.
 35
 36
 37
 38
 39

40
The evolution of hierarchical gene regulatory networksNature Reviews Genetics 10:141–148.https://doi.org/10.1038/nrg2499

41
Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesisJournal of Cell Science 61:31–70.

42
Community detection in graphsPhysics Reports 486:75–174.https://doi.org/10.1016/j.physrep.2009.11.002
 43
 44

45
A systematic analysis of the gap gene system in the moth midge clogmia albipunctataDevelopmental Biology 344:306–318.https://doi.org/10.1016/j.ydbio.2010.04.019
 46
 47
 48
 49
 50

51
Evolvability as the proper focus of evolutionary developmental biologyEvolution & Development 9:393–401.https://doi.org/10.1111/j.1525142X.2007.00176.x

52
Dynamical patterning modules in plant development and evolutionThe International Journal of Developmental Biology 56:661–674.https://doi.org/10.1387/ijdb.120027mb
 53
 54

55
Differential Equations, Dynamical Systems, and an Introduction to Chaos (Third Edition)Amsterdam, Netherlands: Elsevier.
 56
 57
 58
 59

60
Mutual interaction in network motifs robustly sharpens gene expression in developmental processesJournal of Theoretical Biology 252:131–144.https://doi.org/10.1016/j.jtbi.2008.01.027
 61
 62

63
Known maternal gradients are not sufficient for the establishment of gap domains in Drosophila melanogasterMechanisms of Development 124:108–128.https://doi.org/10.1016/j.mod.2006.11.001

64
The gap gene networkCellular and Molecular Life Sciences 68:243–274.https://doi.org/10.1007/s000180100536y

65
The inheritance of process: a dynamical systems approachJournal of Experimental Zoology Part B: Molecular and Developmental Evolution 318:591–612.https://doi.org/10.1002/jez.b.22468

66
Life’s AttractorsIn: O. S Soyer, editors. Evolutionary Systems Biology, 751. New York: Springer. pp. 93–119.https://doi.org/10.1007/9781461435679_5

67
Bioattractors: dynamical systems theory and the evolution of regulatory processesThe Journal of Physiology 592:2267–2281.https://doi.org/10.1113/jphysiol.2014.272385

68
Evolutionary Systems Biology 2.0Dynamical modularity of the genotypephenotype map, Evolutionary Systems Biology 2.0, Berlin, Springer.

69
On the concept of mechanism in developmentIn: 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
A spectrum of modularity in multifunctional gene circuitsMolecular Systems Biology 13:925.https://doi.org/10.15252/msb.20167347

71
The Origins of Order: SelfOrganization and Selection in EvolutionOxford: Oxford University Press.
 72
 73
 74

75
Cell signaling as a cognitive processThe EMBO Journal 36:568–582.https://doi.org/10.15252/embj.201695383
 76
 77

78
Bicoid occurrence and Bicoiddependent hunchback regulation in lower cyclorrhaphan fliesEvolution & Development 10:413–420.https://doi.org/10.1111/j.1525142X.2008.00252.x

79
Maternal activation of gap genes in the hover fly EpisyrphusDevelopment 137:1709–1719.https://doi.org/10.1242/dev.046649
 80
 81
 82
 83
 84

85
The PoincareBendixson theorem for monotone cyclic feedback systemsJournal of Dynamics and Differential Equations 2:367–421.https://doi.org/10.1007/BF01054041
 86

87
Canalization of gene expression and domain shifts in the Drosophila blastoderm by dynamical attractorsPLOS Computational Biology 5:e1000303.https://doi.org/10.1371/journal.pcbi.1000303
 88

89
Reusable building blocks in biological systemsJournal of the Royal Society Interface 15:20180595.https://doi.org/10.1098/rsif.2018.0595

90
A connectionist model of developmentJournal of Theoretical Biology 152:429–453.https://doi.org/10.1016/S00225193(05)803911
 91

92
Mutants highlight the modular control of butterfly eyespot patternsEvolution and Development 5:180–187.https://doi.org/10.1046/j.1525142X.2003.03029.x

93
Comparative insights into questions of lepidopteran wing pattern homologyBMC Developmental Biology 6:52.https://doi.org/10.1186/1471213X652
 94

95
A generalized model of the repressilatorJournal of Mathematical Biology 53:905–937.https://doi.org/10.1007/s0028500600359
 96
 97

98
Dynamical patterning modules: a "pattern language" for development and evolution of multicellular formThe International Journal of Developmental Biology 53:693–705.https://doi.org/10.1387/ijdb.072481sn
 99

100
Gene regulatory network controlling embryonic specification in the sea urchinCurrent Opinion in Genetics & Development 14:351–360.https://doi.org/10.1016/j.gde.2004.06.004
 101

102
The fintolimb transition as the reorganization of a turing patternNature Communications 7:11582.https://doi.org/10.1038/ncomms11582

103
Degradation rate uniformity determines success of oscillations in repressive feedback regulatory networksJournal of the Royal Society Interface 15:20180157.https://doi.org/10.1098/rsif.2018.0157
 104

105
A gene regulatory motif that generates oscillatory or multiway switch outputsJournal of the Royal Society Interface 10:20120826.https://doi.org/10.1098/rsif.2012.0826
 106

107
Genetic variation in pleiotropy: differential epistasis as a source of variation in the allometric relationship between long bone lengths and body weightEvolution; International Journal of Organic Evolution 62:199–213.https://doi.org/10.1111/j.15585646.2007.00255.x

108
Function does not follow form in gene regulatory circuitsScientific Reports 5:13015.https://doi.org/10.1038/srep13015
 109

110
Reverse engineering the gap gene network of Drosophila melanogasterPLOS Computational Biology 2:e51.https://doi.org/10.1371/journal.pcbi.0020051
 111
 112

113
FlyEx, the quantitative atlas on segmentation gene expression at cellular resolutionNucleic Acids Research 37:D560–D566.https://doi.org/10.1093/nar/gkn717
 114

115
The Shape of Life: Genes, Development, and the Evolution of Animal FormChicago, IL, United States: University of Chicago Press.
 116

117
Mechanism of eve stripe formationMechanisms of Development 49:133–158.https://doi.org/10.1016/09254773(94)00310J
 118

119
Segmentation gene expression in the mothmidge clogmia albipunctata (Diptera, psychodidae) and other primitive dipteransDevelopment Genes and Evolution 209:145–154.https://doi.org/10.1007/s004270050238

120
Gene networks capable of pattern formation: from induction to reactiondiffusionJournal of Theoretical Biology 205:587–603.https://doi.org/10.1006/jtbi.2000.2092
 121
 122
 123
 124
 125
 126
 127
 128
 129

130
The architecture of complexityProceedings of the American Philosophical Society 106:467–482.

131
LargeScale metabolic models: from reconstruction to differential equationsIndustrial Biotechnology 9:179–184.https://doi.org/10.1089/ind.2013.0003

132
Constructing kinetic models of metabolism at genomescales: a reviewBiotechnology Journal 10:1345–1359.https://doi.org/10.1002/biot.201400522
 133
 134

135
Switchable genetic oscillator operating in quasistable modeJournal of the Royal Society Interface 7:1071–1082.https://doi.org/10.1098/rsif.2009.0487
 136
 137
 138

139
Characterization of the Drosophila segment determination morphomeDevelopmental Biology 313:844–862.https://doi.org/10.1016/j.ydbio.2007.10.037
 140
 141

142
Gene cooption in physiological and morphological evolutionAnnual Review of Cell and Developmental Biology 18:53–80.https://doi.org/10.1146/annurev.cellbio.18.020402.140619
 143

144
Dynamic maternal gradients control timing and ShiftRates for Drosophila Gap Gene ExpressionPLOS Computational Biology 13:e1005285.https://doi.org/10.1371/journal.pcbi.1005285
 145
 146

147
Modularity in animal development and evolution: elements of a conceptual framework for EvoDevoThe Journal of Experimental Zoology 285:307–325.https://doi.org/10.1002/(SICI)1097010X(19991215)285:4<307::AIDJEZ2>3.0.CO;2V
 148
 149
 150
 151

152
The pleiotropic structure of the genotypephenotype map: the evolvability of complex organismsNature Reviews Genetics 12:204–213.https://doi.org/10.1038/nrg2949

153
Multiple functions of a feedforwardloop gene circuitJournal of Molecular Biology 349:501–514.https://doi.org/10.1016/j.jmb.2005.04.022
 154
 155
 156

157
Enhancer responses to similarly distributed antagonistic gradients in developmentPLOS Computational Biology 3:e84.https://doi.org/10.1371/journal.pcbi.0030084
Decision letter

Lee AltenbergReviewing Editor; The KLI Institute, United States

Patricia J WittkoppSenior 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 anteroposterior 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 anteroposterior 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 nonautonomous 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 evodevo 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 nonautonomous 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 evodevo 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 nonmodularity 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 yesno 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 nonautonomous 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 nonautonomous dynamical systems in other networks that are not as well studied.
https://doi.org/10.7554/eLife.42832.022Author 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 nonmodularity 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 yesno 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 nonmodularity 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 nonautonomous 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 timevariance or nonautonomy more gently. Since “nonautonomy” is a potentially confusing term, we now use “timevariant” along with it throughout the manuscript. This both clarifies the concept and maintains continuity with our earlier publications on nonautonomous 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 nonautonomous 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.023Article and author information
Author details
Funding
MINECO (BFU200910184/BFU201233775/SEV20120208)
 Johannes Jaeger
European Commission (FP7KBBE20115/289434 (BioPreDyn))
 Johannes Jaeger
MECEMBL
 Johannes Jaeger
La Caixa Savings Bank
 Berta Verd
KLI Klosterneuburg
 Berta Verd
Wissenschaftskolleg zu Berlin
 Johannes Jaeger
MaxPlanckGesellschaft
 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érezCarrasco 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
 Patricia J Wittkopp, University of Michigan, United States
Reviewing Editor
 Lee Altenberg, The KLI Institute, United States
Publication history
 Received: January 14, 2019
 Accepted: June 5, 2019
 Accepted Manuscript published: June 6, 2019 (version 1)
 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

 3,691
 Page views

 687
 Downloads

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