The role of spatial averaging in the precision of gene expression patterns

Embryonic development is driven by orderly, spatial patterns of gene expression that assign each cell in the embryo its particular fate. Experiments in recent years have vividly demonstrated that gene expression is often highly stochastic Elowitz et al. (2002). Yet, embryonic development is exceedingly precise. A vivid example is the Bicoid-Hunchback system in the early Drosophila embryo, which has become a paradigm for understanding the formation of spatial gene-expression patterns. Shortly after fertilization, the morphogen protein Bicoid (Bcd) forms an exponential concentration gradient along the anterior-posterior axis of the embryo, which provides positional information to the differentiating nuclei. One of the target genes of Bcd is hunchback (hb), which is expressed in the anterior half of the embryo. The posterior boundary of the hb expression domain is very sharp: by cell cycle 13, the position of the boundary varies only by about one nuclear spacing Houchmandzadeh et al. (2002); Gregor et al. (2007a); He et al. (2008). This precision is higher than the best achievable precision for a time-averaging based read-out mechanism of the Bcd gradient Gregor et al. (2007a). Intriguingly, the study of Gregor et al. revealed that the Hb concentrations in neighboring nuclei exhibit spatial correlations and the authors suggest that this implies a form of spatial averaging enhancing the precision of the posterior Hb boundary Gregor et al. (2007a). However, the mechanism for spatial averaging remained unclear.

In this manuscript, we analytically and numerically study the Bcd-Hb system. Our analysis reveals a simple, yet powerful mechanism for spatial averaging, which is based on the diffusion of Hb itself. We show analytically that Hb diffusion between neighboring nuclei reduces the super-Poissonian part of the noise in its concentration, with a factor which depends on the diffusion length of Hb and the dimensionality of the system. In essence, diffusion reduces the noise by washing out bursts in gene expression. This mechanism is generic, and applies not only to any developmental system, but also to any biochemical network in general. For example, if a signaling protein is activated at one end of the cell and then has to diffuse to another place to activate another system, e.g. the messenger ${\rm CheY}$ in bacterial chemotaxis, then our results show that the non-Poissonian noise in the activation of the signaling protein is washed out by diffusion; for this reason it may be beneficial to spatially separate the in- and output of a signaling pathway. Our analysis also reveals that while Hb diffusion reduces the noise, it also lessens the steepness of its expression boundary Gregor et al. (2007a). The interplay between these two antagonistic effects leads to an optimal diffusion constant of $D\simeq 0.1\operatorname{\mu m}^{2}/\operatorname{s}$ that maximizes the precision of the hb expression domain.

We consider a minimal model of the Bcd-Hb system, which contains the necessary ingredients for understanding the mechanism of spatial averaging. This model includes the stochastic and cooperative activation of Hb by Bcd in each nucleus, and the diffusion of Hb between neighboring nuclei. To model the cooperative binding of Bcd to the hb promoter, we assume that Bcd proteins can bind sequentially to five binding sites on the hb promoter. When all five sites are occupied, the promoter is active (activity $n=1$) and Hb proteins are produced stochastically with rate $\beta$; otherwise the promoter is inactive (activity $n=0$). Hb proteins are degraded stochastically with rate $\mu$. To obtain a lower bound on the precision of the hb expression domain, we assume that Bcd binds to the promoter with a diffusion limited on-rate $k_{\rm on}=4\pi\alpha D_{\rm B}/V=24\alpha D_{\rm B}/d^{3}=8.4\times 10^{-5}\operatorname{s}^{-1}$ where $\alpha\simeq 3\operatorname{nm}$ is the dimension of a binding site, $d\simeq 6.5\operatorname{\mu m}$ the diameter of a nucleus and $D_{\rm B}\simeq 0.32\operatorname{\mu m}^{2}/\operatorname{s}$ the Bcd diffusion constant Gregor et al. (2007b). Since the on-rate is assumed to be diffusion-limited, cooperativity of Hb activation is tuned via the Bcd off-rate $k_{\rm off}(j)=a/b^{j}$ which decreases with increasing number $j$ of promoter-bound Bcd. Using $a=12.3\operatorname{s}^{-1}$ and $b=6$ the average promoter activity as function of the number of Bcd molecules $B$ in the nucleus approximately follows the Hill-function

with $h=5$ and $K_{\rm eq}=690$. The Hill-coefficient $h$ was inferred from the relation between Hb and Bcd concentration, and the activation threshold $K_{\rm eq}$, where $\langle n(K_{\rm eq})\rangle=0.5$, is the average number of Bcd molecules, $\langle B\rangle$, in a nucleus at the Hb domain boundary Gregor et al. (2007a).

To describe the formation of the spatial Hb pattern, we place $N=64\times 64=4096$ nuclei on a square lattice with spacing $\ell=8.5\operatorname{\mu m}$ Gregor et al. (2007a). Hb diffuses over the lattice with diffusion constant $D$. With reflecting boundary conditions in one and periodic in the other direction, our model is a cylinder which is symmetric around the anterior-posterior $x$-axis. Differential activation of Hb is induced by the Bcd profile

which decreases exponentially with the distance $x$ from the anterior pole at $x=0$. The decay length is $\lambda_{\rm B}=119.5\operatorname{\mu m}$ Gregor et al. (2007a) and the threshold position for Hb activation is set to $x_{t}=0.5L$; $L=64\ell=544\operatorname{\mu m}$ is the length of the embryo. Diffusion of Bcd between nuclei induces fluctuations in the Bcd copy number on the time-scale $\tau_{\rm d}=\ell^{2}/(4D_{\rm B})\simeq 10^{2}\operatorname{s}$. Because $\tau_{\rm d}$ is much smaller than the time-scale for promoter binding, $k_{\rm on}^{-1}\simeq 10^{4}\operatorname{s}$, Bcd copy number fluctuations are efficiently averaged out by slow binding of Bcd to the promoter. We therefore assume that the total number of Bcd in a nucleus is constant and given by Eq. (2). The state of the system is described by the number of promoter bound Bcd and the number of Hb in every nucleus. The state changes by reactions inside a nucleus and diffusion of Hb between nuclei. The dynamics is described by a reaction-diffusion master equation which we solve numerically Elf and Ehrenberg (2004).

Fig. 1a shows simulation results for the average number of free Bcd molecules, $\langle B\rangle$, the average promoter activity, $\langle n\rangle$, and the average number of Hb molecules, $\langle H\rangle$, as function of the position $x/L$. Without Hb diffusion $\langle H(x)\rangle$ is proportional to $\langle n(x)\rangle$. For a finite Hb diffusion constant, however, the shape of the Hb profile is determined not only by $\langle n(x)\rangle$, but also by the diffusion length of Hb, $\lambda=\sqrt{D/\mu}$—with increasing $\lambda$ the profile becomes less steep. Fig. 1b shows instantaneous Hb profiles at different times. Fluctuations of promoter activity induce large fluctuations of $H$. These in turn lead to an uncertainty $\Delta x$ in the position at which $H$ crosses the threshold $\langle H(x_{t})\rangle=0.5\beta/\mu$.

We now consider the Hb boundary width when the Hb diffusion constant $D$ is zero. Fig. 2a shows that this decreases with increasing lifetime of the Hb protein. To understand this behavior, we note that the boundary width is, to a good approximation, given by

where $\sigma(x_{t})$ is the standard deviation of the Hb copy number and $|\langle H(x_{t})\rangle^{\prime}|$ is the magnitude of the Hb gradient at the boundary position $x_{t}$ Tostevin et al. (2007); Gregor et al. (2007a). When the Hb lifetime $\tau$ is much longer than the time-scale $k_{\rm off}(5)^{-1}$ of the promoter state fluctuations, then the noise in Hb copy number becomes Poissonian and the variance is given by the mean: $\sigma^{2}(x_{t})=\langle H(x_{t})\rangle=\langle n(x_{t})\rangle\beta/\mu=0.5\beta\tau$. On the other hand, the steepness of the boundary, ${|\langle H(x_{t})\rangle^{\prime}|}$, increases linearly with $\tau$ when $D=0$ and the synthesis rate $\beta$ is kept constant. Eq. (3) thus predicts that the boundary width $\Delta x$ decreases as $\tau^{-1/2}$ for large $\tau\gg k_{\rm off}(5)^{-1}$, which is indeed observed in the simulations (see Fig. 2a). Hence, the hb expression boundary could be made arbitrarily precise if the lifetime of the Hb protein could be increased indefinitely.

In practice the averaging time cannot be made arbitrarily long Gregor et al. (2007a). Ultimately, it is limited by the nuclear division time $\tau_{\rm nd}$, which at cell cycle 13 is roughly 30 minutes. Fig. 2a shows that for $\tau\simeq\tau_{\rm nd}$ and vanishing Hb diffusion constant $D=0$, the boundary width would be $\Delta x\simeq 4$ nuclei, which is larger than the precision measured experimentally Gregor et al. (2007a). The limited averaging time thus puts strong constraints on the precision that can be achieved via the mechanism of time-averaging alone.

Fig. 2a reveals, however, that the precision of the hb expression domain can be enhanced significantly by increasing the Hb diffusion constant to a finite value. This may seem surprising, since the steepness of the boundary, $|\langle H(x_{t})\rangle^{\prime}|$, decreases with increasing Hb diffusion constant (Fig. 1) and this—as Eq. (3) shows—tends to increase the boundary width. However, increasing the Hb diffusion constant also reduces the noise in the Hb copy number. This is the mechanism of spatial averaging, which we now study analytically.

To elucidate the mechanism of spatial averaging, we first analyze how the steepness of the Hb boundary, $|\langle H(x_{t})\rangle^{\prime}|$, depends on the Hb diffusion length $\lambda$, and then how the noise in the Hb copy number at the boundary, $\sigma(x_{t})$, depends on $\lambda$ (see Eq. (3)). To calculate the Hb profile, we note that a nucleus at a position $x_{i}$ produces Hb with an average rate $\beta\langle n(x_{i})\rangle$ (Fig. 3a). In one dimension, the steady state number of Hb molecules in a nucleus of size $\ell$ at $x$ is approximately

The total Hb profile is the sum of $\langle H(x)\rangle_{\rm ps}$ over all point sources $x_{i}=0,\ell,2\ell,\dots$. To do the sum, we approximate the activity $\langle n(x_{i})\rangle$ from Eqs. (1) and (2) by a piecewise linear function which behaves as $\langle n(x_{i})\rangle=0.5+m(x_{t}-x_{i})$ around the boundary for $x_{t}-1/2m<x_{i}<x_{t}+1/2m$ and which is unity below and zero above the boundary. The slope $-m=-h/4\lambda_{\rm B}$ is the derivative of $\langle n(x_{i})\rangle$ at $x_{i}=x_{t}$. At the boundary $x_{t}$ the slope of $\langle H(x)\rangle$ is

It increases linearly with the number of Hb molecules $\beta/\mu$ and decreases as the Hb diffusion length increases. Due to the rotational symmetry, the slope in the full, two-dimensional system is also given by Eq. (5). Fig. 3b plots Eq. (5) as function of $\lambda$ and compares the approximation with two-dimensional simulation results. It is seen that the agreement is very good.

Next, we study the effect of Hb diffusion on the variance in the Hb copy number, $\sigma^{2}$ (see Eq. (3)). For clarity, we first do the calculations for one-dimensional diffusion and extend our approximation to general dimensions afterwards. To compute $\sigma^{2}$, we exploit the observation made above that each nucleus acts as a point source of Hb (Fig. 3a), and that the expression of Hb in each nucleus is an independent stochastic process. We denote the variance of the Hb copy number in a nucleus of size $\ell$ at a distance $x$ from the point source as $\sigma_{\rm ps}^{2}(x)$. In a uniform space, the total Hb variance in a nucleus is the sum $\sigma^{2}=\sum_{i}\sigma_{\rm ps}^{2}(x_{i})$ over point sources at all $x_{i}$. To calculate $\sigma^{2}_{\rm ps}(x)$, we assume that the Hb dynamics is fast on the time-scale of promoter switching, such that the Hb concentration switches between zero and $\langle H(x)\rangle_{\rm ps}/\langle n\rangle$, with $\langle H(x)\rangle_{\rm ps}$ given by Eq. (4). In this limit,

where the second term is due to promoter switching with variance $\langle 1-n\rangle\langle n\rangle$, and the first term describes the Poisson noise coming from the production, diffusion and decay of Hb when the promoter is active.

Fig. 4a shows $\sigma^{2}_{\rm ps}(x)$ for different values of the diffusion length $\lambda$. For small $\lambda$, $\sigma_{\rm ps}^{2}(x)$ is large for small $x$, but reduces quickly with increasing $x$. Hence, the variance at a given nucleus is determined by few nuclei in the immediate neighborhood but each individual contribution is large. For increasing $\lambda$ the number of nuclei contributing to the variance increases, but the individual contributions of nearby nuclei are smaller; this is because diffusion washes out the bursts of Hb production at the source. To see which of the two opposing effects dominates, we sum Eq. (6) over $x_{i}$ to obtain the total variance:

Compared to the variance in an isolated nucleus, the overall effect of diffusion in a one-dimensional, uniform space is a reduction of the non-Poissonian noise by $(4\lambda/\ell)^{-1}$. The Poissonian part cannot be reduced because diffusion itself is Poissonian. Indeed, spatial averaging can reduce the effect of noise in gene expression, but only if this is super-Poissonian. Simulations also allow us to compute $\sigma_{\rm ps}^{2}(x)$ for slower Hb dynamics. In this case, the non-Poissonian part of the variance is smaller. However, the qualitative dependence of $x$ remains the same: the non-Poissonian part is washed out by diffusion over a distance proportional to $\lambda$ (inset Fig. 4a).

The finite range of non-Poissonian fluctuations implies the following simple model. $N$ nuclei contribute equally to a well-stirred reservoir containing $M$ Hb molecules, with average $\langle M\rangle=N\langle H\rangle$ and variance $\sigma^{2}_{M}=N\sigma_{0}^{2}$, where $\sigma^{2}_{0}$ is the noise of $H$ in a nucleus if there was no Hb diffusion. Each nucleus then samples Hb molecules in a binomial fashion from this reservoir. The noise in each nucleus is formally given by $\sigma^{2}=E[V[H|M]]+V[E[H|M]]$, where $E[X]$ denotes the expectation and $V[X]$ the variance of $X$. The first term describes the noise in the sampling process, and is given by $E[V[H|M]]=\langle H\rangle(1-N^{-1})$. The second term describes the variance in $H$ due to fluctuations in the reservoir, and is given by $V[E[H|M]]=\sigma_{0}^{2}/N$. Hence,

This simple expression elucidates that it is the super-Poissonian part of the noise, $\sigma_{0}^{2}-\langle H\rangle$, which is reduced by spatial averaging. Since $N$ scales with the dimensionality $d$ as $N\sim\lambda^{d}$, this reduction becomes more efficient in higher dimensions. Interestingly, this expression also reveals that if hb expression would be sub-Poissonian Pedraza and Paulsson (2008), diffusion would increase the noise.

Fig. 4b shows the non-Poissonian part of the variance in a uniform, two-dimensional system as function of the diffusion length; the exact, numerical result is compared to the prediction of Eq. (8) with $N$ chosen to be the number of nuclei within a distance $2\lambda$ of a given nucleus along the edges of a square lattice, $N(\lambda)=8(\lambda/\ell)^{2}$. For sufficiently large $\lambda$ for which $N(\lambda)>1$, the approximation is excellent. Fig. 4b also compares the results for the uniform system to those of the full, non-uniform system in which a Bcd gradient activates Hb. Clearly, the well-stirred approximation describes the variance at the Hb boundary very well; only for $\lambda>5\ell$ is the variance at the boundary significantly influenced by nuclei far away from the boundary with smaller variance.

Finally, we can combine Eqs. (5) and (8) with Eq. (3) to predict the boundary width $\Delta x$. Fig. 2b shows that the prediction agrees very well with the simulations. This figure also shows that the two antagonistic effects of Hb diffusion—reducing the slope but also the variance of the Hb concentration—lead to a diffusion constant that optimizes the boundary width for a fixed Hb lifetime. Interestingly, this minimal width is less than one nuclear spacing, as found experimentally Gregor et al. (2007a). This suggests that for achieving the necessary precision, mechanisms based on multiple gradients or interactions between multiple Bcd targets are not required, although these may provide robustness against embryo-to-embryo variations Howard and ten Wolde (2005); Manu et al. (2009a, b). Since the exact mechanism for Hb transport is unknown, it is tempting to speculate that the precision achieved via spatial averaging could be increased further by separating the two effects of Hb diffusion using an anisotropic transport of Hb: slow transport along the anterior-posterior axis would allow for steep spatial profiles, while fast diffusion along the perpendicular direction would allow for effective spatial averaging. The benefit of anisotropic transport could imply the existence of active mechanisms for Hb transport.

Spatial averaging can only be beneficial when the noise in Hb production has a super-Poissonian component. The observation of spatial correlations in Hb by Gregor et al. indicates that this is the case Gregor et al. (2007a). The question arises whether bursts in gene expression are inevitable and spatial averaging a prerequisite for achieving a precise hb expression domain. The Hb lifetime is limited by the nuclear division time, while the time-scale for promoter state fluctuations is limited by the diffusion of Bcd. The magnitude of the bursts could be reduced by reducing the promoter strength, but this would lower the Hb copy number and hence the steepness of the Hb boundary; the net result would, in fact, be a decrease of its precision. In a model with explicit hb mRNA, the effect of promoter state fluctuations could be alleviated by reducing the transcription rate. However, to achieve a sufficiently steep and precise Hb boundary with a low number of hb mRNA, the translation rate had to be increased such that Hb production becomes super-Poissonian again. It thus appears that spatial averaging is a fundamental mechanism for generating precise patterns of gene expression. Clearly, it will be of interest to study how other modes of transport, e.g. sub-diffusion or active transport, affect the mechanism of spatial averaging.

We thank Filipe Tostevin for a critical reading of the manuscript. This work is supported by FOM/NWO (TE, PRtW) and The Royal Society (MH).