Dynamic Maternal Gradients Control Timing and Shift-Rates for Drosophila Gap Gene Expression

Pattern formation during development is a highly dynamic process. In spite of this, few experimental and modelling approaches take into account the explicit time-dependence of the rules governing regulatory systems. We address this problem by studying dynamic morphogen interpretation by the gap gene network in Drosophila melanogaster. Gap genes are involved in segment determination during early embryogenesis. They are activated by maternal morphogen gradients encoded by bicoid (bcd) and caudal (cad). These gradients decay at the same time-scale as the establishment of the antero-posterior gap gene pattern. We use a reverse-engineering approach, based on data-driven regulatory models called gene circuits, to isolate and characterise the explicitly time-dependent effects of changing morphogen concentrations on gap gene regulation. To achieve this, we simulate the system in the presence and absence of dynamic gradient decay. Comparison between these simulations reveals that maternal morphogen decay controls the timing and limits the rate of gap gene expression. In the anterior of the embyro, it affects peak expression and leads to the establishment of smooth spatial boundaries between gap domains. In the posterior of the embryo, it causes a progressive slow-down in the rate of gap domain shifts, which is necessary to correctly position domain boundaries and to stabilise the spatial gap gene expression pattern. We use a newly developed method for the analysis of transient dynamics in non-autonomous (time-variable) systems to understand the regulatory causes of these effects. By providing a rigorous mechanistic explanation for the role of maternal gradient decay in gap gene regulation, our study demonstrates that such analyses are feasible and reveal important aspects of dynamic gene regulation which would have been missed by a traditional steady-state approach. More generally, it highlights the importance of transient dynamics for understanding complex regulatory processes in development.

Biological systems depend on time. Like everything else that persists for more than an 2 instant, there is a temporal dimension to their existence. This much is obvious. What is 3 less obvious, however, is the active role that time plays in altering the rules governing 4 biological processes. For instance, fluctuating environmental conditions modify the 5 selective pressures that drive adaptive evolutionary change [1,[3][4][5], time-dependent 6 inductive signals or environmental cues trigger and remodel developmental 7 pathways [6,7], and dynamic morphogen gradients influence patterning, not only across 8 space but also through time [8][9][10][11][12][13][14][15][16]. In spite of this, many current attempts at 9 understanding biological processes neglect important aspects of this temporal 10 dimension [17]. For practical reasons, experimental studies often glance over the 11 detailed dynamics of a process, and focus on its end product or output pattern instead. 12 Similarly, modelling studies frequently restrict themselves to a small-enough time 13 window allowing them to ignore temporal changes in the rules governing the system. 14 Accuracy is sacrificed and the scope of the investigation limited for the sake of 15 simplicity and tractability. Although reasonable, and often even necessary, such 16 simplifications can lead us to miss important aspects of biological regulatory dynamics. 17 We set out to tackle explicitly time-dependent aspects of morphogen interpretation The gap gene system is one of the most thoroughly studied developmental gene 30 regulatory networks today. For our particular purposes, we take advantage of the fact 31 that it has been extensively reverse-engineered using data-driven modelling. This 32 approach is based on fitting dynamical models of gap gene regulation, called gene 33 circuits, to quantitative spatio-temporal gene expression data [19-27, 29, 34, 35]. 34 Dynamical models capture how a given regulatory process unfolds over time. They 35 are frequently formulated in terms of ordinary differential equations (ODEs) with 36 parameter values that remain constant over time. Such equations represent an 37 autonomous dynamical system. Central to the analysis of such dynamical systems is the 38 concept of phase space and its associated features (S1 Figure A). Phase (or state) space 39 is an abstract space that contains all possible states of a system. Its axes are defined by 40 the state variables, which in our case represent the concentrations of transcription 41 factors encoded by the gap genes. Trajectories through phase space describe how a 42 system's state changes as time progresses. The trajectories of a gap gene circuit 43 describe how transcription factor concentrations change over time. All trajectories taken 44 together constitute the flow of the system. This flow is shaped by the regulatory 45 structure of the underlying network-the type (activation/repression) and strength of 46 interactions between the constituent factors-which is given by the system's parameters. 47 Since these parameters are constant over time in an autonomous system, the 48 trajectories are fully determined given a specific set of initial conditions. Once the 49 system's variables no longer change, it has reached a steady state. Steady states can be 50 stable-such as attractors with converging trajectories from all directions defining a 51 basin of attraction-or unstable-such as saddles; where trajectories converge only 52 along certain directions and diverge along others. The type and arrangement of steady 53 states, and their associated basins of attraction define the phase portrait of the system 54 (S1 Figure A). There exist powerful analytical tools to analyse and understand the 55 phase portrait and the range of dynamic behaviours determined by it. Geometrical 56 analysis of the phase portrait enables us to build up a rich qualitative understanding of 57 the dynamics of non-linear autonomous systems without solving the underlying 58 equations analytically [36]. 59

3/30
The application of dynamical systems concepts and phase space analysis to the 60 study of cellular and developmental processes has a long history (see [37][38][39] for recent 61 reviews). In particular, it has been successfully applied to the study of the gap gene 62 system. Manu and colleagues [22,23,40] examined the dynamics and robustness of gap 63 gene regulation in D. melanogaster using diffusion-less gene circuits fit to quantitative 64 expression data. These models have a four-dimensional phase space, where the axes 65 represent the concentrations of transcription factors encoded by the trunk gap genes hb, 66 Kr, gt, and kni. The analysis of these phase portraits yields a rigorous understanding of 67 the patterning capabilities of the system.

68
The analysis by Manu et al. [23] corroborated and expanded upon earlier genetic 69 evidence [41] indicating that the regulatory dynamics responsible for domain boundary 70 placement in the anterior versus the posterior of the embryo are very different. In the 71 anterior, spatial boundaries of gap gene expression domains are positioned statically, 72 meaning that they remain in place over time [42]. Stationary boundaries are regulated 73 in two distinct ways [23]. basins of attraction (attractor selection) ( Fig. 2A). In both of these cases, patterning is 78 largely governed by the position of attractors in a multi-stable phase space.

79
In contrast, gap domain boundaries in the posterior of the embryo shift anteriorly 80 over time [25,42]. In this region, the system always remains far from steady state, and 81 the dynamics of gene expression are transient. Therefore, trajectories here are fairly 82 independent of precise attractor positions. The model by Manu et al. [23] shows that 83 posterior gap gene expression is governed by an unstable manifold ( Fig. 2A). An 84 unstable manifold is the trajectory connecting a saddle to an attractor (S1 Figure A). 85 The authors demonstrate that this manifold has canalising properties since it 86 compresses many incoming neighbouring trajectories into an increasingly smaller 87 sub-volume of phase space over time [23]. This explains the observed robustness of

93
Despite its explanatory power, the analysis by Manu et al. [23] is limited in an 94 important way. In order to simplify phase space analysis, the authors implement 95 simplified dynamics of maternal morphogens Bcd and Cad in their model ( Fig. 2A).

96
They use a time-invariant exponential approximation to simulate the Bcd gradient and 97 Cad is assumed to reach a steady-state profile about 20-30 minutes before 98 gastrulation [22,23]. This steady-state profile is used for model analysis. (Based on this, 99 we will refer to this formulation as the static-Bcd gene circuit model in what follows).

100
Although reasonable, these simplifications affect the accuracy of the model, since Bcd 101 and Cad have their own expression dynamics on a similar time scale as gap proteins.

102
The Bcd gradient decays and Cad clears from much of the posterior trunk region 103 towards the end of the blastoderm stage ( Fig. 2B) [42]. This means that the 104 autonomous analysis of the static-Bcd model is not well suited to investigate the 105 dynamic interpretation of morphogen gradients. In particular, assuming autonomy 106 makes it impossible to isolate and study the explicitly time-dependent effects of 107 changing gradient concentrations on gap gene regulation and pattern formation.

108
For this reason, we consider the dynamics of maternal morphogens explicitly in our 109 model. We have obtained gap gene circuits that incorporate realistic time-variable 110 maternal gradients of Bcd and Cad (Fig. 2B) [26]. These gradients are implemented as 111 4/30 external inputs to gap gene regulation (see Models and Methods section). They are not 112 influenced by any of the state variables and, thus, are parameters of the system. This 113 means that our gap gene circuits become fully non-autonomous [54], since certain 114 parameter values now change over time. While non-autonomous equations are not 115 significantly more difficult to formulate or simulate than autonomous ones, phase space 116 analysis is far from trivial. As model parameters change, so does the geometry of the 117 phase portrait, and consequently system trajectories are actively shaped by this 118 time-dependence. Separatrices and attractors can change their position (geometrical 119 change), and steady states can be created and annihilated through bifurcation events 120 (topological change) (S1 Figure B). In autonomous systems, bifurcations can only occur 121 along the spatial axis of the model. In non-autonomous systems, they also occur in time, 122 implying that trajectories can switch from one basin of attraction to another during a 123 simulation run. We can think of time-variable phase portraits as embedded in parameter 124 space. We call the combination of phase and parameter space the configuration space of 125 the system. The configuration space on non-autonomous models hence encodes a much 126 richer repertoire of dynamical mechanisms of pattern formation than autonomous phase 127 space alone. This can complicate analysis and interpretation of the system considerably. 128  [23]. An exponential function fit to the Bcd profile at cycle C13 was used to calculate trajectories and phase portraits. All Cad profiles until time class T6 were considered for simulating trajectories, but phase portraits were calculated using the profile at T6 only. This gene circuit displays the following mechanisms for boundary formation: patterning between 35-51% A-P position takes place in a multi-stable regime close to steady state. The Gt boundary is established as the relevant attractor moves from high to low Gt concentration in more posterior nuclei. The Hb-Kr interface forms as the maternal Hb gradient places more anterior nuclei in the basin of an attractor with high Hb concentration, and more posterior nuclei in the basin of an attractor at high Kr concentration. Between 51 and 53% A-P position a saddle-node bifurcation takes place, and the dynamics become transient in nuclei posterior of 52%. These nuclei are all in the basin of the same attractor and approach it by first converging towards an unstable manifold. Anterior shifts in these posterior gap gene domains emerge from a coordinated succession of trajectories in more posterior nuclei approaching the unstable manifold. See [23] for details. (B) In the non-autonomous gap gene circuit analysed here, Bcd and Cad gradient profiles are included for every time point. They are used to calculate trajectories of the system and instantaneous phase portraits as discussed in the main text. Plots in (A) and (B) show % A-P position along the x-axes, and protein concentrations (in arbitrary units, au) along the y-axes as in Fig.1B). Using a simple model of a genetic toggle switch, we have established a methodology 129 for the characterisation of transient dynamics in non-autonomous systems (S1 130 Figure B), based on the analysis of instantaneous phase portraits [43,45]. Such portraits 131 are generated by fixing the values of system parameters starting at a given point in 132 time, and then determining the geometrical arrangement of saddles, attractors, and 133 their basins under these "frozen" conditions. The overall non-autonomous trajectory of 134 the system is given by a series of instantaneous phase portraits over time. With 135 sufficiently high temporal resolution, this method yields an accurate picture of the 136 non-autonomous mechanisms of pattern formation implemented by the system. These 137 mechanisms can be classified into four broad categories [43]: (1) transitions of the 138 system from one steady state to another, (2) pursuit of a moving attractor within a 139 basin of attraction, (3) geometrical capture, where a trajectory crosses a separatrix, and 140 (4) topological capture, where a trajectory suddenly falls into a new basin of attraction 141 due to a preceding bifurcation event (S1 Figure B). This classification scheme can be 142 used to characterise the dynamical repertoire of non-autonomous models in a way 143 analogous to phase space analysis in autonomous dynamical systems.

144
In this paper, we present a detailed analysis of a non-autonomous gap gene circuit. 145 Specifically, we use the model to address the effect of non-autonomy, i. e. the effect of 146 time-variable maternal gradient concentrations, on gap gene regulation (Fig. 2). To   [21], modified to include time-variable external regulatory 158 inputs as previously described [26,34]. Gene circuits are hybrid models with discrete cell 159 divisions and continuous gene regulatory dynamics. The basic objects of the model  The state variables of the system consist of the concentration levels of proteins 167 produced by the trunk gap genes hb, Kr, gt, and kni. We denote the concentration of 168 gap protein a in nucleus i at time t by g a i (t). Change in protein concentration over time 169 is given by the following set of ODEs: where R a , D a and λ a are rates of protein production, diffusion, and decay, respectively. 171 Diffusion depends on the distance between neighbouring nuclei, which halves at nuclear 172 division; thus, D a depends on the number of preceding divisions n. φ is a sigmoid regulation. It is defined as follows: parameters implies non-autonomy of the dynamical system (see Introduction and [54]). 182 Interconnectivity matrices W and E define interactions among gap genes, as well as 183 regulatory inputs from external inputs, respectively. The elements of these matrices, 184 w ba and e ma , are called regulatory weights. They encode the effect of regulator b or m 185 on target gene a. These weights may be positive (representing an activating regulatory 186 input), negative (representing repression), or near zero (representing the absence of a 187 regulatory interaction). h a is a threshold parameter that represents the activation state 188 of target gene a in the absence of any spatially and temporally specific regulatory input. 189 This term incorporates the regualtory influence of factors that are not expressed in a 190 spatially specific manner (for example, the pioneer factor Zelda [31]). Equation (1) 191 determines regulatory dynamics during interphase. In order to accurately implement 192 the non-instantaneous duration of the nuclear division between C13 and C14A, the 193 production rate R a is set to zero during a mitotic phase, which immediately precedes 194 the instantaneous nuclear division. Mitotic schedule as in [26]. 195 Model fitting and selection 196 We determine the values for parameters R a , λ a , W , E, and h a using a 197 reverse-engineering approach [19,25,26,34]. For this purpose, we numerically solve gene 198 circuit equations (1) across the region between 35 and 92% A-P position using a 199 Runge-Kutta Cash-Karp adaptive step-size solver [26]. Models are fit to a previously 200 published quantitative data set of spatio-temporal gap protein expression [26,42,46]  Model fitting was performed using a global optimization algorithm called parallel Lam 203 Simulated Annealing (pLSA) [47]. We use a weighted least squares cost function as 204 previously described [26].

205
To enable comparison of our results to the static-Bcd gene circuit analysis by Manu 206 et al. [23], we keep model formalism and fitting procedure as similar as possible to this 207 earlier study. Manu and colleagues fitted gene circuits including a diffusion term, but 208 analysed the model with diffusion rates D a set to zero [23]. This diffusion-less approach 209 reduces the phase space of the model from hundreds of dimensions to 4 by spatially 210 uncoupling the equations and considering each nucleus independently from its 211 neighbours. Dimensionality reduction is essential for geometrical analysis of phase space. 212 Unfortunately, setting diffusion to zero in our best 3 (of a total of 100) non-autonomous 213 gene circuits fitted to data with non-zero diffusion terms leads to severe patterning 214 defects (see S2 Figure for common patterning defects). This is likely due to numerical, 215 not biological issues, since we do find circuit solutions that correctly reproduce gap gene 216 8/30 patterns both in the presence and absence of diffusion using an alternative fitting 217 approach: that fixes diffusion parameters D a to zero during optimization (see below).

218
To further facilitate comparison with the static-Bcd model, we constrained the signs of 219 regulatory weights to those reported in Manu et al. [23]. In previoius work, we have 220 verified this network structure extensively against experimental data [18,25,26,34]. 221 Optimization was performed on the Mare Nostrum supercomputer at the Barcelona 222 Supercomputing Centre (http://www.bsc.es). One optimization run took approximately 223 35 min on 64 cores.

224
The purpose of our reverse-engineering approach is not to sample parameter space 225 systematically, but instead to discover whether there are specific model-fitting solutions 226 that are consistent with the biological evidence and reproduce the dynamics of gap gene 227 expression correctly. Global optimization algorithms are stochastic heuristics without 228 guaranteed convergence, which means that for complex non-linear problems many 229 optimization runs will fail or end up at sub-optimal solutions (see also discussions 230 in [24,26,33]). In order to find the best-fitting solution, we therefore select solutions 231 from 200 initial fitting runs as follows: (1) we discard numerically unstable circuits;

232
(2) we only consider solutions with a root-mean-square (RMS) score less than 20.0 as 233 most circuits with scores above this threshold show gross patterning defects; (3) we use 234 visual inspection to detect remaining gross patterning defects among selected circuits 235 (missing or bimodal domains, and disconnected boundaries. See S2 Figure) as previously 236 described [34]. Out of the resulting 7 highest scoring circuits, only 3 recover the shifting 237 dynamics of posterior gap domains. In order to rule out diffusion as a 238 pattern-generating mechanism in these circuits, we compared their performance in the 239 presence and absence of diffusion (see above). For this purpose, we used values of 240 diffusion rates D a obtained by fitting our non-autonomous models with diffusion. All 241 three circuits produce satisfactory gap gene patterns (including anteriorly shifting 242 posterior trunk domains) whether diffusion is present or not. The best fit among these 243 was selected for detailed analysis (see S1 Table, for parameter values).

244
The residual error of our best-fitting diffusion-less circuit (RMS = 10.73) lies at the 245 lower end of the range of residual errors for fully-non-autonomous circuits with diffusion, 246 which range from RMS scores of 10.43 to 13.32 [26]. This lends further support to the 247 notion that diffusion is not essential for gap gene patterning. Moreover, our previous 248 work also shows that circuits which were fit without weighting the data show somewhat 249 lower RMS scores of 8.71 to 10.11 despite exhibiting more patterning defects at late 250 stages [26]. The RMS score of the static-Bcd model (fit without weights) is higher, at 251 10.76 [22]. Taken together, this implies a slightly better quality-of-fit of our fully 252 non-autonomous diffusion-less model compared to the static-Bcd diffusion-less circuits 253 of Manu et al. [22].

254
Gap gene circuit analysis 255 We characterise the time-variable geometry and topology of phase space in our fully 256 non-autonomous gap gene circuit for every nucleus in a sub-range of the fitted model 257 between 35 and 71% A-P position. This restricted spatial range allows us to simplify 258 the analysis by excluding the influence of terminal gap genes tll and hkb on patterning 259 (similar to the approach in [22]). We aim to identify those features of configuration 260 space that govern the placement of domain boundaries, and thus the patterning 261 capability of the gap gene system. We achieve this by generating instantaneous phase 262 portraits for the model [43,45]  Newton-Raphson method [48,49] as implemented by Manu et al. [23]. We classify 268 steady states according to their stability, which is determined by the corresponding 269 eigenvalues (see S1 Figure A). 270 Nuclei express a maximum of three trunk gap genes over developmental time, and 271 only two at any given time point. Therefore, we project four-dimensional phase 272 portraits into lower-dimensional representations to visualise them more easily. This 273 yields a graphical time-series of instantaneous phase portraits for each nucleus, which 274 allow us to track the movement, creation, and annihilation of steady states (typically 275 attractors and saddles) by bifurcations. The transient geometry of phase space governs 276 the non-autonomous trajectories of the system. We classify the dynamic behaviours 277 exhibited by these trajectories into transitions, pursuits, and captures according to our 278 previously established methodology (see Introduction and S1 Figure B) [43]. non-autonomous gene circuit through instantaneous phase portraits (S1 Figure B) [43], 294 analogous to the autonomous phase-space analysis presented by Manu and 295 colleagues [23] (Fig. 2). This type of analysis requires diffusion-less gap gene circuits to 296 keep the dimensionality of phase space at a manageable level. We obtained fully non-autonomous gap gene circuits that lack diffusion through 298 model fitting with diffusion parameters D a fixed to zero and interaction signs 299 constrained to those of previous works (as described in "Models and Methods"). This 300 resulted in a set of three selected, well-fitted circuits. The network topology of these 301 gene circuit models correspond to that shown in Fig. 3A. The following analysis is based 302 on the best-fitting model with a root mean square (RMS) residual error of 10.73, which 303 constitutes a slight overall improvement in quality-of-fit compared to static-Bcd models 304 (see "Models and Methods" and [22,26]). Its regulatory parameter values are listed in 305 S1 Table. 306 This diffusion-less non-autonomous gene circuit accurately reproduces gap gene 307 expression (Fig. 3B). In particular, it exhibits correct timing and relative positioning of 308 domain boundaries. Together with the fact that it fits the data equally well as 309 equivalent circuits with diffusion (see "Models and Methods", and [26]), this confirms 310 earlier indications that gap gene product diffusion is not essential for pattern formation 311 by the gap gene system [23,25]. Interestingly, previously published diffusion-less 312 static-Bcd circuits show rugged patterns with abrupt "on/off" transitions in expression 313 levels between neighbouring nuclei [23]. In contrast, diffusion-less fully non-autonomous 314 circuits produce smooth spatial expression patterns with a graded increase or decrease 315 in concentration levels across domain boundaries. This is because non-autonomy, with 316 its associated movement of attractors and separatrices over time, provides increased 317 flexibility for fine-tuning expression dynamics over time compared to models with 318 constant phase-space geometry (see below). In biological terms, it suggests that the 319 expression of smooth domain boundaries does not strictly require diffusion. Although 320 diffusion undoubtedly contributes to this process in the embryo, its role may be less 321 prominent than previously thought [23,25]. 322

11/30
Maternal gradient decay affects the level and timing of gap gene 323 expression 324 We used our non-autonomous gap gene circuit to assess the effect of maternal gradient 325 decay on gap gene regulation. One way to isolate this effect is to compare the output of 326 the fully non-autonomous model-with decaying maternal gradients-to simulations 327 using the same model parameters, but keeping maternal gradients fixed to their We do observe, however, that maternal gradient dynamics significantly affect the 334 levels of gap gene expression throughout the trunk region of the embryo (Fig. 4, shaded 335 areas). While early expression dynamics are very similar in both models (time classes 336 C12-T2), they begin to diverge at later stages. The fully non-autonomous model 337 reaches peak expression at T2/T4, but the autonomous model without maternal 338 13/30 gradient decay overshoots observed expression levels in the data between T4 and T8.

339
This indicates that maternal gradient decay leads to decreasing activation rates at the 340 late blastoderm stage, thereby regulating the timing and level of peak gap gene 341 expression. Such a limiting regulatory effect of maternal gradients has been proposed 342 before [25,42], but has never been tested explicitly. without maternal gradient decay (Fig. 4, asterisk). In the posterior, we observe 348 increased levels of Kni and Gt across large parts of their respective expression domains 349 (Fig. 4, arrows). These effects are asymmetric: both posterior Kni and Gt domains  We asked whether the differing effects of maternal gradient decay in the anterior and the 361 posterior of the embryo depend on the presence of different regulatory mechanisms in 362 these regions [23]. To validate this hypothesis, we need to understand and characterise 363 the dynamic mechanisms underlying gene regulation in our non-autonomous model. We 364 achieve this through analysis of the time-variable phase spaces of nuclei across the trunk 365 region of the embryo using the methodological framework presented in the Introduction 366 (S1 Figure B; see [43] for details). To briefly reiterate, this analysis is based on the 367 characterization of the changing phase space geometry that shapes the trajectories of 368 the system. The shape of a trajectory indicates typical dynamical behaviors, that can 369 be classified into four distinct categories-transitions, pursuits, as well as geometrical 370 and topological captures-each showing particular dynamic characteristics. These 371 categories provide mechanistic explanations for the dynamic behavior of the system. For 372 every nucleus, we then compare these non-autonomous mechanisms to the autonomous 373 mechanisms of pattern formation found in the static-Bcd model [23]. This direct 374 comparison allows us to identify the causes underlying the observed effects of maternal 375 gradient decay on the temporal dynamics of gap gene expression.

376
In agreement with Manu et al. [23], we find different patterning modes anterior and 377 posterior to 52% A-P position. Just like in static-Bcd models, anterior expression 378 dynamics are governed by convergence of the system towards attractors in a multi-stable 379 regime. In contrast, our model differs from that of Manu et al. [23] concerning posterior 380 gap gene regulation. We find that a monostable spiral sink drives gap domain shifts in 381 the posterior of the embryo; this differs markedly from the unstable manifold observed 382 in static-Bcd gap gene circuits [23]. An in-depth analysis and biological discussion of 383 spatial pattern formation driven by this mechanism goes beyond the scope of this study. 384 It is provided elsewhere [44]. Here, we focus on temporal aspects of gene regulation and 385 pattern formation, namely the regulation of the velocity of gap domain shifts by  The posterior border of the anterior Gt domain forms between 35 and 40% A-P 394 position (Fig. 5A). Of all the trunk gap genes, nuclei in this region of the embryo only 395 express hb and gt. Gap gene expression dynamics are governed by the same attractor 396 across different nuclei (Fig. 5B). Each trajectory starts at non-zero (maternal) Hb 397 concentration and initially converges towards the attractor located at high Hb and Gt 398 concentrations (Fig. 5B) shaping these trajectories (Fig. 5B, grey trajectories). At some point, the attractor  (Fig.6A). In the non-autonomous model, this boundary interface is set up by two 429 different regulatory mechanisms (Fig. 6B,C). Phase portraits of nuclei between 41 and 430 45% A-P position (Fig. 6A, yellow vertical bar) show the following dynamics: for most 431 of the time, system trajectories converge towards an attractor located at high Hb and 432 high Kr concentration (Fig. 6B, grey trajectory). However, these trajectories are state. At T7, two simultaneous saddle-node bifurcations give rise to two new attractors, 435 one at high Hb and the other at high Kr concentration (Fig. 6B). Two new saddles are 436 also created. System trajectories are caught in the basin of the attractor with high Hb 437 levels. This only has a noticeable effect in more posterior nuclei (e. g. at 43% A-P 438 position in Fig. 6B), where there is a drastic (but late) change in the direction of the 439 trajectory. At T8, two additional saddle-node bifurcations occur, which annihilate the 440 high Hb/high Kr attractor, as well as the newly created attractor at high Hb. This 441 leaves only the attractor at high Kr concentrations (Fig. 6B). Trajectories of the system 442 are once again caught in a different basin of attraction. However, this second round of 443 bifurcations occurs too late to still have a substantial effect on expression dynamics. A 444 16/30 non-autonomous trajectory being caught in a new basin of attraction due to a preceding 445 bifurcation event is called a topological capture (S1 Figure B) [43]. 446 In the region between 46 and 52% A-P position (Fig. 6A, green vertical bar), we 447 observe a different kind of dynamical behaviour. Similar to more anterior nuclei, these 448 instantaneous phase portraits have an attractor at high Hb and high Kr levels and 449 trajectories converge to this steady state at early stages (Fig. 6C). In contrast to more 450 anterior nuclei, however, there is a saddle located on the Hb-Kr plane. Between time 451 class T2 and T6, the position of this saddle moves towards higher Hb levels. When a 452 saddle moves on a phase portrait, it drags the associated separatrix with it (S1 453 Figure A). These concerted movements change the location of the boundaries between 454 existing basins of attraction. When a separatrix "overtakes" a trajectory in phase space, 455 a geometrical capture occurs (S1 Figure B) [43]. This can be observed in the nucleus at 456 47% A-P position (Fig. 6C, grey trajectory). Here, the trajectory gets captured by the 457 moving separatrix between T2 and T5, and later starts to converge towards the  (Fig. 7A) [25,42]. Surprisingly, we find that these shifting posterior gap domains 488 are governed by quite different phase space geometries in our model compared to those 489 previously reported. Manu et al. [23] found that posterior gap gene expression dynamics 490 are controlled by an unstable manifold embedded in a multi-stable phase space 491 geometry in their static-Bcd model. In contrast, our fully non-autonomous gap gene 492 circuit features no such manifold: the phase portraits of posterior nuclei lack saddle 493 points since they are monostable throughout the blastoderm stage and only contain a 494 single attractor (Fig. 7B). This attractor is not a regular point attractor. Its complex 495 eigenvalues reveal that it is a spiral sink (also known as a focus; see S1 Figure A) ( [36]). 496 Like regular point attractors, spiral sinks are stable, in that they draw trajectories 497 asymptotically towards them. Unlike regular point attractors, these trajectories do not 498 approach the steady state in a straight line, but rather spiral inward towards the sink. 499 Since sinks are a consistent feature of the phase spaces of nuclei in the posterior of the 500 embryo, it is likely that they are important for the spiral-shaped geometry of the 501 trajectories observed in this region (Fig. 7B). The spiral geometry in turn is responsible 502 for the ordered succession of transient gap gene expression governing dynamic domain 503 shifts. A full characterization of this patterning mechanism is presented elsewhere [44]. 504 For the purpose of our present analysis, we conclude that the non-autonomous 505 mechanism patterning the posterior region corresponds to a pursuit, where the system 506 follows but never reaches a moving attractor (S1 Figure B).

19/30
The correct geometry of transient trajectories in the posterior of the embryo depends 508 crucially on maternal gradient decay. As we can see in Fig. 7B (black trajectories), 509 simulations without dynamic gradient concentrations show much less tightly wound 510 spirals. This means that the transition between the expression of successive gap genes 511 in this region is delayed. For example, the nucleus at 59% A-P position shows a delayed 512 down-regulation of Kr, while kni keeps on accumulating. This provides a 513 straightforward explanation of the "overshoot" of Kni and Gt domain shifts observed in 514 the simulation without maternal gradient decay (Fig. 7A). In biological terms, it boundaries. As a consequence, gradient decay stabilises spatial gap gene patterns before 533 the onset of gastrulation. An effect of maternal gradient decay on gap gene expression 534 rates has been suggested before-based on the analysis of quantitative expression 535 data [25,42]. However, only mechanistic dynamical models-such as the 536 non-autonomous gap gene circuits presented here-can provide specific mechanisms and 537 quantitative causal evidence for this aspect of gap gene regulation.

538
Our analysis suggests that maternal gradient decay-specifically, the disappearance 539 of Cad from the abdominal region of the embryo-has an important role in regulating 540 the timing of gap gene expression as well as limiting the rate and extent of gap domain 541 shifts in the posterior of the embryo. This result is consistent with experimental data 542 indicating that Cad affects gap domain shifts. Mutants lacking maternal cad, which 543 show a reduced level of Cad protein throughout the blastoderm stage [28], show a delay 544 in the shift of the posterior domains of kni and gt [32,44]. However, Cad does not seem 545 to act exclusively. An indirect role of Bcd in regulating gap domain shifts through 546 altering gap-gap interactions was suggested by a modelling study [30]. It remains 547 unclear whether Cad is also involved in mediating this effect. Finally, a recent study of 548 Bcd-dependent regulation of hb postulated an additional mechanism for gap gene 549 down-regulation that acts before maternal gradient decay occurs [2]. This could have an 550 indirect effect on the timing of late (Bcd-independent) hb regulation, which may mediate 551 the direct effect of Bcd decay on late hb expression we are observing in our models.

552
To better understand the mechanistic basis for the observed differences in patterning 553 between the anterior and the posterior, we analysed the time-variable phase portraits in 554 our non-autonomous model [43]. In agreement with a previous study based on 555 autonomous phase space analysis of static-Bcd gap gene circuits [23], we find that two 556 distinct dynamical regimes govern gap gene expression anterior and posterior to 52% 557 A-P position (Fig. 8). Stationary domain boundaries in the anterior are governed by 558 regulatory mechanisms that are equivalent in static-Bcd and fully non-autonomous 559 models (our work and [23]): they take place in a multi-stable dynamical regime where 560 the posterior boundary of the anterior Gt domain is set by the movement of an attractor 561 in phase space, and the posterior boundary of the anterior Hb domain is set by attractor 562 selection (i. e. the capture of transient trajectories in the non-autonomous case) (Fig. 8, 563 left). Attractor movement in fully non-autonomous models leads to smooth expression 564 boundaries, which are absent in the static-Bcd case. In contrast, static-Bcd and 565 non-autonomous models suggest different mechanisms for gap domain shifts in the 566 posterior of the embryo. While these shifts are controlled by an unstable manifold in 567 the static-Bcd gene circuit model [23], we find a pursuit mechanism featuring a 568 monostable spiral sink to govern their behaviour in our fully non-autonomous analysis 569 (Fig. 8). The spiralling geometry of transient trajectories imposes temporal order on the 570 21/30 progression of gap genes being expressed. If arranged appropriately across nuclei in the 571 posterior of the embryo, this temporal progression from Kr to kni to gt to hb leads to 572 the emergence of the observed kinematic domain shifts [44]. 573 It is important to note that similar regulatory principles can be found in all three 574 solutions of our fully non-autonomous model that reproduce gap-gene patterning 575 correctly both in the presence and absence of diffusion. We have chosen the most 576 structurally stable solution for detailed analysis. The other two circuits show more 577 variability of regulatory features both across space and time. Still, both of these models 578 consistently exhibit multi-stability in the anterior, and spiral sinks as well as transiently 579 appearing and disappearing limit cycles in the region posterior to 52% A-P position.  It is important to note that non-autonomy of the model is not strictly required for 585 the spiral sink mechanism to pattern the posterior of the embryo. Simulations with fixed 586 maternal gradients demonstrate that domain shifts can occur in an autonomous version 587 of our gap gene circuit (see Figs. 4 and 7). The reason why earlier models [22,23] do not 588 feature spiral sinks remains unknown although one possibility is that fitting in the 589 absence of diffusion somehow benefits characterisations of posterior pattern formation in 590 terms of oscillatory behaviours. In spite of this, there are two reasons to consider the 591 mechanism proposed here an important advance over the unstable manifold proposed by 592 Manu et al. [23]. The first reason is technical: non-autonomous gap gene 593 circuits-implementing correct maternal gradient dynamics-are more accurate and 594 stay closer to the data than the previous static-Bcd model. The fact that the quality of 595 a reverse-engineered model usually depends on the quality of its fit to data implies that 596 our model provides more accurate and rigorous predictions than previous efforts. The 597 second reason is conceptual: although it is difficult to interpret an unstable manifold in 598 an intuitive way, it is straightforward to understand the spiral sink as a damped 599 oscillator patterning the posterior of the embryo. The presence of an oscillatory 600 mechanism in a long-germband insect such as D. melanogaster has important functional 601 and evolutionary implications, which are discussed elsewhere [44].

602
Analysis of an accurate, non-autonomous model is required to isolate and study the 603 explicitly time-dependent aspects of morphogen interpretation by the gap gene system. 604 Here, we have shown that such an analysis is feasible and leads to relevant and specific 605 new insights into gene regulation. Other modelling-based studies have used 606 non-autonomous models before (see, for example, [16,26,34,[50][51][52][53]). However, none of 607 them have directly addressed the proposed role of non-autonomy in pattern 608 formation [17]. Our analysis provides a first step towards a more general effort to 609 transcend this limitation in our current understanding of the dynamic regulatory 610 mechanisms underlying pattern formation during animal development.     Figure S2. The three most commonly observed patterning defects in fully non-autonomous diffusion-less gap gene circuits. Commonly observed defects in fully autonomous D. melanogaster gap gene circuits fitted to data without diffusion. Circuits showing any of these gross patterning defects were excluded from further analysis, even if their RMS score was low. Arrows indicate patterning defects as named in the panel headings (A-C). Horizontal axes represent %A-P position (where 0% is the anterior pole). Vertical axes show relative protein expression levels (Rel. Prot. Expr.) in arbitrary units (au). T4/6 indicate time classes C14-T4 and T6, respectively.