Abstract
Background
Extensive variation in early gap gene expression in the Drosophila blastoderm is reduced over time because of gap gene cross regulation. This phenomenon is a manifestation of canalization, the ability of an organism to produce a consistent phenotype despite variations in genotype or environment. The canalization of gap gene expression can be understood as arising from the actions of attractors in the gap gene dynamical system.
Results
In order to better understand the processes of developmental robustness and canalization in the early Drosophila embryo, we investigated the dynamical effects of varying spatial profiles of Bicoid protein concentration on the formation of the expression border of the gap gene hunchback. At several positions on the anteriorposterior axis of the embryo, we analyzed attractors and their basins of attraction in a dynamical model describing expression of four gap genes with the Bicoid concentration profile accounted as a given input in the model equations. This model was tested against a family of Bicoid gradients obtained from individual embryos. These gradients were normalized by two independent methods, which are based on distinct biological hypotheses and provide different magnitudes for Bicoid spatial variability. We showed how the border formation is dictated by the biological initial conditions (the concentration gradient of maternal Hunchback protein) being attracted to specific attracting sets in a local vicinity of the border. Different types of these attracting sets (point attractors or one dimensional attracting manifolds) define several possible mechanisms of border formation. The hunchback border formation is associated with intersection of the spatial gradient of the maternal Hunchback protein and a boundary between the attraction basins of two different point attractors. We demonstrated how the positional variability for hunchback is related to the corresponding variability of the basin boundaries. The observed reduction in variability of the hunchback gene expression can be accounted for by specific geometrical properties of the basin boundaries.
Conclusion
We clarified the mechanisms of gap gene expression canalization in early Drosophila embryos. These mechanisms were specified in the case of hunchback in well defined terms of the dynamical system theory.
Background
Development is surprisingly robust to environmental stress, intrinsic fluctuations, and genetic variability in populations. These facts, together with the observation that cell type is a discrete rather than continuous property, led C. H. Waddington to propose that developmental processes have innate errorcorrection properties, which he called "canalization" [1]. Waddington visualized error correction in terms of an "epigenetic landscape," in which the developmental state of an organism is analogous to a ball rolling down a sloping landscape containing multiple "hills" and "valleys": as development progresses, cells take different paths down this landscape and so adopt different fates. Uncontrolled differentiation does not occur because the hills act as barriers and the state remains near a valley floor. This picture has natural corollaries in terms of genetic variability and evolution. Although the shape of the landscape may alter slightly in the face of genetic variation, the tendency of the system to stay near the valley floor will buffer the phenotypic consequences of these variations. Even under larger changes over evolutionary time, the tendency of the system to stay close to the valley floors will preserve and buffer developmental pathways in the face of evolutionary change.
We recently demonstrated the existence of canalization at the molecular level in the segment determination system of the fruit fly Drosophila melanogaster [2]. We showed that with respect to the gap gene system, canalization was a consequence of gap gene crossregulation [3], and furthermore that it was associated with error correction by dynamical attractors [4], a precise mathematical formulation of Waddington's ideas about hills and valleys [5]. That work remains incomplete for reasons involving both broad biological considerations and specific mathematical points. With respect to large scale biological issues, the fact that canalization is under the control of natural selection means that all regulatory systems in an organism are in some sense selected for canalizing properties. In a general evolutionary context, canalization tends to appear as a multigenic trait involving the buffering of underlying genetic variation [68]. In a specific example concerning the segmentation system, it was shown that the differing proportional placement of evenskipped stripes in three lines and two species of Drosophila [9] depended on differences in the maternally expressed genes of flies in these lines, rather than the zygotic gap genes, the crossregulation of which we have shown to ensure proportional spacing [3] because of the alternating arrangement of strongly mutually repressing gap gene expression domains [10].
In comparing our results to those of Lott and coworkers, we note that we considered the canalizing behavior of a specific zygotic component of the segmentation system, while the study of Lott et al. looked at the canalizing behavior of the full set of organismal genes, many subsets of which are presumably engaging in their own form of canalizing behavior. Closing the gap between the genetic control of canalization at the level of populations versus the well defined small networks of genes considered in developmental genetics will require finding the specific genes responsible for population effects as well as taking the complementary step of incorporating additional mechanisms and genes into the well characterized systems arising in a developmental genetics context.
In [3,4], we used the "gene circuit" approach to demonstrate that the reduction in variance of gap gene expression was a consequence of gap gene cross regulation. Gene circuits [1113] are dynamical models that can reproduce observed gene expression patterns by reconstituting the required set of genetic interactions in silico. Our model accounts for the expression of four mutually interacting gap genes, hunchback (hb), Krüppel (Kr), giant (gt), and knirps (kni), and takes into account the expression of the genes bicoid (bcd), caudal (cad), and tailless (tll) as external inputs. Parameter values in the model were calculated by fitting solutions to mean time dependent expression levels from zygotic and maternal/zygotic genes and bcd expression from a single embryo. The resulting model, when run with Bcd gradients from many individual embryos, correctly predicted the variance in position of six gap gene borders. Analysis of model behavior at the numerical level showed that the observed variance was a consequence of gap gene cross regulation [3]. To further elucidate the general nature of the mechanisms controlling variance, we turned off diffusion and analyzed the circuit in individual nuclei using ideas from dynamical systems theory.
This analysis showed that the observed reduction in variation of gap gene expression patterns is a consequence of the action of robust attracting states [4]. The formation of borders of gap gene expression domains could be understood in terms of three qualitative dynamical mechanisms: (1) The movement of attractors; (2) Selection of attractors; (3) Selection of states on a one dimensional attracting manifold. The last of the three mechanisms also causes the domain shifts of the gap genes.
There were two limitations to this analysis. First, the dynamical interactions underlying the observed reduction in variance were elucidated by performing the dynamical analysis on one particular circuit controlled by the median Bcd gradient used for the fit. As a consequence, the variance reduction analysis described in [4] was shifted towards considering only how attractors canalize the variance of initial conditions in the dynamical system, without characterization of Bcd dependence of these attractors. It is possible that other Bcd gradients may entail different dynamical mechanisms of pattern formation. A second limitation is that the individual Bcd gradients used in the analysis had certain systematic scaling errors that exaggerated the variation in threshold location [3,14]. In this paper we extend the dynamical analysis to multiple Bcd gradients, revealing additional dynamical mechanisms which nevertheless work in a coordinated manner to reduce variance. Moreover, we extend the analysis to a system in which the systematic exaggeration of Bcd variance has been removed.
The key idea of this extended analysis is to express gene expression variability in terms of such basic objects of the dynamical systems theory as attractors and attraction basins. The canalization will be explained by specific geometrical properties of these objects. We demonstrate the applicability of this approach using hb border formation as an example. We analyzed how this border forms in terms of the phase portrait of the gap gene dynamical system for various Bcd profiles. We find two specific mechanisms responsible for canalization of the variance, both connected to the phase portrait geometry.
Methods
The ensemble of Bcd concentration profiles
Two sets of spatial profiles of the Bcd concentration were obtained by the numerical processing of raw Bcd data from individual embryos in two different ways. One set was identical to that previously described [3]: the Bcd profiles were retrieved from 89 embryos in the FlyEx database [1517], with the background removed by basic normalization as described [18]. This method as applied to Bcd profiles assumes that the Bcd profile is exponential and that the background profile is quadratic, a point independently supported by staining in null mutants. One embryo was rejected because of a nonexponential profile [3], thus in total 88 profiles were used in the study. For crosschecking purposes, we applied an alternative normalization method to the same raw Bcd data set as described [14]. The method adjusts both concentration scale and levels of constant background to minimize variance in the ensemble. This way of renormalizing data is useful because it is clear from comparison to in vivo work that our profiles exaggerate variance in Bcd amplitude. However, because there is no independent reason to believe that background should be adjusted to minimize variance, the in situ data processed in this manner should be viewed as a lower limit on embryo to embryo variation of the Bcd profile.
The normalized Bcd profiles were approximated by exponential functions v^{Bcd}(x) = A exp(lx) with x varying along the AP axis of the embryo. We obtained in this way a set of 88 (89 for the alternative normalization method) Bcd parameters {A, l} (Additional file 1: Figure S1). The Bcd concentrations in different nuclei at the AP axis were calculated as , where x_{i }was the position of ith nucleus. We selected a "median" Bcd profile by picking an embryo with parameters {A, l} closest to the centroid point in the set of all {A, l} values. This median profile was used to fit the gap gene circuit as described [3].
Additional file 1. Parameters A and l of exponential approximation of individual Bcd profiles for two normalization methods.
Format: PDF Size: 15KB Download file
This file can be viewed with: Adobe Acrobat Reader
The gap gene circuit
We modeled the expression of the network of four gap genes hb, Kr, gt, and kni with the following equations [3,4,13,19]:
where is the concentration of protein encoded by gene a (1 ≤ a ≤ N , N = 4) in nucleus i (1 ≤ i ≤ M(n)) along the AP axis of the embryo. The model incorporates nuclear divisions, number of nuclei M(n) and diffusion coefficients D^{a}(n) depend on the cleavage cycle number n. The function has a sigmoidal graph and describes regulated scaling for the maximal rate R^{a }synthesis of the ath protein. The argument u of g contains inputs from various transcriptional regulators. The first input is a linear combination of all from the network with parameters T^{ab}, which thus quantify the regulatory interactions between the genes. The Bcd concentration profile is a time invariant maternal input to the system. The term is a given timedependent external input from the transcription factors Caudal (Cad) and Tailless (Tll). The function χ(t) equals one during interphase and zero during mitosis, accounting for the fact that synthesis shuts down during this period. Coefficient equals zero if i = 1 (i = M) and one otherwise, preventing protein diffusion outside of the spatial domain. Coefficient λ^{a }is the rate of protein degradation, and constant h^{a }adjusts the threshold of the regulation function.
The initial conditions in the model consist of a spatial gradient of maternally expressed Hunchback (Hb) protein concentration and zero concentrations of the other three proteins. These conditions correspond to the state of the gap gene system at the start of cleavage cycle 13. The solutions in the model are biologically meaningful until the end of cleavage cycle 14A, at which time the midblastula transition occurs and many properties of the embryo change. The nuclei range along the AP axis within the spatial domain from 35% to 92% of the embryo length (EL) where the gap genes express [3,4,13]. This spatial domain includes 30 nuclei in cleavage cycle 13 and 58 in cycle 14A.
All gene expression levels are on a scale of 0255 chosen to maximize dynamic range in the experimental data without saturation. These levels are called relative concentration units throughout the paper. Time dependent inputs for and were obtained from the data by averaging Cad and Tll expression patterns over individual embryos at various time points, as described in details in earlier work [2]. The background was preliminary removed from the individual Cad and Tll patterns as described [18]. The parameters for the model with Bcd normalized as described [18] were those reported in (See Supplementary Material in [3]). For the alternative normalization procedure, we selected a new median Bcd profile from the renormalized data and fitted the model with that profile to the same time dependent averaged gap gene data used previously by either serial or parallel Lam simulated annealing [11,20].
Simplified equations
We considered a simplified version of the model without diffusion, called the "shorted model" as in previous work [4], keeping the parameter values equal to those found by the fitting procedure in the full model (1). Neglect of the diffusion term in (1) decouples nuclei and reduces 4 × M model equations to M independent systems of four equations, each system corresponding to one nucleus. We investigated only a restricted region of the AP axis surrounding the posterior border of the anterior hb expression domain. The region ranges from 37% to 57% EL, about 11 nuclei in width in cleavage cycle 13 and 22 in cleavage cycle 14A. Tll does not act in this region and can be omitted [4]. We can then analyze shorted equations [4] given by
where index x indicates the parametric dependence of solutions. It stems from the dependence of Bcd and Cad protein concentrations on spatial position as determined from experimental data. Such dependence can come directly from experimental data in nuclei, for which x = i, or as a continuous realvalued interpolation x.
The earliest observable indications of the midblastula transition, which include decay of the Bcd gradient, are evident after time class 6 (T6) of cycle 14A, about ten minutes prior to the onset of gastrulation [2,4]. For this reason, we represent the time dependence of v^{Cad}(x, t) from the beginning of cleavage cycle 13 to the end of time class T6. Thereafter we take v^{Cad}(x, t) = v^{Cad}(x, T6) and χ(t) = 1 for times t >T6. Thus, the shorted equations (2) are nonautonomous until T6 and autonomous later.
Attracting sets, basins of attraction, and hb border positions
We analyzed the dynamical system (2) for spatial position x from 37%57%EL region in two mutually complementary directions. We performed a bifurcational analysis in the autonomous version of the system (2), in which χ(t) ≡ 1 and v^{Cad}(x, t) was replaced by value v^{Cad}(x, T6). The bifurcation structure in the system was studied on the BcdCad plane, which is the plane with coordinates v^{Bcd }and v^{Cad }(with v^{Cad }corresponding to the Cad concentration at t = T6), by means of the AUTO package [21]. In this way, all equilibria and domains of their existence on the BcdCad plane were calculated and all bifurcations separating these domains were elucidated. From another direction, at each spatial position of eleven nuclei in cycle 13 and for each of the 88 Bcd profiles (89 for the alternative normalization method), we calculated basins of attraction for each point attractor. This was done by evaluating the equations until late times when the solution is stabilized, with 10 000 random initial conditions uniformly distributed in the biologically relevant subspace of the initial conditions Ω = {0 ≤ v^{Hb }≤ 100, v^{Kr }= v^{Gt }= v^{Kni }= 0} (Ω is a part of the Hb axis in the 4D phase space of the dynamical system). Each point attractor is the asymptotic limit at late times of dynamics starting from the initial conditions which are grouped in a certain part of Ω, and we call this part the basin of attraction for the attractor. To get a more spatially refined picture, we also performed a calculation of attraction basins at thirty spatial positions in the range 37%57%EL for the Bcd ensemble normalized by the basic method.
For each Bcd profile and fixed spatial position x, the attraction basin of each attractor can be represented as an interval (c_{1}(x), c_{2}(x)) on the Hb axis, as all other protein concentrations are zero in Ω. The values v^{Hb }= c_{1 }and v^{Hb }= c_{2 }are lower and upper boundaries for the basin, respectively. We used linear spatial interpolation for these values in order to study the attraction basin boundaries as continuous functions of x for each Bcd profile. Some attraction basins consist of more than one disjoint interval on the Hb axis. In this case, the basin boundaries comprise the boundaries of each connected part of the basin.
We calculated approximations for 1D unstable manifolds of saddles S having a single eigenvalue with positive real part by solving the simplified model equations from two initial conditions S ± w, where w is a scaled eigenvector (corresponding to the unstable eigenvalue) of the Jacobian at the saddle.
In order to calculate the hb border positions in solutions exhibited at the onset of gastrulation, we used spatial interpolation of first and third order for solutions in the shorted and full models, respectively. The orders are different for the two models since the solution in the shorted model is less smooth in space than in the full one because of the absent diffusion term, and the use of high order interpolation schemes could lead to artifacts. As a consequence, we computed the border positions in different ways for the two models. The hb border position in the spatially interpolated solution of the shorted model was defined as the point in the spatial domain at which Hb concentration reached its halfmaximal value, and in the full model as the local inflection point.
Results
We extend the dynamical analysis reported previously for the median Bcd profile [4] to the entire ensemble of Bcd profiles first introduced in [3]. We placed exponential approximations of Bcd profiles from individual embryos in the model (2) with the parameter values calculated by Manu et al. (See Tables S1S2 of Supplementary Information in [3]) and investigated characteristic features of the phase portraits in the model at various spatial positions.
In visualizing the results of our analysis, it is useful to consider two equivalent representations of the spatial information coded in Eqs. (2), either an explicit parameter x denoting AP positions in the range 37%57% EL or a point on the BcdCad plane (Figure 1). We further use the BcdCad plane to show the bifurcation diagrams and discrete spatial positions for presenting attraction basins.
Figure 1. Two representations of the spatial information in the model. (A): The first representation is provided by the explicit use of the spatial position x at the AP axis of the embryo. The panel shows the Cad spatial profile from time class 6 of cleavage cycle 14A (red line) and exponentially approximated individual Bcd profiles in the ensemble (black lines). The vertical dashed lines indicate positions of eleven nuclei in cleavage cycle 13 and in the given spatial range. (B): The second representation is provided by the use of 2D parameter (v^{Bcd}, v^{Cad}), specifying values of Bcd and Cad concentrations in Eqs. (2). This parameter defines a point on the BcdCad plane. The dots in the panel represent the points whose Cad components come from the intersection points between the dashed lines and the Cad profile in (A), and the Bcd components from the intersections of dashed lines with the Bcd profiles in (A). Therefore, these dots describe the actual values of the external input (v^{Bcd}, v^{Cad}) in Eqs. (2) at late times and for the selected eleven spatial positions. The solid line in the panel is the curve parameterized by x from the spatial range, where is the median Bcd profile and v^{Cad}(x) is the Cad profile from (A), showing how the variation of spatial position x in the model with the fixed Bcd profile is read on the BcdCad plane.
For the median Bcd profile, the gap gene expression patterns generated by Eqs. (2) from 35% to 71% EL occur in the same order and locations as those generated by the full model equations (1). The only exception is that the borders become very sharp and domains tend to be mutually exclusive (See Figure 2 in [4]). Over the full ensemble of Bcd profiles, we found three classes of behavior associated with the qualitatively different expression patterns of genes Kr and gt in the anterior vicinity (Figure 2). These classes were visible, but blurred by diffusion, in our previous study of the numerical behavior of the full model (See Figure 3A in [3]).
Figure 2. Three classes of solutions. Examples of solutions at time class 6 (t = T6) in the full (dashed lines) and simplified (solid lines) models for three different Bcd profiles corresponding to three classes IIII (AC, correspondingly). The simplified model equations were obtained from the full model equations by neglecting the diffusion term and the influence from Tll (see Methods). (A): In class I, all borders are present and in the correct order, as is the case with the median Bcd profile. (B): In class II, the anterior Kr border is missing in the full model but present in the simplified model, typically in association with ectopic anterior expression of Kr in place of gt, possibly with some ectopic gt expression. (C): In class III, the anterior Kr border and anterior gt domains are absent in both the full and simplified models.
Attracting sets and various mechanisms of border formation
For each Bcd gradient, we performed an analysis of how border formation was driven by dynamical attractors. Because different Bcd profiles lie on different portions of the BcdCad plane, we first characterized which combinations of attractors are present in different parts of this plane by performing a bifurcational analysis (Figure 3) of the shorted model (2). We then calculated the basins of attraction for all attractors at a discrete set of eleven positions corresponding to cycle 13 nuclei in the BcdCad plane with initial conditions varied in Ω = {0 ≤ v^{Hb }≤ 100, v^{Kr }= v^{Gt }= v^{Kni }= 0}.
There are four stationary attractors (A_{1}A_{4}) in the portion of the BcdCad plane shown in Figure 3. The attractors can be coded with quadruples consisting of 0, X, or 1 for each attractor component by inspecting whether the corresponding protein concentrations have small, intermediate, or large values at the attractors (Table 1; see the existence domains for attractors corresponding to the alternative normalization method in Additional file 2: Figure S2; examples of spatial dependence of attractors are shown in Additional files 3 and 4: Figures S3 and S4 for the two normalization methods). For example, A_{3 }= 0100 means that Kr is highly expressed at this attractor with the other genes staying repressed. A_{1 }is the only attractor that continuously changes its code with varying Bcd and Cad concentrations. There are six domains on the BcdCad plane in which various combinations of A_{1}A_{4 }exist (Figure 3). Attractor A_{1} exists inside the entire portion of the considered plane, while the other attractors are involved in the bifurcations at the borders of their existence domains. There are three types of bifurcations found affecting attractors: saddlenode, Hopf, and BogdanovTakens bifurcations (see more detailed description of all bifurcations in Additional file 5: Protocol S1 and Additional file 6: Figure S5). A limit cycle appears at the Hopf bifurcation, which is a nonstationary oscillating attracting state. However, we have not found any oscillatory attractors when the initial conditions are varied in Ω, and, therefore, we exclude this type of attractor from the analysis. The gray dots in Figure 3 show the late time values of the external input in Eqs. (2) used in the calculations of attraction basins for eleven nuclei positions and all of the Bcd profiles. The location of the grey dots on the bifurcation diagram indicates the existence of specific attractors in the phase space for a given nucleus. However, not all of these attractors can be reached under the given biological initial conditions in Ω, as attraction basins of some attractors are disjoint with Ω (Additional file 7: Figure S6).
Figure 3. The bifurcation diagram for the simplified model on the BcdCad plane. The abscissa axis shows v^{Bcd }values in the equations, and the ordinate axis is for v^{Cad }values corresponding to the Cad concentrations at times after T6, when the model equations are autonomous. The colored regions are domains of existence for point attractors A_{1}A_{4}, described symbolically in Table 1. The boundaries (black solid lines) delimiting these domains represent the positions of all bifurcations affecting the attractors. The black dashed lines show the positions of bifurcations affecting only saddle equilibria. The white dashed and solid lines are loci of points where the Gt concentration at attractor A_{1 }equals 50 and 150, respectively. Therefore, to the left of the white dashed line, to the right of the white solid line, and between the lines. The gray dots are the same as in Fig. 1B.
Table 1. The symbolic codes for all attractors in the study calculated for the Bcd and Cad ranges from Figure 3.
Additional file 2. The bifurcation diagram for the new parameter values.
Format: PDF Size: 158KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 3. The spatial dependence of attractors in the shorted model for the median Bcd profile.
Format: PDF Size: 17KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 4. The spatial dependence of attractors in the shorted model with the new parameter values and for the median Bcd profile.
Format: PDF Size: 18KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 5. Detailed description of the bifurcations in the model (the text contains reference to Figure S5).
Format: PDF Size: 23KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 6. The bifurcation diagram with more details.
Format: PDF Size: 29KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 7. The existence domains on the BcdCad plane for attractors A_{1}A_{4 }following from calculations at discrete spatial positions.
Format: PDF Size: 153KB Download file
This file can be viewed with: Adobe Acrobat Reader
As shown elsewhere [4], the dynamics in the model for the median Bcd profile is qualitatively different in parts of the AP axis which are anterior and posterior to the position of the bifurcation annihilating A_{3}^{. }This position appears as a diagonal line running from the upper left to lower right of Figure 3. The anterior and posterior dynamical regimes are characterized by different types of attracting sets governing the solution dynamics. Solutions at the end of cycle 14A (t = τ) are very close to point attractors in the anterior regime and to attracting manifolds in the posterior one (examples of these attracting sets are shown in Figure 4). We found that these dynamical regimes are preserved across all the individual Bcd profiles.
Figure 4. Dynamical mechanisms of hb border formation for various Bcd profiles. Each panel shows the HbKrKni projections of the 4D phase portrait fragments from the hbexpressing and hbnonexpressing border nuclei simultaneously. The hbexpressing and hbnonexpressing border nuclei are the ones placed just anterior and just posterior to the hb border position, respectively. The axes labels for all panels are as for (A) and show the corresponding protein concentrations. The blue balls are the attractors from the hbexpressing border nucleus phase portraits, and the brown ones from the hbnonexpressing nucleus portraits. The red balls are the saddles from the nuclei specified below for each panel. The cyan tubes are branches of the unstable manifolds for the hbexpressing border nuclei saddles, and the magenta ones for the hbnonexpressing nuclei saddles. The green trajectories are the solutions for the biological initial conditions from the hbexpressing and hbnonexpressing border nuclei. The initial conditions consist of the maternal Hb concentration and are zero for all other proteins. The trajectories gradually turn to red as time approaches the end of cleavage cycle 14A (t = τ) and become blue for later times. (A): The picture is for Bcd profile #1, with the hbexpressing border nucleus at 49%EL and the hbnonexpressing border nucleus at 51%EL . (B): Bcd profile #71; the border nuclei at 49%EL and 51%EL ; the red saddle is from the hbnonexpressing border nucleus. (C): Bcd profile #6; the border nuclei at 47%EL and 49%EL ; the red saddle is from the hbexpressing border nucleus. (D): Bcd profile #37; the border nuclei at 51%EL and 53%EL ; the red saddle is from the hbexpressing border nucleus.
To study the dynamical mechanisms of hb border formation, we examined the phase portraits in the shorted model at spatial positions on either side of the border. The mechanism determining the border was found by inspecting which attracting sets govern the solution dynamics from the biological initial conditions in the two nuclei, a hbexpressing border nucleus placed just anterior to the hb border position and a hbnonexpressing border nucleus just posterior to that. Border formation is interpreted in these terms as a switch of solution between attracting states in these two nuclei.
We found four qualitatively different mechanisms of hb border formation for Bcd profiles from the ensemble, depending on the type of attracting sets approached by the solutions in the two nuclei by t = τ (Figure 4). For example, Figure 4B shows that for corresponding Bcd profile the border forms by a solution switching from attractor A_{2 }= 1100 in the hbexpressing border nucleus to a hbOFF state at the unstable manifold of one of the saddles in the hbnonexpressing border nucleus. Other panels in the figure can be interpreted in a similar way.
There are 66 Bcd profiles associated with the phase portraits which exhibit the attractorattractor switch mechanism of border formation (Figure 4A) and 20 profiles corresponding to the attractormanifold switch mechanism (Figure 4B), which make these two mechanisms predominant for the Bcd ensemble. The other mechanisms shown in Figure 4C, D are rare and occur for only one Bcd profile each. The distribution of the four mechanisms over the three solution classes from Figure 2 is summarized in Additional file 8: Table S1. For all Bcd profiles, only attractors A_{1}, A_{2}, and A_{3 }participate in hb border formation. A_{4 }contains a nonempty attraction basin in Ω for only a few nuclei which lie in a portion of the BcdCad plane that never contains the Hb border (Additional file 7: Figure S6).
Additional file 8. The distribution of Bcd profiles over the four mechanisms of hb border formation and over solution classes IIII.
Format: PDF Size: 9KB Download file
This file can be viewed with: Adobe Acrobat Reader
Spatial configuration of attraction basins
To further investigate how the border formation mechanisms imply the observed hb border variance, we calculated the spatial and Bcd dependence of attraction basin boundaries. We visualize the way that the biological initial conditions choose a basin in different nuclei and try to connect it to the appearance of the hb border in the solution by the time t = τ. We construct graphs combining the initial Hb profile and the attraction basins for all attractors at different spatial positions, and highlight the range of positional variance for hb. Typical examples of these graphs are shown in Figure 5A,B.
Figure 5. Spatial configuration of attraction basins in the simplified model. (A,B): Basins of attraction at the eleven spatial positions (nuclei in cleavage cycle 13) for the two Bcd profiles from the ensemble (profiles (A) #1 and (B) #71, corresponding to the phase portraits from Fig. 4A and 4B, respectively), together with the initial Hb concentration profile (black curve) and highlighted positions of the hb borders in the model (dashed magenta lines). Each thick vertical line presents the domain Ω, which is the set of initial values 0 < v^{Hb }<100 at the Hb axis, in which the attraction basins for attractors A_{1}A_{3 }are marked by corresponding colors. The basins were calculated and shown at spatial positions ranging from 37% EL to 57% EL with 2% bins, and the width of the vertical lines containing these basins is arbitrarily chosen for better visibility. (C): Intersection of the initial Hb profile (black curve) and the attraction basin boundaries (blue and red curves) for all but one Bcd profiles classified in the two families in (D), as described in the text. The blue and red colors in both (C) and (D) correspond to the first and second Bcd families, respectively. For the first Bcd family, the initial Hb profile crosses a boundary separating either the attraction basins of A_{1 }and A_{3 }(see (A) for an example) or the basins of A_{2 }and A_{3}, and all these boundaries are shown as blue curves in (C) (see Methods for description of basin boundary calculation). For all Bcd profiles from the second family, the basin boundaries shown in (C) as red curves separate the basin of A_{2 }and the lower part of the basin of A_{1} at various spatial positions (see (B) for an example). The two line segments at the bottom of (C) show the ranges of the hb border positions in solutions of the simplified model equations for Bcd profiles from the two families. The inset in (D) shows the correlation between the hb border positions in the model (vertical axis) and the positions of intersection points between the initial Hb profile and the basin boundaries from (C) (horizontal axis), for Bcd profiles from the two families.
For all but one of the Bcd profiles, hb border formation is associated with the transition of the initial Hb profile from the basin of one point attractor to the basin of another point attractor. This can be seen by comparing the attraction basins of the initial Hb concentration values in the two nuclei surrounding the Hb boundary position. For example, the transition in Figure 5A is from the basin of A_{1 }in the hbexpressing border nucleus to the basin of A_{3 }in the hbnonexpressing border nucleus (the A_{1 }→ A_{3 }transition); the transition in Figure 5B is from the basin of A_{2 }to the basin of A_{1 }(A_{2 }→ A_{1 }transition). The basintobasin transitions can also be detected in Figure 4, where the solutions in the two nuclei go to different point attractors. Bcd profile #54 provides the only exception from this basintobasin transition rule. For this Bcd profile, the solution trajectories in both hbexpressing and hbnonexpressing border nuclei eventually end up at the same point attractor A_{2}. The hb border forms in this case according to the mechanism from Figure 4B (the attractormanifold switch) with the only difference that the solution trajectory in the hbnonexpressing border nucleus is attracted to a hbOFF state on an unstable manifold which connects a saddle and attractor A_{2 }(data not shown). This Bcd profile was excluded from the further analysis. For all cases with the A_{1 }→ A_{3 }transition, attractor A_{1 }is in the state for the hbexpressing border nucleus. The A_{1 }→ A_{3 }and A_{2 }→ A_{3 }transitions are similar in the sense that they both correspond to the switch from an attractor with hbON state to an attractor with hbOFF state. We found 69 Bcd profiles in the ensemble associated with either of these two transitions. We call these profiles Family I in what follows. These transitions can hypothetically correspond to any of the four mechanisms possible for the hb border formation, because the attracting invariant manifolds may also participate in the border formation in this case (see, e.g., Figure 4C). On the other hand, the A_{2 }→ A_{1 }transition from Figure 5B describes the transition of the initial Hb concentration between the basins of attractors which both have the hb component in an 'ON' state. The hb border forms in this case by the attracting manifold which provides the necessary hbOFF state at t = τ (see Figure 4B,D). We found a total of 18 Bcd profiles leading to the A_{2 }→ A_{1 }transition. We refer to these profiles as Family II. Among these profiles, 16 are associated with the state and 2 profiles with the state for the hbnonexpressing border nucleus. The basin of A_{3 }in Figure 5A stops to exist for nuclei located posterior to 51%EL. The system is in the basin of A_{1 }for these positions and approaches a hbOFF state by t = τ with the help of an attracting manifold. This participation of an attracting manifold in attaining a hbOFF state is the distinctive feature of the aforementioned posterior dynamical regime.
Table 2 shows how the types of initial Hb transitions between basins are distributed with respect to the border formation mechanisms described above.
Table 2. The distribution of Bcd profiles over the four mechanisms of hb border formation and two families described in the text.
Variability of basin boundaries and hb border positions
The fact that the hb border forms according to the basintobasin transition of the initial Hb profile implies that this profile crosses a boundary between two basins. Therefore, we can relate the border position to the position of the intersection point between the initial Hb profile and the spatial profile of the basin boundary (see Methods for a description of how this profile is calculated). One can speculate that these two positions are the manifestation of two levels of positional information readout in the gap gene circuit. In this subsection, we compare the positional variability for hb predicted by the model for these two levels. We calculated the spatial profiles of the basin boundaries for all Bcd gradients. The ranges of the border positions corresponding to the two Bcd families are shifted with respect to each other, and the same shift is observed at the level of the basin boundaries (Figure 5C,D). As expected, the positions of the hb border and basin boundary intersection points exhibit strong correlation for both Bcd families (see the inset in Figure 5D). Note that the steplike form of the correlation curve in the figure follows from two factors, the discrete spatial positions used in calculations and the absence of diffusion in the model equations (i.e., the solution of the simplified model more sharply depends on the spatial position).
This correlation means that the intersections between the initial Hb profile and specific basin boundaries encode the hb border positions. As shown in Table 3, the canalization of Bcd variation takes place both at the level of the basin boundaries and at the level of hb border in the simplified model. The variabilities of the intersection points and the hb border positions are approximately the same for each Bcd family and between the families. On the other hand, the Bcd positional variance is significantly different in the two families.
Table 3. Positional variance for Hb and Bcd.
A nonlinear curve of hb border response to Bcd variation
The Bcd dependence of hb border can also be investigated by inspecting the curve representing the response of hb border position to Bcd concentration levels in a vicinity of this position. We studied this response prescribed by the model for Bcd profiles from the ensemble (Figure 6).
Figure 6. The nonlinear curve of hb border response to Bcd variation in (A) the simplified and (BE) full models. (A): The value of Bcd concentration at the AP position where the initial Hb profile intersects the corresponding basin boundary from Fig. 5C versus this position, for all but one Bcd profiles from the ensemble (Bcd profile #54 was excluded from the analysis as explained in the text). The colors mark the three solution clusters IIII from Fig. 2: green corresponds to class I, blue to II, and red to III. The yellow stars indicates the median Bcd case. The inset in (A) shows the Bcd dependence of the attraction basin boundary for the attractor A_{3 }at fixed AP position equal to 47% EL. As explained in Methods, the attraction basin has the form of interval (c_{1},c_{2}) at the Hb axis; the Bcd concentration at 47% EL versus the value [Hb] = c_{2} at the same position is plotted in the inset for all Bcd profiles from the first family. (B): The same response curve as in (A) but obtained from calculations in the full model (1) and with the hb border positions used instead of the intersection point positions. The inset shows all Bcd profiles used in the calculations. Colors in (B) correspond to those in (A). (C): The response curve expressed in terms of 'masses' of the Bcd patterns and Hb solutions. The mass is calculated as the integral of a pattern (a Bcd profile or a Hb solution at the end of cleavage cycle 14A) over the spatial region 3564%EL. The Hb pattern mass positively correlates with the hb border position (data not shown). The dots represent results of calculations in the full model, while the triangles correspond to calculations with the same equations but with the gap gene cross regulation of hb turned off (T^{hb←a }set to zero in (1) for each gap gene a). The green, blue, and red symbols indicate cases with a small, intermediate, and large mass of the Bcd patterns, respectively. The Bcd and Hb patterns used in calculations of the colored dots are shown in (D) with respective colors, and the same for triangles in (E). The Bcd patterns are dashed lines, and the Hb solutions are solid ones.
The response curve exhibits a specific nonlinear form which is the same in the simplified model (Figure 6A) and in the full one (Figure 6BC). It persists even if we recalculate the curve by using the normalized individual Bcd profiles instead of their exponential approximations (Additional file 9: Figure S7), indicating that the results are independent of the nuclear noise present in the nonapproximated Bcd data.
Additional file 9. The response curve for the normalized individual Bcd profiles instead of their exponential approximations.
Format: PDF Size: 14KB Download file
This file can be viewed with: Adobe Acrobat Reader
The nonlinear form of the curve is the result of the reaction to Bcd variation from the whole gap gene network. As a consequence, we have a linearly increasing response curve when gap gene cross regulation of hb is removed from the full model equations (Figure 6C). In that case, the gradually increasing amplitude (or a mass) of Bcd profile across the ensemble results in gradual posterior shift of the hb border (Figure 6E). This happens because Bcd is an activator of hb and the spatial profiles of both Bcd and Hb concentrations are monotonously decreasing in the vicinity of hb border position. In contrast, there can be reversed (anterior) shifts of the hb border for some range of Bcd concentration in the presence of gap gene cross regulation (Figure 6D).
The inset in Figure 6A shows an example for one spatial position of how the attraction basin boundary depends on Bcd concentration. This dependence locally represents the response curve at the level of attraction basins in the model. The specific form of the curve in the inset underlies the nonlinear form of the response curve for the intersection points and border positions.
Results for the alternative Bcd normalization method
As the Bcd variance can be exaggerated in the data normalized by the basic method, we considered the same ensemble of Bcd profiles but normalized by the alternative artificial normalization method, which minimizes variance in the ensemble (see Methods). We studied the model on the newly normalized Bcd profiles in order to crosscheck our results. The corresponding parameters A and l of exponential approximation for the Bcd profiles are shown in Additional file 1: Figure S1. We calculated a new set of parameter values (Additional file 10: Table S2) in the full model (1) with new median Bcd profile by fitting to the same gap gene expression data that were used previously [3,4]. Main qualitative features of these parameter values and corresponding solution of the simplified equations (2) (Additional file 11: Figure S8) are very close to those in the above study [3,4].
Additional file 10. New parameter values in the model obtained by optimization with a median Bcd profile from the Bcd data normalized by the alternative method.
Format: PDF Size: 15KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 11. The solutions of the full model equations and their simplified version for the new parameter values.
Format: PDF Size: 18KB Download file
This file can be viewed with: Adobe Acrobat Reader
Due to the specific statement of the alternative normalization method, the Bcd positional variance in the ensemble became much lower: the standard deviation σ = 1.2%EL. The corresponding positional variance for hb in the full model equals to σ = 0.5%EL.
We calculated a new bifurcation diagram on the BcdCad plane for the simplified model (Additional files 2, 12, 13: Figures S2, S9, and Protocol S2). There are two new attractors, A_{5 }= 0101 and A_{6 }= 1000, in addition to attractors A_{1}A_{4 }which are analogs of the corresponding attractors to the above study. Attractor A_{6 }exists in a very narrow domain on the BcdCad plane, where it coexists with attractor A_{1 }at the state (Additional file 2: Figure S2). Thus, the attractor 1000 is unique in this domain. The calculations of attraction basins at discrete spatial positions did not detect A_{4}, since the newly normalized Bcd concentrations never entered the existence domain for this attractor (Additional file 14: Figure S10). As in the above study, only attractors A_{1}A_{3 }are important for the process of hb border formation. The spatial dependence of these attractors is mostly similar to the previous results, except that the bifurcations changed their positions at the AP axis (see, e.g., the case of the median Bcd profile in Additional file 4: Figure S4, in comparison with previously discussed Additional file 3: Figure S3). Despite the difference in the bifurcation structure, new calculations also revealed the existence of an attracting invariant manifold governing the dynamics in a posterior spatial region, which is the analog of the previously found manifold [4]. This manifold has similar geometrical properties, reflecting the fact that it is responsible for the observed anterior shifts of the posterior gap domains taking place in cycle 14A [4].
Additional file 12. The bifurcation diagram for the new parameter values with more details.
Format: PDF Size: 19KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 13. Detailed description of the bifurcations for the new parameter values (the text contains reference to Figure S9).
Format: PDF Size: 14KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 14. The existence domains on the BcdCad plane for attractors A_{1}A_{6 }in the model with the new parameter values following from calculations at eleven spatial positions.
Format: PDF Size: 58KB Download file
This file can be viewed with: Adobe Acrobat Reader
Three of the four mechanisms shown in Figure 4 exist for the new parameter values and for various newly normalized Bcd profiles: the attractorattractor (Figure 4A), manifoldattractor (Figure 4C), and manifoldmanifold (Figure 4D) switches. Analyzing how the initial Hb profile chooses the attraction basin in the nuclei neighboring the hb border, we detected 85 Bcd profiles leading to the A_{1 }→ A_{3 }transition, 2 Bcd profiles leading to the A_{2 }→ A_{3 }transition, and 2 profiles which were not associated with the change of basins across these nuclei (they corresponded to the switch between different states on the unstable manifold of a saddle) (see Additional file 15: Table S3 for a distribution of these cases with respect to the solution classes IIII). The model demonstrates a picture of the initial Hb profile crossing the attraction basin boundaries similar to the first family case in Figure 5C, but with a smaller variance (Additional file 16: Figure S11). The Bcd response curve for the hb border in the model with the new parameter values is also not linear (Additional file 17: Figure S12).
Additional file 15. Classification results for the Bcd profiles in the case of the alternative normalization method.
Format: PDF Size: 13KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 16. The spatial configuration of attraction basins in the model with the new parameter values.
Format: PDF Size: 24KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 17. The response curve for the new parameter values.
Format: PDF Size: 63KB Download file
This file can be viewed with: Adobe Acrobat Reader
Discussion
Mechanisms of border formation
We presented the dynamical analysis of the simplified model of the gap gene network on the ensemble of early Drosophila embryos. The main goal was to decode the mechanistic basis of the gap gene border formation and stability under the Bcd morphogen variance. The hb border formation mechanisms were described in terms of attracting sets and their attraction basins calculated in the nuclei surrounding the border position.
The results reveal that the border formation can be associated with the event of intersection between a boundary separating the attraction basins of two different point attractors and the initial Hb profile presenting the input from the maternally expressed hb gene. Attracting sets of another type, the unstable manifolds of saddle equilibria, actively participate in the adjustment of the border position. They do so by attracting the solution trajectories in the nuclei surrounding this position. The model predicts that these attracting manifolds can be involved in the border formation for some Bcd profiles.
The hb border correctly forms in the model by the onset of gastrulation for all individual Bcd profiles. For about a half of these profiles, however, the Kr and Gt patterns in the solutions exhibit defects in the anterior part of the spatial domain (solution classes II and III). It turns out that the hb border formation mechanism involving the attracting manifolds is mostly associated with these cases. This may lead to the conclusion about restricted applicability of this mechanism in the case of hb expression. However, this mechanism exists and plays an important role for the gap domain borders in a posterior part of the embryo, where the domains form and vary in time under the control of an unstable manifold [4]. To analyze canalization for the posterior borders, the variation for external inputs from Cad and Tll should be taken into account, where these transcription factors are among the key regulators, and a modified model should be considered including an input from the terminal gene huckebein [22].
Mechanisms of canalization
As previously reported, the model exhibits a significant filtration (canalization) of the Bcd positional variability at the level of hb border formation [24]. Our results show how this filtration stems from the stable behavior of the attraction basin boundaries. It was shown in [3,10] that the mutual regulatory repression between the gap genes accounts for the observed variance reduction, thus presenting a buffering mechanism for canalization [23,24]. We translated this buffering mechanism to the level of attractors and their attraction basins. As the hb border position is well encoded by the intersection between the initial Hb profile and corresponding attraction basin boundaries, the stability of hb border predicted by the model can be explained by inspecting the geometrical properties of these attraction basins.
From this inspection, we can elucidate the following two mechanisms responsible for the observed robustness. First, the initial Hb profile is a monotonously decreasing function of AP position, while the basin boundary to be crossed is a monotonously increasing one (Figure 5), i.e., these curves have opposite dependencies on the AP position. This purely geometrical fact evidently prescribes a smaller variation of the intersection point when the basin boundary is changing due to the variance of Bcd concentration, as opposed to the case if the curves would jointly rise or jointly fall along the AP axis (we illustrated this mechanism of canalization in Additional file 18: Figure S13).
Additional file 18. Schematic illustration of the first canalization mechanism.
Format: PDF Size: 16KB Download file
This file can be viewed with: Adobe Acrobat Reader
The second mechanism is associated with the specific nonlinear form of the response curve from Figure 6. The gap gene cross regulation of hb bends the response line exhibited in absence of this regulation (Figure 6C). This bending effectively reduces the Hb positional variance by about half. In terms of attractors, this bending is controlled by the fact that a basin boundary responsible for the hb border formation does not change monotonously, but oscillate in the state space with the changing Bcd profile.
The results show that the full range of the hb positional variance is broken down into two almost equal parts, the anterior and posterior ones (see the line segment in Figure 5C). These parts are associated with two families of the Bcd individual profiles (Family I and Family II, respectively) and two different mechanisms of hb border formation. The Bcd profiles from Family I lead to the hb border formation as a switch from a hbON attractor in a hbexpressing nucleus to a hbOFF attractor in a hbnonexpressing nucleus, while for Family II the border forms with the help of an attracting invariant manifold in a hbnonexpressing nucleus. Since the difference between the two families is in the amplitude of the Bcd profiles, we conclude that Bcd profiles of high amplitude canalize by a dynamical mechanism different from those of lower amplitude. Each dynamical mechanism provides only half of the full variance for the hb border, but in two adjacent spatial domains. Therefore, the change of the dynamical mechanism that happens with rising Bcd amplitude effectively doubles the variance.
The hb border positions from the more posterior range are placed posterior to the spatial position of a bifurcation annihilating attractor A_{3}. This bifurcation position delimits the anterior and posterior dynamical regimes in the model, as described previously [4]. Therefore, the Bcd profiles from the second family shift the hb border to the posterior dynamical regime, which is characterized by an active role of an attracting invariant manifold in the pattern formation.
The results indicate that the posterior range of hb positional variation is almost equal to the anterior one only due to smaller variation of the Bcd profiles in Family II compared to Family I. This suggests that the solutions in the anterior and posterior dynamical regimes have quite different sensitivity rates to variation of the Bcd concentration. For Family I, the standard deviation for the hb border position is 2.6 times less than for the Bcd threshold position, while it is only 1.4 times less in the case of Family II. This difference can be explained by an observation that Bcd profiles of higher amplitude correspond to the linear part of the response curve from Figure 6, and this is a consequence of specific regulatory interactions in the gap gene circuit as explained further.
We have used the model (1) to study the canalization mechanisms based on the assessment that the model provides one of the best spatiotemporal precision for the description of gap gene expression [25]. This model is an approximation to a more general model of gene regulation, which should be grounded on the statisticalmechanical formalism. One possible limitation is the linear approximation for the argument of the nonlinear regulation function g. The canalization mechanisms described in terms of attractors and attraction basins generally depend on the structure of the model that predicts these attracting states. Therefore, an important direction for future investigations should be verification of the proposed mechanisms in a phase space of a more general model.
Response curve
The nonlinear nature of the Bcd readout by the gap gene circuit is clearly represented in a specific nonlinear form of the response curve showing the Bcd dependence of the hb border position in the model. The nonlinear part of the curve can be explained by the regulatory actions on hb from the other gap genes. In particular, a regulatory analysis in the full model revealed that the regulatory interactions between hb, gt, and Kr underlie the folding part of the response curve (Additional file 19: Protocol S3). The gap gene crossregulation also participate in the linear parts of the response curve by tuning the incline of these parts.
Additional file 19. The regulatory analysis of the response curve (the Protocol contains Figures S14S17).
Format: PDF Size: 81KB Download file
This file can be viewed with: Adobe Acrobat Reader
It was previously pointed out that the gt and Kr expression borders in the anterior part of the AP axis show large variation in the model in response to Bcd variation because the model is missing some regulators in this part [3]. For example, these gt and Kr borders are absent in the solutions from class III. This fact raises doubts on the specific folding part that the response curve exhibits in the middle range of the Bcd concentration values. On the other hand, the folding part exists only for the Bcd profiles associated with the solutions from class I, with all expression borders formed correctly, which means that an essential portion of the artificial variation of the gt and Kr borders can be excluded from the consideration without affecting the folding form of the curve.
Comparison of results for two normalization methods
We investigated the model on the ensemble of Bcd profiles normalized by the alternative method, which provided lower Bcd variance [14]. One used this method as an artificial limit case, in which we dealt with the ensemble possessing minimal Bcd variance, and applied it for the crosschecking purposes.
We have not found any essential discrepancy in the mechanisms of hb border formation and canalization for the two normalization methods. A distinct bifurcation structure in the model with the new parameter values does not lead to changes in the solutions during the biologically important time. The model preserves an attracting invariant manifold related to the posterior dynamical regime. The same border formation mechanisms appear except the one associated with the attractormanifold transition. It is important that, even though the second family of Bcd profiles does not appear in the alternative normalization case, the invariant manifolds still play their role in adjusting the border position. The model also demonstrates an essentially nonlinear response curve for the hb border. Therefore, our conclusions formulated above are robust with respect to the choice of the normalization method, and, in more general terms, they should be valid for different estimates of the actual Bcd variance.
This correspondence can be explained by the fact that the parameters A and l obtained for the alternatively normalized Bcd profiles form a subset in similar parameters obtained in the case of the basic normalization method (see Additional file 1: Figure S1). Roughly speaking, we can associate the alternatively normalized Bcd profiles with Family I. In particular, this means that the Bcd data rescaled according to the alternative algorithm support the conclusion formulated above about different dynamical mechanisms of canalization for Bcd profiles of different amplitude.
There is an important issue concerning the comparison of the Bcd variance filtration rates. The calculations reveal that, for the basic normalization method, the Hb positional variation of 1.3%EL in the model output follows from the Bcd positional variation of 4.5%EL, thus implying that more than 70% of the positional variance has been filtrated. The same calculations for the alternative normalization method give the filtration rate of approximately 60%. Therefore, the filtration still happens in the model even if we normalize Bcd profiles according to the precisionist hypothesis [14]. This result is quite expected since the reported dynamical mechanisms underlying the processing of the Bcd variation in the model are valid irrespective of the absolute variation range. Whatever actual variation the Bcd morphogen exhibits, the nonlinear model response translates it to a smaller variation of the target gene patterns.
Conclusions
The formation of hb border is coded by the intersection between the maternal Hb gradient and a boundary between attraction basins in the gap gene dynamical system. Small positional variance for hb border can be explained by the geometrical properties of this basin boundary and its nonmonotonic dependence on the Bcd concentration. Main features of the phase portraits underlying the canalization mechanisms do not depend on the normalization method for Bcd.
Authors' contributions
VVG, M, and JR conceived and designed the course of the study. VVG performed calculations of attractors and attraction basins, model optimization, and analysis of data and results. LP carried out the bifurcational analysis and calculated the phase portraits. EMM applied alternative normalization to the Bcd data and analyzed the data. MGS and JR performed biological analysis of the results. VVG and LP made the figures. AMS, MGS, and JR supervised the work. VVG, LP, JR, MGS, and AMS wrote the paper. All authors read and approved the final manuscript.
Acknowledgements
The study was supported by NIH Grant RR07801, CRDF GAP award RUB11578ST05, and RFBR grants 080100315a, 080400712a, 100100627a, 110100573a, and 110401162a. We thank Konstantin N. Kozlov for the help with the parameter optimization code and Johannes Jaeger for valuable feedback on the manuscript. The optimization calculations were partly performed at the Joint Supercomputer Center of the Russian Academy of Sciences.
References

Waddington CH: Canalization of development and the inheritance of acquired characters.
Nature 1942, 150:563565. Publisher Full Text

Surkova S, Kosman D, Kozlov K, Manu , Myasnikova E, Samsonova A, Spirov A, VanarioAlonso CE, Samsonova M, Reinitz J: Characterization of the Drosophila Segment Determination Morphome.
Developmental Biology 2008, 313(2):844862. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Manu , Surkova S, Spirov AV, Gursky V, Janssens H, Kim A, Radulescu O, VanarioAlonso CE, Sharp DH, Samsonova M, Reinitz J: Canalization of gene expression in the Drosophila blastoderm by gap gene cross regulation.
PLoS Biology 2009, 7:e1000049. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Manu , Surkova S, Spirov AV, Gursky V, Janssens H, Kim A, Radulescu O, VanarioAlonso CE, Sharp DH, Samsonova M, Reinitz J: Canalization of gene expression and domain shifts in the Drosophila blastoderm by dynamical attractors.

Thom R: Topological Models in Biology.
Topology 1969, 8:313335. Publisher Full Text

Gibson G, Wagner G: Canalization in evolutionary genetics: a stabilizing theory?
BioEssays 2000, 22:372380. PubMed Abstract  Publisher Full Text

Hansen TF: The Evolution of Genetic Architecture.
The Annual Review of Ecology, Evolution, and Systematics 2006, 37:123157. Publisher Full Text

Meiklejohn CD, Hartl DL: A single mode of canalization.
Trends in Ecology and Evolution 2002, 17(10):468473. Publisher Full Text

Lott S, Kreitman M, Palsson A, Alekseeva E, Ludwig M: Canalization of segmentation and its evolution in Drosophila.
Proceedings of the National Academy of Sciences USA 2007, 104:1092610931. Publisher Full Text

Vakulenko S, Manu , Reinitz J, Radulescu O: Having the right proportions: Interacting interfaces ensure robust spatial patterning.
Physical Review Letters 2009, 103:168102. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Reinitz J, Sharp DH: Mechanism of eve stripe formation.
Mechanisms of Development 1995, 49:133158. PubMed Abstract  Publisher Full Text

Reinitz J, Sharp DH: Gene Circuits and Their Uses. In Integrative Approaches to Molecular Biology. Edited by Collado J, Magasanik B, Smith T. Cambridge, Massachusetts, USA: MIT Press; 1996:253272.

Jaeger J, Surkova S, Blagov M, Janssens H, Kosman D, Kozlov KN, Manu , Myasnikova E, VanarioAlonso CE, Samsonova M, Sharp DH, Reinitz J: Dynamic control of positional information in the early Drosophila embryo.
Nature 2004, 430:368371. PubMed Abstract  Publisher Full Text

Gregor T, Tank DW, Wieschaus EF, Bialek W: Probing the limits to positional information.
Cell 2007, 130:153164. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Poustelnikova E, Pisarev A, Blagov M, Samsonova M, Reinitz J: FlyEx Database. [http://urchin.spbcas.ru/flyex] webcite
2005.

Poustelnikova E, Pisarev A, Blagov M, Samsonova M, Reinitz J: A database for management of gene expression data in situ.
Bioinformatics 2004, 20:22122221. PubMed Abstract  Publisher Full Text

Pisarev A, Poustelnikova E, Samsonova M, Reinitz J: FlyEx, the quantitative atlas on segmentation gene expression at cellular resolution. [http://nar.oxfordjournals.org/cgi/content/abstract/gkn717v1] webcite
Nucleic Acids Research 2008, 37:D560D566. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Myasnikova E, Samsonova M, Kosman D, Reinitz J: Removal of background signal from in situ data on the expression of segmentation genes in Drosophila.
Development, Genes and Evolution 2005, 215:320326. Publisher Full Text

Mjolsness E, Sharp DH, Reinitz J: A connectionist model of development.
The Journal of Theoretical Biology 1991, 152:429453. Publisher Full Text

Chu KW, Deng Y, Reinitz J: Parallel simulated annealing by mixing of states.
The Journal of Computational Physics 1999, 148:646662. Publisher Full Text

Doedel E, et al.: AUTO: Software for continuation and bifurcation problems in ordinary differential equations. [http://indy.cs.concordia.ca/auto] webcite

Ashyraliyev M, Siggens K, Janssens H, Blom J, Akam M, Jaeger J: Gene Circuit Analysis of the Terminal Gap Gene huckebein.
PLOS Computational Biology 2009, 5:e1000548. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Waddington CH: Organisers and Genes. Cambridge, U.K.: Cambridge University Press; 1940.

Waddington CH: The Strategy of Genes. London: George Allen & Unwin; 1957.

Jaeger J: Modelling the Drosophila embryo.
Molecular BioSystems 2009, 5:15491568. PubMed Abstract  Publisher Full Text