## Abstract

Morphogen gradients infer cell fate as a function of cellular position. Experiments in *Drosophila* embryos have shown that the Bicoid (Bcd) gradient is precise and exhibits some degree of scaling. We present experimental results on the precision of Bcd target genes for embryos with a single, double or quadruple dose of *bicoid* demonstrating that precision is highest at mid‐embryo and position dependent, rather than gene dependent. This confirms that the major contribution to precision is achieved already at the Bcd gradient formation. Modeling this dynamic process, we investigate precision for inter‐embryo fluctuations in different parameters affecting gradient formation. Within our modeling framework, the observed precision can only be achieved by a transient Bcd profile. Studying different extensions of our modeling framework reveals that scaling is generally position dependent and decreases toward the posterior pole. Our measurements confirm this trend, indicating almost perfect scaling except for anterior most expression domains, which overcompensate fluctuations in embryo length.

## Introduction

Bicoid (Bcd) is a well studied morphogen involved in the patterning of the anterior–posterior (AP) axis of the *Drosophila* embryo (Driever and Nüsslein‐Volhard, 1988a, 1988b). Zygotic downstream genes read out this gradient and their expression domains determine the basic body plan of the embryo along this axis. The positions of these domains are remarkably insensitive to fluctuations in the external environment (Houchmandzadeh *et al*, 2002; Crauk and Dostatni, 2005; Lucchetta *et al*, 2005) and their relative proportions are maintained across embryos of different sizes. The latter feature is referred to as *scaling* and occurs within a single species (Houchmandzadeh *et al*, 2002; Lott *et al*, 2007) and also across different species (Gregor *et al*, 2005, 2008; Lott *et al*, 2007).

Recent experiments also show that the Bcd gradient itself is rather precise and that its length scale correlates to some extent to the embryo size (Gregor *et al*, 2007a; He *et al*, 2008). These new findings suggest that the precision and scaling of Bcd target genes may, at least in part, be attributed to that of the morphogen gradient itself. Therefore, we focus here on single morphogen models that aim at explaining the precision and scaling of Bcd target genes at the level of the morphogen gradient. However, it is likely that other mechanisms (e.g. Bcd interactions with the *staufen* gene (Aegerter‐Wilmsen *et al*, 2005), gap genes interactions (Jaeger *et al*, 2004a, 2004b, 2007; Jaeger and Reinitz, 2006; Manu *et al*, 2009a, 2009b), Bcd interactions with maternal Hunchback and the terminal system (Ochoa‐Espinosa *et al*, 2009) or bistability (Lopes *et al*, 2008)) also contribute to further increase robustness in AP patterning.

In this work, we assess the spatial precision of expression domains for the gap genes Krüppel (Kr), Giant (Gt) and Hunchback (Hb), as well as the pair‐rule gene Even‐skipped (Eve) in more than 150 staining images of embryos with a single, double and quadruple dose of *bcd*. Our data indicate that the precision is maximal at mid‐embryo as well as *position* dependent rather than *gene* dependent. This provides independent support for the conclusions drawn from direct measurements of Bcd (Gregor *et al*, 2007a; He *et al*, 2008) that the morphogen gradient itself is the main contributor to the precision of the target genes. It motivates our subsequent analytical investigation of noise propagation during Bcd gradient formation within a single morphogen modeling framework. Our analysis shows that fluctuations in the parameters affecting morphogen production, degradation, diffusion as well as nuclear trapping can give rise naturally to highest precision at mid‐embryo, provided the Bcd gradient is decoded at *pre‐steady*‐state. Finally, we investigate the scaling of the morphogen gradient in a variety of models that couple the embryo size to parameters affecting gradient formation. To single out the most likely scenario, we measure scaling at the level of the target genes. We find that expression domain boundaries posterior to ∼40% embryo length (*L*) exhibit almost perfect scaling, whereas in the more anterior region we observe over‐compensation to fluctuations in embryo size. This effect appears to be *position* dependent rather than *gene* dependent. Interestingly, it is in good agreement with a model that explicitly includes the nuclear trapping of Bcd. We conclude that the formation of the Bcd gradient itself is likely to be a main contributor to robust patterning along the AP axis and that pre‐steady‐state decoding and nuclear trapping are efficient means to increase robustness.

## Results

### The gap and pair‐rule gene expression domains are precise

We studied staining images of 154 *Drosophila melanogaster* embryos at cleavage cycle 14 (see Methods in Bergmann *et al*, 2007). Using image processing tools described in ‘Materials and methods’, we measured the relative positions *x*/*L* of the protein domains of Gt, Hb, Kr and Eve. We screened both wild‐type embryos with two copies of *bcd* and mutant strains with one or four *bcd* copies resulting in shifted expression domains. We observed that in mid‐embryo, the standard deviation of the relative domain localizations, σ(*x*/*L*), is between 1 and 2% and it tends to be higher (2–4%) toward the anterior and posterior poles (Figure 1).

These results are in good agreement with previous results on the Hb domain (Houchmandzadeh *et al*, 2002). Importantly, they show that the precision of the target genes has the same magnitude and positional trend as that of Bcd according to Gregor *et al* (2007a). Moreover, our experimental analysis shows that precision is more *position* dependent than *gene* dependent. Indeed, when expression domains are shifted because of different *bcd* mRNA dosages, the precision of the domain seems to change according to the new domain position. Overall, all the investigated domains follow similar precision trends. This suggests that a major contribution to the precision of the Bcd target genes can be attributed to the morphogen gradient itself. We note that position dependent precision may also arise as an experimental artifact when analyzing embryos with different orientations and varying ages (times classes T5–T8 of nuclear cycle 14) or because of imperfect scaling (cf. Supplementary Text S3). However, estimating the maximal contributions of such effects we still find that higher precision at mid‐embryo is statistically significant (cf. Supplementary Text S1; Supplementary Figure S4).

### Modeling precision in a single morphogen framework

Given the aforementioned pieces of evidence for achieving precision at the level of the gradient formation, we consider a French‐flag model (Wolpert, 1969) in which a single morphogen gradient induces expression domain boundaries in a concentration‐dependent manner (see Materials and methods). During the first stages of embryonic development, 13 nuclear divisions occur about every 10 min. Gregor *et al* (2007b) have recently shown that nuclear concentrations stabilize at cycle 10. From a robustness perspective, this means that Bcd read‐out is based on a similar number of molecules and so the associated stochastic noise does not change significantly from cycle 10 on. Therefore, we focus here on the noise propagated from the *external* Bcd gradient and investigate the fluctuations in domain localization when perturbing the various parameters characterizing the gradient formation (i.e. production, degradation, diffusion and nuclear trapping rates).

We want to assess the variation in domain position Δ*x*, which is induced by embryo‐to‐embryo fluctuations of magnitude Δ*q* in a parameter *q* affecting the morphogen gradient formation. Δ*x* can be estimated analytically to first order in Δ*q*

where we used that the threshold concentration is fixed, yielding d*M*=(∂*M*/∂*x*)d*x*+(∂*M*/∂*q*)d*q*=0. This estimate can be computed for any morphogen distribution *M*(*x*, *t*) with explicit dependence on *q*. Here, we consider the time‐dependent solution of a model in which Bcd is produced at the anterior pole with production rate *s*_{0}, diffuses according to a uniform diffusion constant *D*, is degraded at uniform rate α and is trapped and released by the nuclei in the embryo at rates *k*_{n}, *k*_{−n}, respectively (see Materials and methods). In Figure 2, we plotted the imprecision measure Δ*x* for small fluctuations (5%) in each of these parameters (Figure 2A–D) as well as a combination of fluctuations in the production and diffusion rates (Figure 2E) (numerical simulations gave very similar results, data not shown). We choose *t*_{i}=100 min as the time when patterning is initiated, but qualitatively our results only depend on the relation between *t*_{i} and the Bcd decay time τ=1/α (for *t*_{i}≪τ the gradient has not yet reached steady‐state). We find that fluctuations in *s*_{0} (Figure 2A) give rise to decreasing Δ*x* toward the posterior pole if the pre‐steady‐state gradient is read out (i.e. for small α). In contrast, for fluctuations in *D* or α, Δ*x* increases toward the posterior pole (Figure 2B and C). Thus, in the case of Bcd pre‐steady‐state decoding (bluish curves) and fluctuations both in the production and degradation or diffusion rates (Figure 2E), maximal precision of the Bcd gradient around mid‐embryo can arise naturally, as was indeed observed directly by Gregor *et al* (2007a) and indirectly by our analysis (cf. Figure 1). In contrast, for high decay rates (reddish curves) the gradient rapidly equilibrates such that fluctuations in the production rate give rise to uniform noise in Bcd concentrations, preventing the possibility of minimal noise in the central region. Interestingly, fluctuations in the nuclear trapping rate alone also yield higher precision around mid‐embryo (Figure 2D) if the Bcd gradient is decoded early in pre‐steady‐state. However, the position with highest precision is sensitive to the decoding time, shifting more toward the posterior pole for later decoding.

Our modeling approach provides a proof of principle that maximal precision at mid‐embryo can arise if the gradient is decoded before steady‐state. However, we cannot rule out that other mechanisms may yield such a pattern of precision even for a steady‐state Bcd gradient (e.g. in the wing disc, it was shown that cell‐to‐cell variability in the production, diffusion and degradation rates can yield higher precision around mid‐field, Bollenbach *et al*, 2008).

In our analysis, we assumed the classical French‐flag model (Wolpert, 1969) in which domain boundaries are determined from critical morphogen concentrations. This is a simplification and recent work by the Reinitz group (Jaeger *et al*, 2004a, 2004b, 2007; Jaeger and Reinitz, 2006; Surkova *et al*, 2008; Manu *et al*, 2009a, 2009b) showed that the gap gene expression domains are established by a highly dynamic process characterized by drifting domain boundaries. In our modeling, we neglected the gap gene dynamics and rather focused on how precision can be achieved already at the level of the Bcd gradient, as suggested by our data (as well as those reported in Gregor *et al*, 2007a; He *et al*, 2008). Thus, we assume that the positional information is transmitted relatively early to the gap genes, around cycle 10 (as suggested in Lucchetta *et al*, 2008). The precision and proper scaling of this read‐out is then maintained (and possibly refined) by the gap genes independently of Bcd, because of their cross‐ and auto‐regulation (Bergmann *et al*, 2007). Meanwhile, the Bcd gradient may continue to evolve and eventually decays at the onset of gastrulation (e.g. by activation of its PEST domain, Niessing *et al*, 1999).

### Modeling scaling

We next sought to extend our modeling framework to accommodate also scaling. To quantify scaling analytically at the morphogen level, we define a scaling coefficient *S*(*x*, *t*) for a position *x* and time *t*

Perfect scaling corresponds to *S*=1. In this case, fluctuations in embryo length, d*L*/*L*, are exactly compensated by fluctuations in position, d*x*/*x*, implying perfectly conserved proportions. We use the terms hypo‐ and hyper‐scaling to refer to *S*<1 and *S*>1, respectively. A position that hypo‐scales does not compensate enough for a change in embryo size, meaning that in a bigger embryo the absolute position is not shifted enough posteriorward to keep the correct proportions. In contrast, hyper‐scaling is the tendency to overcompensate for a change in embryo size.

Scaling is a property of the external gradient, which is then transmitted to the nuclear concentrations (cf. equation (7)). Assuming that the threshold concentration is fixed (implying d*M*=(∂*M*/∂*x*)d*x*+(∂*M*/∂*L*)d*L*=0) it follows from equation (2) that the scaling coefficient is

The above definition of scaling is generic and can be computed for any morphogen distribution *M*(*x*, *t*) with explicit dependence on *L*. Here, we consider a model in which the embryo length impacts the generation of the morphogen gradient through the nuclei density. Specifically, assuming that at decoding time all embryos have the same number of nuclei independent of their size (which agrees with the deterministic doubling of nuclei at each cycle), the nuclei density *N* depends on the embryo size like *N*∝*L*^{−n}, where *n*∈[1;3] (*n*=3 corresponds to a uniform distribution of nuclei, whereas *n*=2 is true if nuclei are distributed on a shell with a fixed width, Gregor *et al*, 2007b). In this scenario, we find that scaling is time dependent and position dependent (Figure 3B; see Supplementary Figure S9 for the dependence of scaling on *N*/*K* and *n*). Specifically, anterior domain boundaries hyper‐scale, in particular if decoding occurs relatively late, whereas posterior domains show very good scaling for a wide range of decoding times.

We also investigated three other models that give rise to scaling by coupling the nuclear degradation, cytoplasmic degradation or the morphogen production rate to the embryo length (cf. Supplementary Text S2). We find that their scaling behavior is qualitatively similar (in the sense that *S* decreases toward the posterior pole) but more sensitive to the Bcd decoding time (Supplementary Figure S8).

### The gap and pair‐rule gene expression domains scale with embryo size

The model we investigated above predicts that scaling should be *position* dependent rather than *gene* dependent and close to perfect except for anterior most domains. To test these two predictions experimentally on our collection of staining images, we adapted the continuous definition of the scaling coefficient *S*(*x*, *t*) in equation (2) for discrete measurements as follows:

where is the *estimated* slope from a linear regression *x*=α+β*L* of the domain positions *x*_{i} (with mean value ) onto their respective embryo sizes *L*_{i} (with mean value ) (see Figure 3A). In our data, the relative embryo size fluctuations were of the order of 10–15%. The results of our scaling analysis for Gt, Hb, Kr and Eve in embryos with single, double and quadruple *bcd* dosage are presented in Figure 3B. We observe that anterior domains indeed tend to hyper‐scale, whereas mid‐embryo and posterior domains show good scaling (*P*‐value <0.00021, cf. Supplementary Text S1; Supplementary Figure S5). Moreover, the magnitude of the scaling coefficient depends mainly on the position of the respective domain (boundary) rather than the associated gene. This is in good agreement with the model discussed above, whereas in the other models we investigated (cf. Supplementary Text S2) it is difficult to obtain close to perfect scaling both at mid‐embryo and toward the posterior pole.

### Perspectives

Precision and scaling are important to achieve robust pattern formation. Our new position‐dependent measures provide a unified quantification of scaling for both functional descriptions of gradient profiles derived from models (equation (2)) and experimentally determined expression domains (equation (4)). These measures will be useful to address scaling also in different systems like the wing‐disk. Importantly, our measures clearly identify perfect scaling (*S*=1) corresponding to the correct preservation of proportions, which is not necessarily equivalent to perfect correlation (see Supplementary Text S4 for a detailed discussion).

Interestingly, our experimental analysis using embryos with different *bcd* dosages showed that both precision and scaling are more *position* dependent than *gene* dependent, suggesting that the morphogen gradient itself has an important function to set robust domain boundaries. As our study showed, gradient formation of a single morphogen can naturally ensure maximal robustness only in a limited patterning domain. For Bcd, this appears to be mid‐embryo, where most of its targets are expressed. Outside of this domain, other systems and mechanisms may have evolved to cooperate in maintaining and potentially increasing robustness. Considering precision and scaling as position‐ (and possibly time‐) dependent features will allow to better characterize and develop models for these systems.

## Materials and methods

### Image analysis

We use an interface developed in MATLAB (Supplementary Figure S1) to extract and analyze expression profiles from the staining images of different Bcd target genes, yielding quantitative information on the positions of each protein domain (Bergmann *et al*, 2007). A total of 154 staining images were analyzed using a semi‐automated analysis tool where the position of the anterior and posterior poles were marked 50 times for each embryo. On the basis of this input, the software extracted a rectangular region from the image from which it generated protein concentration profiles and automatically determined the positional information of the gap and pair‐rule gene expression domains. Thus, for each embryo, its total length and the domain localizations were characterized by a mean value and a standard deviation computed from the 50 markings. The right boundary of the anterior Hb expression domain was defined as the position at which the decline in the Hb concentration is the steepest. Gt and Kr domains were described similarly by their left and right boundaries, whereas the position with maximal staining intensity characterized their center. As the Eve stripes are not so broad, we only extracted the position at which the Eve intensity was maximal. (All the staining images are available for download at http://www.unil.ch/cbg/morphogen/Images_paper.zip)

We also developed a second, almost fully automated analysis tool that only requires human input for selecting which of the two poles is the anterior one (while their position was detected automatically). This tool was used to analyze the expression domain positions of a set of 119 embryos, which included 78 (51%) from the first analysis. Precision and scaling results turned out to be very similar for the two analysis approaches, suggesting that measurement errors and inter‐embryo variation (e.g. because of different embryo orientations) can be safely neglected (Supplementary Figure S6). We also note that part of the observed fluctuations in embryo size may be induced by the fixation procedure. However, computing precision and scaling by considering embryos on different slides separately yielded similar results (Supplementary Figure S7), suggesting that fixation or other batch effects have no positional bias toward either of the embryo poles.

### Modeling

We consider a model in which Bcd is produced at the anterior pole with production rate *s*_{0}, diffuses according to a uniform diffusion constant *D*, is degraded at uniform rate α and is trapped and released by the nuclei in the embryo at rates *k*_{n}, *k*_{−n}, respectively. The corresponding coupled partial differential equations for the free morphogen concentration *M*(*x*, *t*) and the nuclear concentration *M*_{n}(*x*, *t*) read (Bergmann *et al*, 2007; Coppey *et al*, 2007; Gregor *et al*, 2007b):

where *N* is the nuclei density and *s*_{0}δ(*x*) the source term localized at *x*=0 (anterior pole). We consider zero‐flux (Neumann) boundary conditions at the posterior pole and that there is no morphogen at the initial time (*t*=0). Assuming that nucleo‐cytoplasmic exchanges occur rapidly (∂*M*_{n}/∂*t*∣_{M}≈0), as demonstrated by the experiments reported in Gregor *et al* (2007b), the following *effective* diffusion equation holds for the free morphogen (Bergmann *et al*, 2007; Coppey *et al*, 2007; Gregor *et al*, 2007b):

where =*D*/(1+(*N*/*K*)), and (with *K*=*k*_{−n}/*k*_{n}). At time *t*=150 min (corresponding to nuclear cycle 14), the morphogen profile is numerically fitted to an exponentially decaying gradient with length scale λ=0.2*L*. Thus, given a degradation rate, the diffusion constant is adjusted accordingly (see caption of Figure 2 for numerical values). We account for the presence of the nuclei by setting *N*/*K*=1 (i.e. to be specific, we assume that the probability that external Bcd is trapped equals the probability that nuclear Bcd is released; yet our qualitative results for precision are robust with respect to the exact choice of *N*/*K*).

According to the nuclear trapping model (equation (5)), the nuclear concentration *c*_{n}, which determines the precision and scaling of Bcd target genes, is given by (see also Coppey *et al*, 2007)

where *ν*_{n} is the nuclear volume. Importantly, according to equation (7) the fluctuations in the nuclear concentrations *c*_{n} are proportional to those in the external gradient, provided we can neglect variability in *ν*_{n} and *K*.

## Conflict of Interest

The authors declare that they have no conflict of interest.

## Supplementary Information

Precision and Scaling in Morphogen Gradient Read‐out

* Image analysis and control experiments:

We discuss the image analysis tools and limitations as well as the validity of our results.

* Scaling models

We discuss other models that could yield scaling of Bicoid as well as the advantages of our measure of scaling compared to correlations. [msb20107-sup-0001.pdf]

Dataset S1

Dataset S1. Numerical values for the quantitative analysis of precision (relative position x/L, imprecision s(x/L), error, gene name, co‐staining (for Eve domains), *bcd* dosage). [msb20107-sup-0002.zip]

Dataset S2

Dataset S2. Numerical values for the quantitative analysis of scaling (relative position x/L, scaling coefficient, 68% CI, gene name, co‐staining (for Eve domains), *bcd* dosage). [msb20107-sup-0003.zip]

## Acknowledgements

We thank Naama Barkai, Ben‐Zion Shilo and Zvi Tamari for careful reading of the manuscript, as well as all the members of the Bergmann *Computational Biology Group* and David Morton de Lachapelle for their support and comments during the whole project. In particular, we thank Toby Johnson, Karen Kapur and Sascha Dalessi for their help with the theoretical analyses and are very grateful to Zoltán Kutalik who implemented the fully automated version of the profile extraction tool. We are also grateful for financial support from the Giorgi‐Cavaglieri Foundation, the Swiss National Science Foundation (Grant #3100AO‐116323/1), SystemsX.ch (through the WingX project) and the European Framework Project 6 (through the AnEuploidy project).

## References

This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.

- Copyright © 2010 EMBO and Macmillan Publishers Limited