Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS One**|**v.5(4); 2010**|**PMC2858077

Formats

Article sections

Authors

Related links

PLoS One. 2010; 5(4): e10275.

Published online 2010 April 21. doi: 10.1371/journal.pone.0010275

PMCID: PMC2858077

Andreas Bergmann, Editor^{}

University of Texas MD Anderson Cancer Center, United States of America

* E-mail: gro.cmhcc@am.nuj

Conceived and designed the experiments: JD JM. Performed the experiments: JD. Analyzed the data: JD WW LJL JM. Contributed reagents/materials/analysis tools: WW LJL. Wrote the paper: JD WW JM.

Received 2010 February 1; Accepted 2010 March 31.

Copyright Deng et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

This article has been cited by other articles in PMC.

Bicoid (Bcd) is a *Drosophila* morphogenetic protein responsible for patterning the anterior structures in embryos. Recent experimental studies have revealed important insights into the behavior of this morphogen gradient, making it necessary to develop a model that can recapitulate the biological features of the system, including its dynamic and scaling properties.

We present a biologically realistic 2-D model of the dynamics of the Bcd gradient in *Drosophila* embryos. This model is based on equilibrium binding of Bcd molecules to non-specific, low affinity DNA sites throughout the *Drosophila* genome. It considers both the diffusion media within which the Bcd gradient is formed and the dynamic and other relevant properties of *bcd* mRNA from which Bcd protein is produced. Our model recapitulates key features of the Bcd protein gradient observed experimentally, including its scaling properties and the stability of its nuclear concentrations during development. Our simulation model also allows us to evaluate the effects of other biological activities on Bcd gradient formation, including the dynamic redistribution of *bcd* mRNA in early embryos. Our simulation results suggest that, in our model, Bcd protein diffusion is important for the formation of an exponential gradient in embryos.

The 2-D model described in this report is a simple and versatile simulation procedure, providing a quantitative evaluation of the Bcd gradient system. Our results suggest an important role of Bcd binding to non-specific, low-affinity DNA sites in proper formation of the Bcd gradient in our model. They demonstrate that highly complex biological systems can be effectively modeled with relatively few parameters.

Development is a robust process that must achieve a precise and reproducible outcome despite inevitable variations among individual embryos [1], [2], [3], [4]. One such variation is the size of embryos. In *Drosophila*, patterning along the anterior-posterior (A–P) axis exhibits scaling properties such that embryonic structures are patterned in proportion to embryo length [5], [6], [7], [8]. Bicoid (Bcd), a morphogenetic protein in *Drosophila*, is responsible for instructing patterning along the A–P axis [9], [10], [11]. How developmental scaling is achieved is a subject of intense interest [12], [13], [14]. Several models have been proposed to explain scaling along the A–P axis in *Drosophila* but they currently remain at a theoretical level [15], [16], [17], [18], [19], [20], [21], [22]. For example, it was proposed recently that, theoretically, scaling could be achieved through a cytoplasmic streaming field [21], but the proposed streaming has not yet been observed experimentally.

Recent experimental studies have significantly advanced our knowledge on Bcd gradient behavior in *Drosophila* embryos. These new findings represent an important foundation for developing and evaluating models of Bcd gradient formation. In particular, a recent live-imaging study demonstrated that Bcd concentrations inside individual nuclei exhibit a striking stability even though the number of nuclei is doubling after each nuclear division in early embryos [23], [24]. These findings are important because they suggest that the positions of specific nuclear Bcd concentration thresholds are maintained even though the embryo itself, during this critical period of development, is undergoing an exponential growth in terms of the number of nuclei. In another study using stained embryos, we observed that the Bcd gradient itself exhibits properties of scaling [7]. In particular, the amount of Bcd in the anterior is correlated with embryo size. Furthermore, Bcd intensity is more precise when measured as a function of normalized A–P position than without such normalization. In embryos from *staufen* (*stau*) females, a loss of Bcd gradient scaling can directly explain the scaling defect of its target expression boundary. These findings suggest that the scaling properties of the Bcd gradient itself are important for scaled patterning along the A–P length in early embryos.

According to a widely-held view [25], Bcd gradient formation can be described as a dynamical process of localized protein production and embryo-wide diffusion and degradation. To fully describe and understand the Bcd gradient system, it is necessary to have a realistic model that considers the geometry of the embryo, the changing diffusion media within which the Bcd gradient is formed, and key features of the *bcd* mRNA from which Bcd is produced. For example, unlike what is assumed in most of the current models, the maternally deposited *bcd* mRNA is not restricted to a single point at the anterior tip [26], [27], [28], [29]. Instead, it forms a cloud-like shape in early embryos with its own volume, shape, density and location, properties that affect where Bcd molecules are produced in the embryo and, thus, how the gradient is formed. How the distribution of *bcd* mRNA may affect the shape and formation of the Bcd gradient remains an important question [30]. This question is further highlighted by a recent report that challenged the diffusion model [31]. Based on the observed dynamic distributions of *bcd* mRNA in early embryos, it was suggested that the shape and formation of the Bcd protein gradient are dictated by *bcd* mRNA distributions.

In this report, we develop a biologically realistic 2-D model to simulate the formation of the Bcd gradient in *Drosophila* embryos. It allows us to evaluate, among other things, the effect of *bcd* mRNA redistribution on Bcd gradient formation, suggesting that, in our model, Bcd protein diffusion is important for the formation of an exponential protein gradient. We also present simulation results that recapitulate several other important features of the Bcd gradient, including the experimentally observed scaling properties of the Bcd gradient and the stability of nuclear Bcd concentrations between different nuclear cycles.

A primary goal of our current study is to establish a biologically realistic model that simulates the dynamics of Bcd protein distribution in early *Drosophila* embryos, focusing on evaluating the effects of the changing diffusion media within which the Bcd gradient is formed and the molecular properties of *bcd* mRNA from which Bcd protein is produced. Since stochastic properties of the Bcd gradient system have been extensive investigated recently [32], [33], [34], they are not the focus of our current work. Here, the fluctuation in the number of Bcd molecules is averaged within a finite volume in our simulations [35], [36], [37], [38] to reveal the “average” properties of the Bcd gradient. In our model, we consider two forms of Bcd molecules: free-diffusing molecules, B_{free}, that are in the cytoplasm, and immobile molecules, B_{bound}, that are bound to low-affinity DNA sites inside the nuclei. The binding/unbinding behavior of Bcd to these DNA sites is characterized by the rate constants *k _{1}* and

(1)

At equilibrium, the ratio of bound to free Bcd molecules is determined by the association constant *K _{A}*=

(2)

(3)

Since only the unbound Bcd molecules in the cytoplasm can diffuse freely in our model, the effective diffusion constant *D _{eff}* can be expressed as:

(4)

where *D* is the diffusion constant for free-diffusing Bcd molecules inside the cytoplasm.

In our 2-D model, we simulate the mid-coronal slice of the *Drosophila* embryo, where Bcd molecules are distributed. An ellipse is selected to represent the boundary of the embryo slice, with major axis *a _{major}* and minor axis

(5)

where [DNA_{sites}]_{10} is the concentration of low-affinity DNA sites for Bcd within the cortical layer at nuclear cycle 10.

In our simulations, we assign a status to each cube in our mesh grid system depending on its center coordinate and developmental time: the cortical layer (I) and the inner part (II); prior to the formation of the cortical layer, all cubes within the embryo boundary are treated as status II as discussed above. We further assign status III to cubes that contain *bcd* mRNA and, thus, produce Bcd protein at a constant rate *J* (mol/s). In the main model simulations (Figs. 1–5),5), the *bcd* mRNA sphere is modeled as a disc (orthographic projection) inside the simulated embryo slice, with a uniform density, a radius of 45 µm and a center coordinate of (75 µm, 0); the shape, location and amount of *bcd* mRNA remain constant at all nuclear cycles in these simulations. The *J* value is identical for all cubes that contain *bcd* mRNA; thus the effective Bcd production rate for the entire embryo is the aggregate *J* (i.e., individual cube *J* value times cube number). At time *t*, Bcd molecules within each cube are in equilibrium with the low-affinity DNA sites, and we use the Euler forward method () to numerically calculate [B_{tot}] in cube (*i*,*j*) at time *t*+Δ*t* as a result of protein diffusion to and from its four immediate neighbors and position-independent protein degradation:

(6)

where *υ _{i,j}* is the volume of cube (

(7)

and

(8)

In our simulations we used Δ*t*=0.2s; we also tested Δ*t*=0.01, 0.1, 0.5 and 1.0s and obtained consistent results based on the values for the three established criteria to quantify simulated Bcd gradient behavior (see below), indicating that our solutions are numerically stable.

The developmental time in all our simulations is based on experimental observations [24], [39]. In our simulations, *t*=0 represents the time of egg deposition and our simulations end at 15 minutes into nuclear cycle 14. Nuclear cycles 10–14 start at 81, 95, 110, 125 and 140 minutes after egg deposition, respectively. The Bcd concentration profiles for nuclear cycles 10–14 shown in this report represent simulated data at 94, 109, 124, 139 and 154 minutes after egg deposition unless otherwise stated. The parameter values used in our simulations are either calculated based on estimates (for [DNA_{sites}]_{10}, see below) or chosen through a systematical sampling within the biologically reasonable ranges to yield Bcd gradient properties that satisfy three established criteria (see below). In our main model (Figs. 1–5),5), parameter values are as follows unless stated otherwise: *D*=2 µm^{2}s^{−1}, *ω*=0.00005 s^{−1}, *K _{A}*=5×10

We establish three different criteria to quantify the behavior of the simulated Bcd gradient. First, we measure the stability of nuclear Bcd concentration [B_{n}] during nuclear cycles 10–14 by computing the space-dependent relative change of [B_{n}] between the two successive nuclear cycles and then average the differences over space and nuclear cycles: . We consider a [B_{n}] gradient profile to be stable if |*g*|<0.1 [24]. Second, we measure the enrichment of Bcd molecules within the cortical layer using the ratio of total Bcd molecules within the cortical layer to those within the inner part of the embryo. We regard Bcd molecules to be enriched in the cortical layer if *Ratio*>1.5 at nuclear cycle 14. Third, we evaluate the shape of the Bcd gradient using its length constant *λ*, and we assume a gradient to be properly shaped if *λ/L* has a value between 0.15 and 0.2 [5], [7], [23]. To evaluate how well a simulated Bcd profile in a given region of the embryo fits an exponential function, we also calculate Adjusted *R*-square (Adjusted *R*^{2}) values [41]. Curve Fitting Tool in Matlab was used for fitting the exponential function .

Fig. 1A represents a simulated 2-D embryo at nuclear cycle 14 showing the local total Bcd concentration [B_{tot}]. These simulation results realistically recapitulate the biological system and reveal several important features that are consistent with experimental observations. First, consistent with experimental observations [5], [7], [9], [24], [42], Bcd is concentrated within the cortical layer (Fig. 1A). Within this layer, the local DNA-bound Bcd concentration forms an exponential gradient (see Fig. 1B,C) with a length constant λ of 92 µm, in agreement with experimentally measured values [5], [7]. Second, unlike other models in which Bcd production is restricted to a single point at the anterior resulting in a Bcd concentration gradient that follows an exponential function throughout the entire A–P length, our simulated results realistically reproduce the experimentally observed “deviation” at the anterior part of the embryo, a region where Bcd concentrations are known not to follow the exponential function ([5], [7], [42]; see Fig. S1 for a comparison between simulated data and experimental data). While the simulated data in the region of *x*/*L*=0.2 to 0.7 exhibit an excellent exponential fit (Adjusted *R*^{2}=0.9998; same value also obtained for the fitting region of *x*/*L*=0.2 to 1), the exponential fitting becomes significantly worse when the simulated data points from the most anterior part of the embryo are included (Adjusted *R*^{2}=0.9296 and 0.9461 for the fitting regions of *x*/*L*=0 to 0.7 and *x*/*L*=0 to 1, respectively). The observed anterior “deviation” of the Bcd gradient profiles in our 2-D model is due to the size and location of *bcd* mRNA (see Fig. S1 for 2-D simulation results where *bcd* mRNA is strictly located at the anterior tip).

One of the most striking features of the Bcd gradient dynamics is the stability of nuclear Bcd concentrations during nuclear cycles 10–14 [24]. To specifically evaluate our model on this critical feature, we simulated Bcd gradient formation and analyzed Bcd concentration profiles at distinct nuclear cycles. We plotted two separate profiles: local DNA-bound Bcd concentration [B_{bound}] and nuclear Bcd concentration [B_{n}], which are shown in Fig. 2A and B, respectively. In these figures, Bcd concentration profiles at nuclear cycles 10–14 are shown. Our results demonstrate that, even as the nuclear number doubles after each nuclear division during this period, the profiles of nuclear Bcd concentrations [B_{n}] remain stable (Fig. 2B). Consistent with the live-imaging data [24], our simulated [B_{n}] profiles maintains <10% variation between nuclear cycles 10–14 throughout the entire A–P length. Similar to nuclear cycle 14 (Fig. 1B,C), the simulated [B_{n}] profiles at other nuclear cycles also follow an exponential function with a length constant *λ* consistent with experimental values and exhibit the proper anterior “deviation” from an exponential function (see Fig. 2 legend for details).

As suggested by our recent studies [7], the scaling properties of embryonic patterning along the A–P axis can be directly traced to the scaling properties of the Bcd gradient itself. To determine whether our model can recapitulate scaling properties of the Bcd gradient, we conducted two different analyses. In our simulations, we assume that the amount of *bcd* mRNA is correlated with the embryo volume. This is a reasonable assumption (but remains to be tested experimentally) if the amounts of maternally-deposited components, including *bcd* mRNA, are proportional to the egg volume, i.e., the concentrations of various maternal components remain as constants among individual eggs even though the total maternally-deposited contents (i.e., egg volume) can vary between eggs. In our first analysis, we simulated two individual embryos that have distinct lengths (550 µm and 600 µm). Fig. 3A and B show the local DNA-bound Bcd concentration profiles at nuclear cycle 14 expressed as a function of, respectively, absolute distance from the anterior *x* (in µm) and fractional embryo length *x*/*L*. Our results show that, while these two profiles in the anterior and middle portions of the embryo are apart from each other in the *x* plot, they converge in the *x*/*L* plot, most notably in the broad, mid-portion of the embryo. When the same *bcd* mRNA amount was applied to these two embryos, such a convergence was not observed (data not shown). In a second analysis, we simulated a group of 50 embryos that differ in size (see Fig. 3 legend). We calculated [B_{bound}] noise (standard deviation divided by the mean) at nuclear cycle 14 for these embryos as a function of either *x* or *x*/*L*. Our results (Fig. 3C) show a significantly lower [B_{bound}] noise in *x*/*L* than in *x*, a finding that is fully consistent with experimental data ([7]; note that all of the simulated [B_{bound}] noise levels are lower than experimental values because parameters, except *L*, do not have variations in our simulations). In contrast, when there is no correlation between the *bcd* mRNA amount and embryo size in our simulations, [B_{bound}] noise actually became higher as a function of *x*/*L* than of *x* (Fig. 3D). These two sets of analyses demonstrate that, with a simple assumption that the amount of *bcd* mRNA deposited into an egg is proportional to egg volume, our model can readily reproduce the experimentally observed scaling properties of the Bcd gradient.

Although the precise values for the parameters *D* and *ω* remain either controversial or unknown at this time [24], [43], [44], the values chosen in our main model simulations described thus far were selected from a systematic sampling to yield Bcd gradient properties consistent with experimental observations as evaluated by the three criteria established above. Here, we present our sampling results in Fig. 4, with individual panels showing the identified parameter space satisfying each of the three criteria. While the space satisfying the cortical enrichment criterion is relatively broad (Fig. 4C), the space that satisfies the other two criteria is narrower (Fig. 4A, B), particularly for the gradient shape criterion *λ* (Fig. 4B). To further illustrate the effects of individual parameters on model performance, we also plot the criterion values as a function of each of the three parameters while holding the other two parameters at set values used in the main model simulation (Fig. S2). Our results show that the operating ranges of both *D* and *ω* values in our main model are constrained to satisfy the criterion of *λ* (Fig. S2E, H), indicating that, as expected, these parameters are particularly important for determining the gradient shape.

Our parameter evaluation analysis also shows that the upper and lower limits of the operating range of *K _{A}* values are constrained to satisfy the [B

In our main model, the enrichment of Bcd molecules within the cortical layer and effective *D* are governed by the binding/unbinding chemical equilibrium of Bcd molecules to low-affinity DNA sites throughout the *Drosophila* genome. This model is consistent with previous experimental findings that the DNA-binding homeodomain itself plays a critical role in nuclear localization [47]. Our main model assumes that the disappearance of the nuclear structure at the mitotic phase does not disrupt the chemical equilibrium and, therefore, it does not consider the mitotic process separately. To specifically evaluate the effects of mitosis, we conducted simulations by incorporating a 3-minute mitotic process during which all Bcd molecules are allowed to diffuse freely (as a result of chromosome condensation that prevents Bcd from binding to non-specific DNA sites) throughout the entire embryo, thus *D _{eff}* =

A recent study challenges the diffusion model of Bcd gradient formation, raising the question of whether this process is dependent on Bcd protein diffusion [31]. We used our 2-D simulation model to specifically evaluate the effects of *bcd* mRNA dynamics on Bcd gradient formation. Our simulations are based on the available experimental data [31] and assume that *bcd* mRNA is distributed uniformly within a sphere located in the anterior during nuclear cycles 1–9 as in our main model but as a gradient-like profile within the cortical layer during nuclear cycles 10–14 (see Fig. 6 legend for further details). Since there is currently insufficient knowledge about the mechanisms leading to *bcd* mRNA redistribution, our simulations simple assume that such redistribution takes place at the onset of nuclear cycle 10 without modeling specifically this process. Fig. 6A shows the local total Bcd concentration [B_{tot}] in a simulated embryo with redistributed *bcd* mRNA [31] and a diffusion constant value (*D*=2 µm^{2}s^{−1}) identical to that used in Fig. 1. The simulated profiles of nuclear Bcd concentrations [B_{n}] remain stable during nuclear cycles 10–14 (Fig. 6B) and exhibit an exponential function (Fig. 6D, blue line; also see below). The length constant *λ* of the nuclear Bcd profiles at nuclear cycle 14 in these simulations is 145 µm; while this value is somewhat higher than experimental measurements, those at earlier nuclear cycles are smaller and closer to experimental values (107, 124, 134 and 141 µm at nuclear cycles 10–13, respectively). These results show that, with modest parameter adjustments, our model can yield, from the redistributed *bcd* mRNA profiles, nuclear Bcd concentration profiles and dynamics properties that are broadly consistent with experimental observations.

To specifically determine whether Bcd protein diffusion might play a role in proper gradient formation, we reduced *D* progressively in our simulations. Our results show that as *D* is reduced, [B_{bound}] profiles at nuclear cycle 14 become progressively less exponential (Fig. 6C,D; profiles are shown only at selected *D* values for clarity). The Adjusted *R*^{2} values of the exponential fit (in the region of *x*/*L*=0.2 to 0.7) decreased from 0.9854 at *D*=2 µm^{2}s^{−1}, to 0.9708 and 0.8987 at *D*=1 and 0.1 µm^{2}s^{−1}, respectively (see Fig. 6 legend for further details). Fig. S4A shows the Adjusted *R*^{2} values as a function of *D* for both the main model and the model with *bcd* mRNA redistribution. While Adjusted *R*^{2} is insensitive to *D* in our main model, this value precipitously drops as *D* decreases in the mRNA redistribution simulation. These results show that, in our model, Bcd protein diffusion is important for the formation of an exponential Bcd protein gradient. Our simulation results reveal two effects of the cortical localization of *bcd* mRNA: 1) it enhances the enrichment of Bcd molecules within the cortical layer relative to the inner part of the embryo (Fig. S4B), and 2) it enlarges the range of *D* that satisfies the stability criterion *g* (Fig. S4C). These findings provide potential biological roles of the recently reported dynamic redistribution of *bcd* mRNA.

In our main model, the size of the nuclei is assumed to be a constant between different nuclear cycles, hence the relationship [B_{n}][B_{bound}]/nuclear number. Furthermore, our main model assumes that the thickness of the cortical layer remains as a constant between different nuclear cycles. To determine whether changes in nuclear volume and cortical layer may affect the [B_{n}] stability, we performed simulations where both the cortical layer thickness and nuclear size vary between nuclear cycles. In these simulations, [B_{n}] calculations were based on estimated nuclear diameter of 10, 10.5, 9.2, 8.2 and 6.5 µm for nuclear cycles 10–14, respectively [24]. Fig. S5B and S5C show, respectively, the simulated [B_{bound}] and [B_{n}] profiles at different nuclear cycles (also see Fig. S5A for a 3-D presentation of [B_{tot}] in the simulated embryo). These results show that our model, with proper parameter adjustments, can produce stable [B_{n}] profiles when incorporating changes in both nuclear size and cortical layer thickness between different nuclear cycles (see Fig. S5 legend for further details).

The 2-D simulation model described in this report captures key properties of the Bcd gradient system in a biologically realistic manner. This model considers changes in both nuclear number and their relative locations in early embryos. Consistent with experimental findings [24], our simulated results show that the profiles of nuclear Bcd concentrations [B_{n}] are stable during nuclear cycles 10–14 (Fig. 2B). Our 2-D model also reproduces the scaling properties of the Bcd gradient in *Drosophila* embryos with a simple assumption that the amount of *bcd* mRNA (thus/or the rate of Bcd production) is proportional to the embryo volume. Our simulation results (Fig. 3) are consistent with the experimentally-observed scaling properties the Bcd gradient [7]. Finally, our model allows us to specifically evaluate other biological features of the system, including the mitotic process and *bcd* mRNA redistribution. Our results show that formation of an exponential Bcd gradient is dependent on Bcd protein diffusion in our model (Fig. 6). They also offer potential biological roles of the recently observed dynamic redistribution of *bcd* mRNA (Fig. S4).

A recent model [48] based on “nuclear trapping” of Bcd molecules has also successfully recapitulated the stability of the nuclear Bcd gradient during development. While both the nuclear trapping model and our current model support an important role of nuclei in Bcd gradient formation [24], [44], they differ from each other in the following fundamental aspects. First, the Coppey et al. model assumes that there is no Bcd-protein degradation in early embryos. Although the mechanisms controlling Bcd stability remain elusive, the experimentally observed rapid disappearance of Bcd after cellularization [9] indicates that Bcd is subject to active protein degradation in early embryos. Our model, while specifically considering Bcd degradation, is able to recapitulate the experimentally demonstrated stability of [B_{n}] profiles during development (see Fig. S2G–I for the effects of *ω* on the three established criteria, with Fig. S2H showing the operating range of *ω* that satisfies the criterion of *λ*). Second, while the Coppey et al. model assumes a nuclear barrier-mediated equilibrium between Bcd molecules inside and outside the nucleus, our main model assumes a chemical equilibrium between free Bcd molecules and those that are bound to non-specific DNA sites inside the nucleus. Since nuclear localization of a protein is insufficient for proper gradient formation in embryos [44], it is evident that additional properties of Bcd, such as its DNA binding function as modeled in our study, are likely to play an important role in this process. Finally, unlike the Coppey et al. model, our model does not restrict Bcd protein diffusion in the inner part (yolk) of the embryo during nuclear cycles 10–14. It has been observed experimentally that the yolk nuclei, similar to those in the cortical layer, are also enriched with Bcd molecules [24]. We conducted simulations in which a fraction of the nuclei at the onset of nuclear cycle 10 remain in the inner part of the embryo, and we observed an enrichment of Bcd molecules in these “yolk nuclei” as seen experimentally (data not shown).

The 2-D model described in this report represents a simple and versatile tool to simulate the Bcd gradient system in a biologically realistic way. In contrast to a recent model that is also sufficient to explain [B_{n}] stability [33], our simulation model takes a distinct approach and requires significantly fewer physical quantities that should be experimentally measurable and verifiable. Effectively, the biologically realistic properties of the Bcd gradient achieved in our simulations are dependent only on three relevant parameters: the diffusion constant of free Bcd molecules in the cytoplasm, the degradation rate of Bcd molecules, and the fraction of bound to free Bcd molecules at nuclear cycle 10 (see above; also see Eq. 2). While the fact that a complex biological system can be modeled with only three effective parameters illustrates the effectiveness and simplicity of our model, the choice of a limited number of parameters also poses inevitable limitations. For example, our current model does not consider the precise mechanisms of nuclear-cytoplasmic shuttling of the Bcd molecules, which is likely to be affected by the structures of nuclear pores and/or cytoplasmic islands [24], [49], [50], [51] and may be important for finer aspects of the system's behavior. We emphasize that, since our current model can recapitulate the key biological properties of the system, these three effective parameters likely represent the principal physical quantities most relevant to the system behavior.

In addition to the biological features evaluated in the current study, the versatility of our simulation model should allow us in the future to evaluate the effects of other biological features and physical or genetic perturbations on Bcd gradient formation, such as *bcd* mRNA localization [29], [52], [53], Bcd protein diffusion or degradation, or developmental clock distortions [5], [54], [55]. As an example, we have recently identified distinct features of the Bcd gradient profiles on the dorsal and ventral sides of the embryo (unpublished results). Our 2-D simulation model based on simple geometric asymmetry of the embryo can fully recapitulate our experimentally observed differences of the Bcd gradient profiles (unpublished results), differences that cannot be adequately modeled by 1-D simulations. Our simple but biologically realistic 2-D model represents a useful platform that should facilitate the formulation of experimentally testable hypotheses. It will expand the toolbox available to our studies of the fundamental mechanisms of developmental processes.

See the Model section in Results and Discussion for details about methods.

Comparison between simulated data and experimental data. A. Shown are simulated [B_{n}] at nuclear cycle 14 (red and green) and experimentally observed Bcd gradient profile also at early nuclear cycle 14 (blue). The experimental data shown here are from He et al. [7], which represent background-subtracted Bcd intensities, with error bars (standard deviation) shown. The two simulated [B_{n}] profiles are obtained from simulations identical to the main model simulations expect the center coordinate of bcd mRNA was fine tuned to yield a [B_{n}] profile that matches the experimentally observed Bcd profile (50 µm, 0; red) or the bcd mRNA is restricted to a single cube at the anterior tip of the embryo (0, 0; green). The Adjusted R2 values for our experimental data within the fitting ranges of x/L=0.2 to 0.7 and 0 to 0.7 are 0.9961 and 0.9851, respectively. We note that both simulated Bcd profiles have lower levels in the posterior part of the embryo than the experimentally observed profile, a difference whose biological relevance will require further experimental and modeling investigations. B. Same as in A, except on ln scale. While the simulated red profile matches well with the experimental data and exhibits the experimentally observed anterior “deviation,” the simulated green profile clearly fails to exhibit this property.

(0.51 MB TIF)

Click here for additional data file.^{(494K, tif)}

Evaluating the effects of parameter values in different simulations. Parameter values for KA (panels A–C), D (panels D–F) and ω (panels G–I) are systematically altered to evaluate model performance. Three criteria are used in these evaluations: [B_{n}] stability as measured by g (panels A, D and G), gradient shape as measured by length constant λ at nuclear cycle 14 (panels B, E, H), and cortical enrichment as measured by the ratio of total Bcd molecules in the cortical layer to those in the inner part of the embryo at nuclear cycle 14 (panels C, F and I). The regions where the Bcd gradient profiles satisfy the established criteria, i.e., |g|<0.1, Ratio>1.5, and 1.5<λ/L<2.0, are shaded. In these analyses, individual parameters are systematically changed in simulations when the other two parameters are at set values at their respective model simulations. Different colors represent different simulation procedures as indicated in the figure.

(0.76 MB TIF)

Click here for additional data file.^{(742K, tif)}

2-D simulation with mitosis. A. A simulated embryo at nuclear cycle 14 showing [B_{tot}] (arbitrary units). In this simulation, the mitotic process is specifically considered at nuclear cycles 10–13, during which all Bcd molecules are allowed to diffuse freely. Other parameter values are identical to those used in the main model. The A–P position is shown as absolute distance x (in µm) from the anterior. At nuclear cycle 14, the ratio of total Bcd molecules in the cortical layer to those in the inner part of the embryo is 1.7975. B. A plot of [B_{bound}] (arbitrary units) within the cortical layer as a function of x/L, at nuclear cycles 10–14. C. Same as in B, except now showing [B_{n}] within the cortical layer at nuclear cycles 10–14. Similar to the main model (Fig. 2B), the mitotic process does not affect [B_{n}] stability (g=0.0143). The simulated [B_{n}] profile at nuclear cycle 14 has a length constant λ=93.3 µm. See Fig. S2 for a model performance comparison.

(0.73 MB TIF)

Click here for additional data file.^{(711K, tif)}

Investigating the effects of bcd mRNA redistribution. A. Adjusted R2 value of the exponential fitting (within the range of x/L=0.2 to 0.7) of simulated [B_{bound}] is plotted as a function of D. For comparative purposes, results obtained from both the main model simulation (red) and the simulation with bcd mRNA redistribution (blue) are shown. Note the difference in Adjusted R2 sensitivity to D. B. The ratio of total Bcd molecules in the cortical layer to those in the inner part of the embryo plotted as a function of D. Color codes are the same as in A. Regions where Ratio >1.5 are shaded. Note the higher Ratio value obtained in the mRNA redistribution simulation (blue). C. A plot of g as a function of D, with color codes being the same as in A. Regions where |g|<0.1 are shaded. Note the blue curve is within the shaded area under all D values tested, suggesting another potential role of the observed bcd mRNA redistribution.

(0.54 MB TIF)

Click here for additional data file.^{(524K, tif)}

Evaluating the effects of nuclear size and cortical layer at different nuclear cycles. A. A simulated embryo at 5 min into nuclear cycle 14 showing [B_{tot}] (arbitrary units). In this simulation, non-specific DNA binding site concentrations were calculated based on the following estimates: the thickness of the cortical layer is 15, 20, 20 23, and 25 µm for nuclear cycles 10–14, respectively, corresponding to a cortical layer volume of ~3.3, 4.2, 4.2 4.7, and 5.0 nl. The concentrations of the non-specific DNA binding sites within the cortical layer are: 5×10^{−7}, 4×10^{−7}, 4×10^{−7}, 3.6×10^{−7}, and 3.4×10^{−7} M, for nuclear cycles 10–14, respectively. The relative volumes of a single nucleus, which were calculated based on the experimental estimates [24], are 3.64, 4.21, 2.83, 2.0, and 1.0 for nuclear cycles 10–14, respectively. Parameter values used in this simulation are: D=6 µm^{2}s^{−1}, ω=0.0004 s^{−1}, KA=2.4×10^{6} M^{−1}. At nuclear cycle 14, the ratio of total Bcd molecules in the cortical layer to those in the inner part of the embryo is 5.0044. B. A plot of [B_{bound}] (arbitrary units) within the cortical layer as a function of x/L, at nuclear cycles 10–14. C. Same as in B, except now showing nuclear Bcd concentrations [B_{n}] at nuclear cycles 10–14. [B_{n}] is calculated from [B_{bound}] at each nuclear cycle based on the nuclear number and the volumes of each nucleus and the cortical layer. As seen in the main model (Fig. 2B), [B_{n}] profiles exhibit stability between different nuclear cycles (g=−0.050) when changes in nuclear volume and cortical layer between different nuclear cycles are incorporated in our simulation. The [B_{n}] profile at nuclear cycle 14 has a length constant λ=105 µm. We note that the results shown in Fig. S5 do not represent an improvement over those obtained in the main model simulations. It is possible that the use of additional parameters may improve the simulation outcome when more biological features are included in our model.

(0.75 MB TIF)

Click here for additional data file.^{(728K, tif)}

We thank members of our groups for discussions and anonymous reviewers for constructive suggestions.

**Competing Interests: **The authors have declared that no competing interests exist.

**Funding: **This work is supported in part by grants from the National Science Foundation (IOS-0843424) and the National Institutes of Health (GM072812, GM78381). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

1. Martinez Arias A, Hayward P. Filtering transcriptional noise during development: concepts and mechanisms. Nat Rev Genet. 2006;7:34–44. [PubMed]

2. Kerszberg M, Wolpert L. Specifying positional information in the embryo: looking beyond morphogens. Cell. 2007;130:205–209. [PubMed]

3. Lander AD. Morpheus unbound: reimagining the morphogen gradient. Cell. 2007;128:245–256. [PubMed]

4. Lewis J. From signals to patterns: space, time, and mathematics in developmental biology. Science. 2008;322:399–403. [PubMed]

5. Houchmandzadeh B, Wieschaus E, Leibler S. Establishment of developmental precision and proportions in the early Drosophila embryo. Nature. 2002;415:798–802. [PubMed]

6. Lott SE, Kreitman M, Palsson A, Alekseeva E, Ludwig MZ. Canalization of segmentation and its evolution in Drosophila. Proc Natl Acad Sci U S A. 2007;104:10926–10931. [PubMed]

7. He F, Wen Y, Deng J, Lin X, Lu J, et al. Probing intrinsic properties of a robust morphogen gradient in Drosophila. Dev Cell. 2008;15:558–567. PMID: 18854140. [PMC free article] [PubMed]

8. Manu, Surkova S, Spirov AV, Gursky VV, Janssens H, et al. Canalization of gene expression in the Drosophila blastoderm by gap gene cross regulation. PLoS Biol. 2009;7:e1000049. [PMC free article] [PubMed]

9. Driever W, Nüsslein-Volhard C. A gradient of bicoid protein in Drosophila embryos. Cell. 1988;54:83–93. [PubMed]

10. Struhl G, Struhl K, Macdonald P. The gradient morphogen bicoid is a concentration-dependent transcriptional activator. Cell. 1989;57:1259–1273. [PubMed]

11. Ephrussi A, St. Johnston D. Seeing is believing. The bicoid morphogen gradient matures. Cell. 2004;116:143–152. [PubMed]

12. Day SJ, Lawrence PA. Measuring dimensions: the regulation of size and shape. Development. 2000;127:2977–2987. [PubMed]

13. Patel NH, Lall S. Precision patterning. Nature. 2002;415:748–749. [PubMed]

14. Ben-Zvi D, Shilo BZ, Fainsod A, Barkai N. Scaling of the BMP activation gradient in Xenopus embryos. Nature. 2008;453:1205–1211. [PubMed]

15. Howard M, Ten Wolde PR. Finding the center reliably: robust patterns of developmental gene expression. Physical Rev Lett. 2005;95:208103. [PubMed]

16. Houchmandzadeh B, Wieschaus E, Leibler S. Precision domain specification in the developing Drosophila embryo. Pysical Rev E. 2005;72:061920. [PubMed]

17. Aegerter-Wilmsen T, Aegerter CM, Bisseling T. Model for the robust establishment of precise proportions in the early Drosophila embryo. J Theor Biol. 2005;234:13–19. [PubMed]

18. McHale P, Rappel WJ, Levine H. Embryonic pattern scaling achieved by oppositely directed morphogen gradients. Phys Biol. 2006;3:107–120. [PubMed]

19. Ishihara S, Kaneko K. Turing pattern with proportion preservation. J Theor Biol. 2006;238:683–693. [PubMed]

20. Bergmann S, Sandler O, Sberro H, Shnider S, Schejter E, et al. Pre-steady-state decoding of the Bicoid morphogen gradient. PLoS Biol. 2007;5:e46. [PMC free article] [PubMed]

21. Hecht I, Rappel WJ, Levine H. Determining the scale of the Bicoid morphogen gradient. Proc Natl Acad Sci U S A. 2009;106:1710–1715. [PubMed]

22. Morishita Y, Iwasa Y. Accuracy of positional information provided by multiple morphogen gradients with correlated noise. Phys Rev E Stat Nonlin Soft Matter Phys. 2009;79:061905. [PubMed]

23. Gregor T, Tank DW, Wieschaus EF, Bialek W. Probing the limits to positional information. Cell. 2007;130:153–164. [PMC free article] [PubMed]

24. Gregor T, Wieschaus EF, McGregor AP, Bialek W, Tank DW. Stability and nuclear dynamics of the bicoid morphogen gradient. Cell. 2007;130:141–152. [PMC free article] [PubMed]

25. Wolpert L. Positional information and the spatial pattern of cellular differentiation. J Theor Biol. 1969;25:1–47. [PubMed]

26. Berleth T, Burri M, Thoma G, Bopp D, Richstein S, et al. The role of localization of bicoid RNA in organizing the anterior pattern of the Drosophila embryo. EMBO J. 1988;7:1749–1756. [PubMed]

27. Frigerio G, Burri M, Bopp D, Baumgartner S, Noll M. Structure of the segmentation gene paired and the Drosophila PRD gene set as part of a gene network. Cell. 1986;47:735–746. [PubMed]

28. St. Johnston D, Driever W, Berleth T, Richstein S, Nüsslein-Volhard C. Multiple steps in the localization of bicoid mRNA to the anterior pole of the Drosophila oocyte. Development. 1989;Suppl.:13–19. [PubMed]

29. Crauk O, Dostatni N. Bicoid determines sharp and precise target gene expression in the Drosophila embryo. Curr Biol. 2005;15:1888–1898. [PubMed]

30. Lipshitz HD. Follow the mRNA: a new model for Bicoid gradient formation. Nat Rev Mol Cell Biol. 2009;10:509–512. [PubMed]

31. Spirov A, Fahmy K, Schneider M, Frei E, Noll M, et al. Formation of the bicoid morphogen gradient: an mRNA gradient dictates the protein gradient. Development. 2009;136:605–614. [PubMed]

32. Wu YF, Myasnikova E, Reinitz J. Master equation simulation analysis of immunostained Bicoid morphogen gradient. BMC Syst Biol. 2007;1:52. [PMC free article] [PubMed]

33. Okabe-Oho Y, Murakami H, Oho S, Sasai M. Stable, precise, and reproducible patterning of bicoid and hunchback molecules in the early Drosophila embryo. PLoS Comput Biol. 2009;5:e1000486. [PMC free article] [PubMed]

34. Tostevin F, ten Wolde PR, Howard M. Fundamental limits to position determination by concentration gradients. PLoS Comput Biol. 2007;3:e78. [PMC free article] [PubMed]

35. Miura S, Shimokawa T, Nomura T. Stochastic simulations on a model of circadian rhythm generation. Biosystems. 2008;93:133–140. [PubMed]

36. Wilkie J, Wong Y. Positivity preserving chemical Langevin equations. Chem Phys. 2008;353:132–138.

37. Breuer HP, Honerkamp J, Petruccione F. A stochasitic approach to complex chemical reactions. Chem Phys Lett. 1992;190:199–201.

38. Gillespie DT. The chemical Langevin equation. J Chem Phys. 2000;113:297–306.

39. Foe VE, Alberts BM. Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesis. J Cell Sci. 1983;61:31–70. [PubMed]

40. Adams MD, Celniker SE, Holt RA, Evans CA, Gocayne JD, et al. The genome sequence of Drosophila melanogaster. Science. 2000;287:2185–2195. [PubMed]

41. Lattin J, Carroll JD, Green PE. Analyzing multivariate data. Pacific Grove, CA: Thompson Books/Cole; 2003.

42. Holloway DM, Harrison LG, Kosman D, Vanario-Alonso CE, Spirov AV. Analysis of pattern precision shows that Drosophila segmentation develops substantial independence from gradients of maternal gene products. Dev Dyn. 2006;235:2949–2960. [PMC free article] [PubMed]

43. Gregor T, Bialek W, van Steveninck RR, Tank DW, Wieschaus EF. Diffusion and scaling during early embryonic pattern formation. Proc Natl Acad Sci U S A. 2005;102:18403–18407. [PubMed]

44. Gregor T, McGregor AP, Wieschaus EF. Shape and function of the Bicoid morphogen gradient in dipteran species with different sized embryos. Dev Biol. 2008;316:350–358. [PMC free article] [PubMed]

45. Affolter M, Percival-Smith A, Muller M, Leupin W, Gehring WJ. DNA binding properties of the purified Antennapedia homeodomain. Proc Natl Acad Sci U S A. 1990;87:4093–4097. [PubMed]

46. Damante G, Tell G, Formisano S, Fabbro D, Pellizzari L, et al. Effect of salt concentration on TTF-1 HD binding to specific and non-specific DNA sequences. Biochem Biophys Res Commun. 1993;197:632–638. [PubMed]

47. Ghaffari M, Zeng X, Whitsett JA, Yan C. Nuclear localization domain of thyroid transcription factor-1 in respiratory epithelial cells. Biochem J. 1997;328:757–761. [PubMed]

48. Coppey M, Berezhkovskii AM, Kim Y, Boettiger AN, Shvartsman SY. Modeling the bicoid gradient: diffusion and reversible nuclear trapping of a stable protein. Dev Biol. 2007;312:623–630. [PMC free article] [PubMed]

49. Quimby BB, Corbett AH. Nuclear transport mechanisms. Cell Mol Life Sci. 2001;58:1766–1773. [PubMed]

50. DeLotto R, DeLotto Y, Steward R, Lippincott-Schwartz J. Nucleocytoplasmic shuttling mediates the dynamic maintenance of nuclear Dorsal levels during Drosophila embryogenesis. Development. 2007;134:4233–4241. [PubMed]

51. Mavrakis M, Rikhy R, Lippincott-Schwartz J. Cells within a cell: Insights into cellular architecture and polarization from the organization of the early fly embryo. Commun Integr Biol. 2009;2:313–314. [PMC free article] [PubMed]

52. Driever W, Nüsslein-Volhard C. The bicoid protein determines position in the Drosophila embryo in a concentration dependent manner. Cell. 1988;54:95–104. [PubMed]

53. Ferrandon D, Elphick L, Nusslein-Volhard C, St Johnston D. Staufen protein associates with the 3′UTR of bicoid mRNA to form particles that move in a microtubule-dependent manner. Cell. 1994;79:1221–1232. [PubMed]

54. Lucchetta EM, Lee JH, Fu LA, Patel NH, Ismagilov RF. Dynamics of Drosophila embryonic patterning network perturbed in space and time using microfluidics. Nature. 2005;434:1134–1138. [PMC free article] [PubMed]

55. Lucchetta EM, Vincent ME, Ismagilov RF. A precise Bicoid gradient is nonessential during cycles 11–13 for precise patterning in the Drosophila blastoderm. PLoS One. 2008;3:e3651. [PMC free article] [PubMed]

Articles from PLoS ONE are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |