Canalization of gene expression and domain shifts in the Drosophila blastoderm by dynamical attractors

Fig. 4A shows the phase portrait of the nucleus at 37% EL from the anterior region. Since the phase space is four-dimensional, a three-dimensional projection onto the Hb-Kr-Gt axes is shown (see Fig. S6 for other projections). There are three attractors, (hb,gt-on), (hb,Kr-on), and (Kr-on) shown in blue. Two saddle equilibria, , are shown in red. Saddle equilibria that are unimportant for the dynamics are not shown (see Protocol S5 and Table S4). In order to show the family of trajectories allowed by the values of Bcd and Cad concentrations in this nucleus, trajectories are plotted with starting points distributed uniformly on the Hb axis between 0–100, corresponding to the observed range of maternal Hb in the modeled region [14]. Time is represented as a color gradient along the trajectories, with the start of cycle 13 green, and gastrulation red. Trajectories are blue after gastrulation to indicate their asymptotic behavior.

The red-to-blue transition on the trajectories occurs very close to the attractors (Fig. 4A) implying that the nucleus is close to equilibrium at gastrulation. It is apparent by visual inspection that the initial variation of the trajectories is reduced dramatically by the onset of gastrulation, that is, trajectories are stable due to the point attractor. To verify this property numerically, we constructed a four-dimensional box such that its dimensions corresponded to the range of gap protein concentrations observed in cycle 13 and its volume represented initial variation. We then calculated the time evolution of the volume of this box using the model (Fig. 4C; see Protocol S8 for details of calculation). Variation in the concentration of each protein shrinks by a factor of ∼20 on average by gastrulation, confirming the stability of the developmental trajectory.

The concentration of maternal Hb, depending on which basin of attraction it lies in, selects a particular attractor state. The boundaries between the basins of attraction of the three attractors are indicated by arrows in Fig. 4A (see Protocol S5 for details). The concentration of maternal Hb in the nucleus at 37% EL puts it in the basin of the hb,gt-on () attractor (Table 1). This analysis correctly accounts for the observed gap gene expression in the nucleus, which lies at the peak of the third anterior gt domain and in the anterior hb domain (Fig. 2F,H). Posterior to this nucleus and up to 45% EL, all nuclei remain in the basin of the attractor (Table 1). changes its position in the phase space from hb,gt-on at 37% EL to hb-on at 43% EL (Table 1 and Fig. 3B). The posterior border of the third anterior gt domain therefore forms by the movement of in the phase space. Since the Bcd and Cad concentrations are the only parameters being varied, the movement of an attractor is one mechanism by which Bcd patterns the anterior.

The nucleus at 47% EL is in the basin of the Kr-on () attractor, unlike the nuclei to its anterior, which are in the basin of hb-on (Table 1). The switch from hb-on to Kr-on causes the formation of the posterior border of the anterior hb domain and the anterior border of the central Kr domain between 45% EL and 47% EL (Fig. 2F,H). This happens due to 1) an increase in the size of the basin of the Kr-on attractor and 2) a decrease in the concentration of maternal Hb (second and fifth columns of Table 1). As is the case with the movement of the attractor, the first effect is due to the changing value of Bcd concentration. Hence, Bcd and maternal Hb pattern the anterior region in tandem [19],[73]; Bcd has the role of setting attractor positions and the extent of basins and maternal Hb selects a particular basin.

The nuclei remain in the basin of the Kr-on attractor up to 53% EL, where a saddle-node bifurcation annihilates Kr-on (Fig. 3A). Although there are other bifurcations in the A–P region considered (Table S3), only this bifurcation affects the dynamics of the gap gene system in a significant way. Posterior to this bifurcation the remaining two attractors hb-on () and hb,Kr-on () persist to the end of the region being analysed. No new attractors that might correspond to gap gene expression domains of the posterior appear, suggesting that the mechanisms of canalization and pattern formation are different from the anterior.

Fig. 4B shows the Hb-Kr-Kni projection (see Fig. S8 for others) of the phase portrait of the nucleus at 57% EL in the posterior region. There are two attractors (hb-on) and (hb,Kr-on). Two saddle equilibria, , are shown, see Protocol S6, Fig. S7, and Table S4 for others. The basin of attraction of is very small ((43.99, 44.04)) and is shown with an arrow. The basin of attraction of is divided into two intervals by the basin of . Trajectories right of the arrow (“direct interval”) reach close to by gastrulation. Trajectories from the “indirect interval” to the left of the arrow are not close to the attractor at gastrulation as can be seen from the red-to-blue transitions. The direct trajectories are nonbiological since the concentration of maternal Hb (Fig. 2C) in the posterior region nuclei places all biological trajectories in the indirect interval.

Although the biological indirect trajectories are not in the vicinity of an attractor at gastrulation, they appear to be converging to a single trajectory in the phase space. This trajectory can be visualized by following the blue (post-gastrulation) segments of the trajectories in Fig. 4B. Therefore the trajectories are demonstrating stability, which we verified numerically by calculating the time evolution of a volume representing initial variation (Fig. 4D). The average variation in the concentration of each protein reduces by a factor of ∼10 by gastrulation. In order to understand this puzzling stability further, we calculated the one-dimensional unstable manifolds of the saddle equilibria . The unstable manifold of , is shown as a translucent magenta tube in Fig. 4B. This manifold is precisely the trajectory to which all indirect trajectories are converging. Furthermore, the trajectories are close to this manifold by gastrulation (red-to-blue transitions). Hence this attracting manifold plays the same role in the posterior as attractors do in the anterior, and is responsible for the stability of trajectories and canalization (see also Protocol S7).

The attracting manifold is also important for pattern formation in the posterior region. The positions and stability of equilibria in the phase space are qualitatively the same in all nuclei of the posterior region. This is a reflection of the relative constancy of Bcd, which is shallow due to exponential decay with position, and Cad, which has almost uniform expression in the posterior region (Fig. 2C; [14]). The invariance of the qualitative properties of the phase space implies that Bcd does not provide positional information in the posterior, a result arrived at by different means in other work [18],[21].

In order to understand how pattern formation occurs in the posterior, it is important to note two properties of the phase space (Fig. 4B). First, the attracting manifold traverses all the posterior region states (Kr-on, kni-on, gt-on, and intermediate values; see Fig. 4E), before reaching . Second, the state achieved by an indirect trajectory at gastrulation has a continuous dependence on initial Hb value. High values of maternal Hb (20–40) in a nucleus lead to a Kr-on state, intermediate levels (12–20) lead to a kni-on state and lower values (4–8) lead to a gt-on state. Since maternal Hb decreases monotonically with A–P position (Fig. 2C), these two properties lead to the formation of gap gene expression patterns in the posterior region. In fact, changing only the maternal Hb concentration in a single posterior region nucleus (Bcd concentration is fixed), produces correct gap gene expression patterns in the posterior region (Fig. 5). This explains the experimental result that maternal Hb alone can pattern the embryo, if the anterior and terminal systems are absent, producing the posterior region expression patterns over the length of the whole embryo [75]. There is a further decrease in variability within the subintervals of maternal Hb concentration that lead to a particular gap gene state at gastrulation (see Protocol S7 and Fig. S9 for details). These results establish that maternal Hb is the morphogen in the posterior and patterns the region by selecting states on the attracting manifold .

Finally, the attraction of trajectories by the manifold provides a mechanism for dynamical shifts of gap gene domains [21]. We illustrate this mechanism with the nucleus through which the Kr posterior border and the Kni anterior border pass as they shift anteriorly (Fig. 6B,D). Fig. 6A shows the trajectory of the nucleus at 59% EL. The trajectory starts on the Hb axis, and is attracted to , reaching close to the kni-on state at gastrulation. However in approaching , it goes through intermediate states. First Kr increases, and then Kni increases with a concomitant reduction in Kr (see Fig. 6C). We note that the trajectory is a result of asymmetric repression of anterior domain genes by posterior domain proteins [21]. A nucleus in the Kr-on state goes to the Kni-on state because Kni's repression of Kr is stronger than that of kni by Kr. A similar mechanism applies to nuclei at the Kni posterior and Gt anterior boundaries (data not shown).

A discrete [6],[7] and buffered response to perturbations is the hallmark of a canalized developmental system. Without recourse to molecular data, Waddington sought to explain these two properties of the response by postulating certain favored stable developmental trajectories which he called chreods. Our results (see Fig. 7 for summary) show that dynamical systems with multiple attracting states possess both of these properties. Small perturbations are damped because of phase space volume contraction driven by attractors. A discrete response to larger perturbations is a consequence of the discontinuous boundaries between the basins of attraction of a multistable system or of bifurcations. Using a model based on gene expression data, we can conclude that the trajectory of the gap gene system is a chreod.

The initial high variation of gap gene expression may arise from early events governed by stochastic laws. Previous observations indicate that the first nuclei in which gap gene transcription is activated are selected probabilistically [68],[76]. Moreover, the gradients of Hb and Cad proteins are formed by translational repression from the Nanos and Bcd gradients respectively [57],[77],[78] under conditions of relatively low molecular number [23], which is likely to lead to intrinsic fluctuations [79]. Our results show that a deterministic description of gap gene dynamics is sufficient to account for the reduction of initial variation regardless of its source. It is evident however that there are at least two other types of variation that the system might be subject to. First, a natural population will have genotypic variation which, in the framework of the model, would be reflected in the variation of its parameters. Second, gap gene expression itself is likely to be a stochastic process rather than a deterministic one. Notwithstanding this fact, there is no evidence in Drosophila for the coupling of molecular fluctuations to phenotypic fluctuations as seen in prokaryotes [80], suggesting that molecular fluctuations are buffered in some sense. We emphasize that an attractor is stable against small perturbations of the model itself [29], and hence is a model property that is preserved to an extent if there is genotypic variation in a population or if errors are introduced by stochastic gene expression. However, further study of both of these aspects of canalization is required in order to more fully understand their role.

With regard to pattern formation in the blastoderm, the prevailing theory is that the border positions of downstream genes are determined at fixed values or thresholds of the Bcd gradient [23],[52]. This idea cannot, however, account for either the low variability of downstream gene border positions [14],[17],[18], or the dynamical shifts of domains in the posterior [14],[21]. Fixed threshold specification also cannot explain precise placement of the borders in the posterior since the low molecular number of Bcd in the nuclei implies a high level of molecular noise [23],[81]. In the dynamical picture (Fig. 7), contrary to the threshold view, Bcd ceases to have a role in positional specification posterior to the peak of the Kr domain since, posterior to this position, the geometry of the phase space does not change qualitatively with A–P position. Instead, maternal Hb acts as a morphogen, obviating the problems arising from a low molecular number of Bcd. Maternal Hb has long been recognized as a morphogen [74],[75] for the posterior region but the mechanism with which it specifies the posterior region pattern was not clear. As is the case with Bcd, a threshold-based theory for positional specification by Hb [82] is incomplete and requires the postulation of thresholds that can be modified by their targets. The qualitative dynamics provides a viable mechanism for posterior patterning. The attracting manifold is the geometric manifestation of asymmetric repression between the gap genes in reverse order of gap gene domains, . The initial Hb concentration determines which neighborhood of the manifold the trajectory traverses as it is reaches the manifold: Kr-on, kni-on, or gt-on. In other words, posterior patterning works by triggering particular feedback loops in the gap gene network based on maternal Hb concentration.

This mechanism also accounts for domain shifts, a property particular to the posterior region, since the trajectories mimic the geometry of the manifold as they approach it.

The dynamical analysis of the gap gene system provides a simple and integrative view of pattern formation in the blastoderm (Fig. 7). The existence of distinct anterior and posterior patterning systems was inferred from the effect of maternal mutations on larval cuticle phenotype and was subsequently characterized in terms of the effects of the Bcd [52],[73],[83],[84] and maternal Hb gradients [56],[57],[77]. But where and how is the control of patterning transferred from Bcd to maternal Hb? Our analysis shows that the hand-off occurs at the A–P position where the Kr-on attractor is annihilated through a saddle-node bifurcation, implying a sharp rather than gradual transfer. With knowledge of the two dynamical regimes, the complex spatiotemporal dynamics of the gap gene system can be understood in the simple terms of three mechanisms: movement of attractors through phase space, selection of attractors by initial conditions, and the selection of states on an attracting manifold (Fig. 7).

Finally, we mention the advantage of having the unexpected mechanism of a one dimensional manifold for canalization and patterning. The Bcd concentration is a bifurcation parameter of the dynamical equations. If there were specific attractors corresponding to each gap gene state, with bifurcations creating and annihilating them successively as the Bcd concentration is varied, the molecular noise in Bcd [23] would give rise to “jitter” or rapid switching between attractors. The manifold with its smooth dependence on maternal Hb is qualitatively robust to such fluctuations. In a connectionist model of cognition [85], one dimensional unstable manifolds connecting a sequence of saddle points have been proposed as a means of representing transient brain dynamics. The gap gene phase space is a low dimensional projection of the high dimensional phase space of all the molecular determinants in the blastoderm. It may well be that the attractors found in our analysis are actually saddle points in the high dimensional phase space and are way points, with manifolds connecting them, rather than final end points.

The methods used to obtain and characterize the quantitative data are as described in earlier work [14]. All gene expression levels are on a scale of 0–255 chosen to maximize dynamic range without saturation. The numerical implementation of the gene circuit equations is as described [18],[21]. The gap gene circuit was fit to integrated gap gene data [14] using Parallel Lam Simulated Annealing (PLSA) [86],[87]. PLSA minimizes the root mean squared (RMS) difference between model output and data. For each nucleus, data were available at nine time points (Table S1). Search spaces, penalty function, and other annealing parameters were as described [22],[88]. The circuit analyzed in detail had an RMS score of 10.76, corresponding to a proportional error in expression residuals of about 4–5%.

Equilibria were determined by the Newton-Raphson method as described in Protocol S3. One-dimensional unstable manifolds of hyperbolic equilibria were calculated by solving the ODEs using the Bulirsch-Stoer [71] method with starting points in the unstable eigenspace of the equilibria [72]. The basin boundaries on the Hb axis were calculated by finding starting points for trajectories that reach saddle points with one positive eigenvalue (Protocol S3). The time evolution of volume phase space was calculated as described (Protocol S8). The methods used to calculate the equilibria branches and to determine the type of bifurcations are described in Protocol S4.

We thank R. K. Lewis, J. P. Gergen, J. Jaeger, K. Krishan, C. Martinez, and A. M. Samsonov for valuable discussions and several anonymous reviewers for helpful comments on the manuscript.