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

The variation in the expression patterns of the gap genes in the blastoderm of the fruit fly Drosophila melanogaster reduces over time as a result of cross regulation between these genes, a fact that we have demonstrated in an accompanying article in PLoS Biology (see Manu et al., doi:10.1371/journal.pbio.1000049). This biologically essential process is an example of the phenomenon known as canalization. It has been suggested that the developmental trajectory of a wild-type organism is inherently stable, and that canalization is a manifestation of this property. Although the role of gap genes in the canalization process was established by correctly predicting the response of the system to particular perturbations, the stability of the developmental trajectory remains to be investigated. For many years, it has been speculated that stability against perturbations during development can be described by dynamical systems having attracting sets that drive reductions of volume in phase space. In this paper, we show that both the reduction in variability of gap gene expression as well as shifts in the position of posterior gap gene domains are the result of the actions of attractors in the gap gene dynamical system. Two biologically distinct dynamical regions exist in the early embryo, separated by a bifurcation at 53% egg length. In the anterior region, reduction in variation occurs because of stability induced by point attractors, while in the posterior, the stability of the developmental trajectory arises from a one-dimensional attracting manifold. This manifold also controls a previously characterized anterior shift of posterior region gap domains. Our analysis shows that the complex phenomena of canalization and pattern formation in the Drosophila blastoderm can be understood in terms of the qualitative features of the dynamical system.

The result confirms the idea that attractors are important for developmental stability and shows a richer variety of dynamical attractors in developmental systems than has been previously recognized.

Canalization refers to the constancy of the wild type phenotype under varying developmental conditions [1]–[4]. In order to explain canalization, C. H. Waddington hypothesized that there must only be a finite number of distinct developmental trajectories possible, since cells make discrete fate decisions, and that each such trajectory, called a chreod, must be stable against small perturbations [5]. One aspect of canalization, the buffering of phenotypic variability against genotypic variability in wild type, has received considerable experimental [2], [6]–[10] and theoretical [11]–[13] attention. The phenomenon of canalization of genotypic and environmental variation was seen by Waddington as a consequence of the underlying stability of developmental trajectories, an idea supported by theoretical analysis [13]. But this central idea of Waddington's has heretofore received little attention in real developmental systems because of a lack of relevant quantitative molecular data. The further investigation of Waddington's hypothesis is of great importance because it provides a scientific connection between the reliability and invariance of the formation of cell types and tissues in the face of underlying molecular variability, as we now explain.

Quantitative molecular data permitting the study of developmental canalization are now available for the segment determination process in Drosophila [14]. The segmented body plan of the fruit fly Drosophila melanogaster is determined when the embryo is a blastoderm [15] by the segmentation genes [16]. Quantitative spatiotemporal gene expression data show that the maternal protein gradients and the early expression patterns of the zygotic gap and pair-rule genes vary a great deal from embryo to embryo [14],[17]. The variation of the expression patterns of the gap and pair-rule genes decreases over time so that it is significantly lowered by the onset of gastrulation at the end of cellularization ([14], Fig. 1). The observed reduction of variability over time in the segmentation gene system suggests that the developmental trajectory of the Drosophila embryo is stable against perturbation. The characterization of the stability properties of the developmental trajectory is central to our understanding of the mechanisms that underlie canalization [3].

In the case of the gap genes, we have shown elsewhere [18] that variation reduction relative to the maternal gradient Bicoid (Bcd) occurs because of gap gene cross regulation. Using a gene circuit model of the gap gene network [18]–[22] we identified specific regulatory interactions responsible for variation reduction in silico and verified their role in canalization experimentally. Importantly, the model reproduces the observed low variation of gap gene expression patterns [18], which provides an opportunity to analyze the properties of the system that give rise to developmental stability.

These results raise two generic problems that occur in the analysis of complex numerical models. First, even if the model describes a natural phenomenon faithfully, understanding the natural phenomenon is only achieved when the model's behavior can be understood as well. The complexity of the model, unsurprising in terms of the underlying complexity of the biological system itself, poses a significant challenge to understanding model function. Second, any model is an approximation to the actual mechanisms operating in an organism. The model's behavior must be robust to perturbation, since organisms develop and function reliably even though the underlying mechanisms are subject to a wide variety of perturbations and stresses. There is extensive molecular variability among cells and embryos ([14], [17], [23]–[26]; reviewed in [27]) and yet there is functional identity between equivalent cell types or conspecific individuals.

René Thom tried to resolve this apparent contradiction between the constancy of biological function and the variability in biological substructure by proposing a qualitative topological view of the trajectories of dynamical models [28]. The term “topology” is used here to refer to properties of developmental trajectories that are invariant under continuous deformation. The preservation of these properties ensures the robustness of model behavior, while a qualitative view often leads to an intuitive understanding of complex mechanisms. One such robust property is an attractor state, or a stable steady state of a dynamical system, that attracts all trajectories in some neighborhood of itself. Attractor states are locally stable under small perturbations of the dynamical model [29], and for this reason it has been proposed that cell fates are attractors [11], [13], [30]–[33].

The presence of an attractor state in the phase space of a system implies that there exists a region of phase space, called the basin of attraction, in which all trajectories approach the attractor asymptotically [34],[35]. This suggests that an attractor is the kind of qualitative robust property that could explain the stability of trajectories, and hence canalization. There are, however, three important considerations to keep in mind when using attractors to describe the Drosophila blastoderm. First, the reduction of variation due to attractors is only guaranteed at late times, but the reduction in the variation of the gap gene expression patterns takes place over about 100 minutes prior to gastrulation. The reduction of variation before gastrulation is biologically essential as the expression patterns of engrailed and wingless, which form the segmentation prepattern, have a resolution of one nucleus and are created by the precise overlap of pair-rule and gap domains [14],[36]. Furthermore, at about the time of gastrulation the embryo undergoes the midblastula transition [37],[38] at which time a qualitative change occurs in the genetic control of the embryo. Second, in general there can be more than one attractor in the phase space [39]–[43]. Thus, the basins of attraction need to correspond to biological initial conditions and be large enough to ensure robustness. Finally, the set of attractors found must succeed not only in explaining canalization but also the morphogenetic properties of the system. One such property is the anterior shift of gap gene domains located in the posterior region [14],[21],[44]. These shifts are biologically significant and are difficult to reconcile with stable point attractors.

In this paper we show that the variation reduction of gap gene expression patterns is a consequence of the action of robust attracting states. We further show that the complex patterning of the gap gene system reduces to the three qualitative dynamical mechanisms of (1) movement of attractors, (2) selection of attractors, and (3) selection of states on a one dimensional manifold. The last of the three mechanisms also causes the domain shifts of the gap genes, providing a simple geometric explanation of a transient phenomenon.

In the Gap Gene Circuits section we briefly describe the gene circuit model; see [18] for a full description. For each nucleus in the modeled anteroposterior (A–P) region, we identified the attractors in the gap gene phase space, calculated the trajectories, the basins of attraction and other invariant sets such as one dimensional attracting manifolds (Stability Analysis of the Trajectories of the Gap Gene System section). The stability of the trajectories was tested by varying the initial conditions within a biological range, based on gene expression data, that represents the variability of early gap gene expression. We plotted the attractors and several trajectories corresponding to different initial conditions to make phase portraits that show the global qualitative behavior of the system. Finally, we studied how the phase portraits changed as A–P position was varied to infer qualitative pattern formation mechanisms. The biological conclusions about canalization and pattern formation arising from the dynamical characterization are presented in the Mechanisms of Canalization and Pattern Formation section.

The gene circuit used in this study models the spatiotemporal dynamics of the protein expression of the gap genes hunchback (hb), Krüppel (Kr), giant (gt), and knirps (kni) during the last two cleavage cycles (13 and 14A) before gastrulation [37] in the Drosophila blastoderm. The protein products of these genes localize to nuclei [45]–[48] so that the state variables are the concentrations of the proteins in a one dimensional row of nuclei along the A–P axis of the blastoderm. The concentration of the protein in the nucleus at time is denoted by . In the model we considered a region, from 35% to 92% egg lenth (EL) along the A–P axis, which corresponds approximately to the region of the blastoderm fated to form the segmented part of the adult body [49],[50].

The gap genes are expressed in broad domains (Fig. 2A,B; [14]) under the control of maternal cues. The anterior maternal system acts primarily through the protein gradient Bcd [51]–[53] which is essentially stationary and has an exponential profile (Fig. 2C; [14],[51],[54]) during the modeled time period. The posterior maternal system is represented by the maternal Hb gradient (Fig. 2C; [55]–[57]). The terminal system regulates gap gene expression by activating tailless (tll) and huckebein (hkb) [58]–[61]. The terminal system is represented in the model by the Tll gradient, which is expressed posterior to 80% EL in the modeled region during cycles 13 and 14 ([14] and Fig. S1B). tll is considered upstream of the gap genes since its expression pattern is unchanged in gap gene mutants [62]. The concentration of Bcd in nucleus is denoted by and was determined using Bcd data from a representative cycle 13 embryo by an exponential fit, so that (see [18] for details). The concentrations of Tll and another upstream regulator, Caudal (Cad) [63],[64], were determined by interpolating average data in time [18]. The concentrations of Tll and Cad are denoted by respectively, with an explicit dependence on time, since these gradients are not stationary (Fig. S1).

where in a gene circuit with genes and nuclei. The first term on the right hand side of Eq. (1) represents protein synthesis, the second one represents protein transport through Fickian diffusion and the last term represents first-order protein degradation. The diffusion constant, varies inversely with squared internuclear distance, and is the degradation rate. The synthesis term is set to zero during the mitosis preceding the thirteenth nuclear division as synthesis shuts down [65]. Following this mitosis, the nuclei are divided and daughter nuclei are given the same state as the mother nucleus.

is the maximum synthesis rate, and is a sigmoidal regulation-expression function. The first term in the argument of represents the transcriptional cross regulation between the gap genes and the genetic interconnectivity is specified by the matrix . Positive elements of imply activation while negative ones imply repression. The regulation of the gap genes by Bcd is represented in the second term and is the regulatory strength. The regulation of the gap genes by upstream time-varying inputs is represented in the third term and is the number of such inputs. There are two such inputs in this model, Cad and Tll, and the elements of the matrix have the same meaning as those of . The last term, , represents the effect of ubiquitous transcription factors and sets the threshold of activation.

The initial conditions for Hb are specified using cleavage cycle 12 data. Cycle 12 data are a good approximation to the maternal Hb gradient since the zygotic expression of hb appears to begin in cleavage cycle 13 [17]. The initial conditions for Kr, Gt, and Kni are taken to be zero, since their protein expression is first detected in cycle 13 [14], [61], [66]–[68].

The gene circuit's parameters were determined by performing a least-squares fit to a time series of averaged gap gene data [14] using the Parallel Lam Simulated Annealing algorithm (see Methods). This time series has nine points (time classes; see Table S1), one in cycle 13 and the rest in cycle 14A. The output of the gene circuit (Fig. 2E,F) fits the data (Fig. 2D) well and its network topology (Fig. 2K) is consistent with previous results (see [18] for discussion and parameters).

In order to characterize the stability of the trajectories of the gap gene system in terms of qualitatively robust features like attractors, we apply the tools of dynamical systems theory [34],[69]. Since the gene circuit has variables (Gap Gene Circuits section) its state is represented as a point in an -dimensional concentration space, or phase space. In general the concentrations of gap proteins change with time, and hence, a solution of the gene circuit is a curve in this phase space. The gene circuit can also have solutions which do not change with time. Such a solution, called an equilibrium or steady state solution, is represented as a single point in phase space. The positions of the equilibrium solutions in phase space and their stability properties determine the stability of a general time varying solution of the gene circuit. The reader not familiar with linear analysis near an equilibrium point should see Protocol S2 for a pedagogical description of equilibria and their stability in two dimensions.

In the gap gene circuit used in this study (Gap Gene Circuits section), there are 58 nuclei and 4 gap genes, giving rise to a dynamical system having 232 dimensions. Such a large number of dimensions pose a significant challenge to visualizing the results of the phase space analysis. In order to make the analysis tractable, we made the approximation that there is no diffusion, in Eq. (1), while keeping all other parameters in Eq. (1) at their original values [18]. This uncouples the nuclei and the system of coupled ODEs reduces to a set of independent systems of ODEs. Eq. (1) thus becomes

where the dependence of the concentrations of gap gene proteins, Bcd and time varying external inputs on the location of nucleus along the A–P axis is denoted by .

In the absence of diffusion, the model still gives the correct sequence of gap gene domains (Fig. 2G,H); and the border positions are close to their values in the presence of diffusion (Table S2). In contrast to the circuit with diffusion however, the domains have sharp boundaries and have little or no overlap (Fig. 2H). The general agreement between the circuits with and without diffusion allows the analysis of the reduced system and supports the result from earlier work [21] that diffusion is not required for making gap gene expression patterns.

In Eq. (2), the terms represent the anterior and terminal maternal systems. They parametrize the set of solutions possible in each nucleus as a function of A–P position. From the set of solutions specified by the concentrations of Bcd, Cad, and Tll, a particular solution is chosen by the initial conditions specified by the concentration of maternal Hb (Gap Gene Circuits section; Fig. 2C). Therefore the posterior maternal system specifies position in a manner distinct from the anterior and terminal systems. We further simplified the analysis by neglecting the effects of Tll in Eq. (2), giving

The gap gene expression patterns (Fig. 2I,J) produced by Eq. (3) agree with a model that includes tll (Eq. 2) until 71% EL, therefore the simplified analysis is valid in the region from 35% EL to 71% EL.

We considered each nucleus in cleavage cycle 13 and its posterior daughter in cycle 14A since the anterior daughter differs only slightly from the posterior one in the concentrations of Bcd and Cad. For simplicity we refer to mother-daughter pairs by the daughter nuclei, which are located at odd numbered positions from 35% to 71% EL. For each nucleus we determined the equilibria and, based on maternal Hb data (Fig. 2C), the family of solutions corresponding to the observed range of initial conditions. The concentration of Cad () varies with time (Fig. S1A), implying that the equilibria in the phase space of a nucleus also change with time. We are interested in studying the properties of the solutions during late cycle 14A, since the variation in gap gene expression levels is least in that time period (Fig. 1). Also, the stationary approximation for Bcd is valid until time class T6 (∼16 mins before gastrulation; Table S1) as Bcd concentration decreases two-fold afterward [14]. For these reasons, we calculated equilibria at time class T6 (Protocol S1: Eq. S3). The trajectories were calculated using time varying Cad data until time class T6 and with T6 Cad data thereafter (Protocol S1: Eq. S1) so that their late time behavior corresponded to the calculated equilibria.

The equilibria were calculated using the Newton-Raphson method [70],[71] and classified according to their stability (Protocol S3). We calculated several trajectories with different initial conditions to test variation reduction and stability, and continued the integration of Eq. (S1) (Protocol S1) to very late times in order to visualize asymptotic behavior. Since the segments are determined by the time of gastrulation, significant variation reduction must occur beforehand for the precise specification of position. Hence, we distinguish the biological behavior from asymptotic behavior graphically. We found that each nucleus has multiple attractors and that its trajectory can potentially approach any one depending on the initial condition. It was necessary, therefore, to characterize the basins of attraction [34],[69] of the attractors. We exploited the fact that only Hb has non-zero initial conditions (Gap Gene Circuits section) to characterize the basins as intervals on the Hb axis (see Protocol S3 for details of calculation). We also calculated one dimensional unstable manifolds [72] of saddle points to better understand the transient behavior of solutions (Protocol S3).

The analysis of the phase space for different nuclei revealed several A–P positions at which the set of equilibria changes. This suggests that the dynamical system is undergoing bifurcations resulting from the changing values of its parameters, the concentrations of Bcd and Cad, as A–P position is varied. Since we used Cad data in Eq. (3), it is also possible that these bifurcations are artifacts caused by fluctuations in the data. We verified that the changes observed in the analysis at discrete odd numbered positions were actual bifurcations by varying the concentrations of Bcd and Cad, parametrized by the A–P position, continuously (Protocol S4). This was done by interpolating the Cad data with a cubic polynomial (Fig. S3). Note that the Bcd concentration profile has already been parametrized as a function of position (; Gap Gene Circuits section), and can be varied continuously. The continuous analysis validated all the bifurcations at discrete locations (Fig. 3A and Fig. S4; circles) save for two. These two bifurcations (Fig. 3A; boxes) could not be observed in the analysis at discrete locations as they occur consecutively between two nuclei and leave the set of equilibria unchanged. Although the bifurcations occur at particular values of the concentrations of Bcd and Cad, subsequently we will refer to them by A–P position for simplicity.

We further characterized the type of bifurcations (Protocol S4). Though many types of bifurcations are possible only one type, the saddle-node (see Fig. S2 for illustration in two dimensions), occurs in the gap gene system. In the four-dimensional gap gene system any creation or annihilation of a pair of equilibria that differ in the sign of one eigenvalue is a saddle-node bifurcation (Fig. 3A and Fig. S4).

Based on the analysis in the previous section the region of interest, from 35% EL to 71% EL, can be divided into an anterior and a posterior region (Fig. 2J) having distinct modes of canalization and pattern formation. The two regions are separated by a saddle-node bifurcation that occurs at 53% EL (Fig. 3A), that is, at the peak of the central Kr domain.

We next demonstrate that in the anterior region (Fig. 2J), which extends from the peak of the third anterior gt domain to the peak of the central Kr domain, the state of a nucleus at gastrulation is close to a point attractor. The trajectories are stable by virtue of being in the basin of attraction of the nucleus's attractor state and hence canalize. Pattern formation occurs by the selection of one state from many in a multistable phase space. The concentrations of the Bcd and Cad gradients control pattern formation in the anterior by determining the sizes of the basins and the positions of the attractors, while maternal Hb concentration selects a particular attractor by setting the starting point in its basin. Previous experimental [73] and theoretical [19] work suggested that Bcd and maternal Hb patterned the anterior of the embryo synergistically; our results identify specific roles for Bcd and Hb in anterior patterning.

The posterior region extends from the peak of the central Kr domain to the peak of the posterior gt domain (Fig. 2J) and its nuclei have phase spaces with very different properties. In this region, the state of the nucleus is far from any attractor state at gastrulation. Instead the state of a nucleus is close to a one-dimensional manifold and canalization is achieved due to attraction by this manifold. Even though the phase space is multistable, the biological range of maternal Hb concentrations in the posterior region place all nuclear trajectories in one basin of attraction. As a consequence, the modes of pattern formation operative in the anterior cannot function in the posterior. Maternal Hb patterns the posterior by determining the position on the attracting manifold which a particular trajectory reaches by the time of gastrulation. These results reveal the mechanism by which maternal Hb acts as a morphogen in the posterior [73]–[75] and also explain the dynamical shifts of gap gene domains [14],[21], a significant biological property of the posterior region.

We begin the presentation of detailed results by describing the phase spaces of typical nuclei in the two regions, highlighting mechanisms for canalization and pattern formation. An equilibrium is labeled by either (point attractor) or (saddle equilibrium), denoted by a superscript, with subscripts denoting the number of eigenvalues having positive or negative real parts. For example denotes the second saddle equilibrium in the modeled region which has one eigenvalue with positive real part and three with negative real parts. Equilibria are also given descriptive names based on which proteins are at high levels (on) ignoring the proteins that are at low levels. For example, if a point attractor is at hb-on, Kr-off, gt-on, and kni-off, it is referred to as the “hb,gt-on” attractor.