Morphogen Profiles Can Be Optimised to Buffer Against Noise

Morphogen profiles play a vital role in biology by specifying position in embryonic development. However, the factors that influence the shape of a morphogen profile remain poorly understood. Since morphogens should provide precise positional information, one significant factor is the robustness of the profile to noise. We compare three classes of morphogen profiles (linear, exponential, algebraic) to see which is most precise when subject to both external embryo-to-embryo fluctuations and internal fluctuations due to intrinsically random processes such as diffusion. We find that both the kinetic parameters and the overall gradient shape (e.g. exponential versus algebraic) can be optimised to generate maximally precise positional information.

Morphogens are signaling molecules which play a vital role in biological development by inducing responses in a concentration-dependent manner [1]. In the standard model of morphogen gradients, morphogen proteins originate from a localized source, diffuse and are degraded, setting up a concentration gradient across the system. This gradient can control patterns of gene expression, where, for example, a gene is switched on when the concentration is above a certain fixed threshold, but is off otherwise. In this work, we focus on morphogen gradients specifying boundaries of gene expression at fixed absolute distances from the morphogen source, a scenario that is frequently realised in developmental biology.
Developmental systems need to be robust to sources of noise in order to generate precise patterns of gene expression. Here, we address a simple question: in the presence of noise, which morphogen profile is most precise in specifying the positions of gene expression boundaries? In principle, any spatially non-uniform profile could be used to position gene expression boundaries; our goal is to understand which profiles might be preferred. Previous theoretical approaches have predominantly examined robustness of morphogen profiles to embryo-to-embryo fluctuations in the morphogen production rate [2,3,4], and analysis suggests that algebraic profiles are most precise [2]. However, even in (hypothetical) embryos with no embryo-to-embryo fluctuations, there would still be variation in positional information due to internal fluctuations in morphogen production, diffusion and degradation [5]. Such fluctuations impose limits on the precision of biochemical signaling [6,7,8], and could in principle alter the shape of the profile with the best precision. Internal fluctuations are particularly large if the morphogen is, directly or indirectly, a transcription factor. In that case, the arrival of a morphogen molecule at the nanometer-scale DNA binding sites on the target genes will be a rare, stochastic event [8,9]. Moreover, internal fluctuations in the morphogen production rate may also play an important role [10].
Recent experiments in the fruit fly Drosophila melanogaster have quantitatively studied the morphogen proteins Bicoid (Bcd) [9,11,12], Decapentaplegic (Dpp) [13] and Wingless [13]. Interestingly, the observed profiles are exponentially decaying in all cases. In order to better understand this finding, we compare three classes of morphogen profiles (linear, exponential, algebraic) to see which is most precise when subjected to the combined effects of both external and internal fluctuations. We find that the kinetic parameters describing morphogen profiles can be optimised to buffer against the combined effects of both sources of noise. By comparing optimised profiles, we then see that the overall shape of the profile (e.g. exponential versus algebraic) can also be optimised. Exponential profiles frequently emerge as the best compromise: such profiles are not particularly robust to either external or internal fluctuations taken singly, but when both types of fluctuation are taken together, exponential profiles can be the most precise. We therefore propose a simple design principle for morphogen profiles, namely that evolution has selected gradients with optimal robustness to the combined effects of embryo-toembryo and internal noise. Depending on which source of noise is most important, qualitatively different morphogen profiles will be selected. Given that very high positional precision can be achieved by morphogens (e.g. a few percent of embryo length in the Bcd system), it seems plausible that optimisation may well be exploited by evolution. However, there may still be other constraints on morphogen systems, for example, on the information capacity of the signaling network [14,15].

Models.
We consider three representative models of morphogen gradients: a freely diffusing morphogen with a source and a sink [16]; diffusion with linear decay; and diffusion with quadratic decay [2]. Why study these three particular models? Diffusion with linear decay leads to exponential morphogen profiles that have been observed experimentally. Non-linear decay mechanisms leading to algebraic profiles are robust to external morphogen production fluctuations [2] and we study quadratic decay (decay via dimerisation) as a representative example [2]. Bone Morphogenetic Protein (BMP) in Xenopus embryos is a possible candidate for a morphogen shaped by such effective non-linear degradation [17]. Finally, the sourcesink model generating a linear profile provides insight into gradients generated by diffusion with localised degradation [16]. This is realisable when the morphogen degradation factors are tightly localised, such as for retinoic acid along the anterio-posterior axis of mice embryos [18].
A three-dimensional system is considered with a planar source of morphogen at x = 0 that produces a flux of proteins (J) which diffuse (D) through the system (length L). In the absence of fluctuations the morphogen concentration is given by the reaction-diffusion equation where ψ = µ[s −1 ] (n = 1) and ψ = α[µm 3 s −1 ] (n = 2). The boundary conditions are D∂ x ρ(x)| x=0 + J = 0 and, for n = 1, 2, D∂ x ρ(x)| x=L = 0. The source-sink model corresponds to ψ = 0 and a sink at the boundary x = L (i.e. ρ(L) = 0). The steady-state solutions to (1) depend only on the linear distance from the source in the x-direction: where ρ s−s,lin,quad correspond to the source-sink, linear decay and quadratic decay models respectively and λ = D/µ, External fluctuations. We first investigate robustness to external fluctuations solely in J as this is believed to be an experimentally relevant scenario [9]. The position x T where the concentration passes through a threshold (ρ T ) varies as a result of embryo-to-embryo fluctuations (δJ) in the morphogen production rate. Keeping the threshold concentration fixed, and expanding around x T to leading order, the width due to external fluctuations is given by . denotes ensemble averaging) are the variations due to external embryo-toembryo fluctuations, and ρ ′ (x T ) = ∂ x ρ| x=xT . The widths for the three models are Internal fluctuations. Within an embryo internal fluctuations are also an important source of noise. We consider morphogen production, diffusion and (if appropriate) degradation to be stochastic processes. The morphogen concentration is sampled in a region of linear size ∆x corresponding to the size of the binding site at target genes. Note that incorporating details of the binding process is unlikely to alter our results [7]. Internal fluctuations in particle density within the sampling volume again cause the position where the gradient passes through ρ T to vary, leading to imprecision in the positional information provided by the gradient. The width due to internal fluctuations ω is again given by ω ≈ ∆ρ(x T )/| ρ ′ (x T ) | [8], but where ∆ρ is now due to internal fluctuations. In the source-sink and linear decay models the statistics of the particle number n(x) due to internal fluctuations are Poissonian, since both models are linear and morphogen production, diffusion and degradation (if present) are all Poisson processes [8,19,20]. We present numerical results below demonstrating that the particle number statistics in the quadratic decay model are also effectively described by a Poisson distribution. This result is not obvious, as non-linear decay processes could, in principle, generate non-Poisson fluctuations. In all cases, we assume that diffusion of morphogen proteins is a purely three dimensional process (there could be one dimensional sliding along DNA but this leads to similar fluctuations as three dimensional diffusion [21]).
Simulations. Stochastic simulations were performed on a three-dimensional lattice containing discrete particles with particle injection on the plane x = 0, diffusion and (if appropriate) degradation. For the quadratic decay model, simulations were performed using a range of parameter values for J, D, P (probability of two particles degrading within a single time step given that they occupy the same lattice site) and system size L. In Fig. 1A we demonstrate that concentrations measured in our simulations agree well with (4) (the finite size ef-fects at large x do not alter our conclusions [22]). To confirm that particle number fluctuations are described by Poisson statistics for the quadratic decay model we calculated f (x) = [ n(x) 2 − n(x) 2 − n(x) ]/ n(x) , where each simulation was averaged for 5000s (Fig. 1B). In three dimensions diffusion is efficient enough to prevent the build-up of non-Poissonian correlations resulting from the quadratic decay reactions. This result is robust to parameter variations in our simulations (data not shown).
Time/spatial averaging. We now consider the effects of time and spatial averaging [8,9], which act to reduce ω. Time-averaging is performed by the down-stream signalprocessing network, where the timescale is typically given by the lifetimes of the transcripts/proteins of the target gene. Over an averaging period τ , there can, at most, be N τ = τ /τ ind independent readings of the morphogen concentration which reduce the measurement width by a factor ∼ N −1/2 τ . Intuitively, τ ind ∼ (∆x) 2 /D 0 , the typical timescale for the sampling volume to empty and then be refilled by new protein particles via diffusion [6,8] (D 0 is 'local' diffusion responsible for movement into the sampling volume, which may not equal D, the bulk diffusion). A constant associated with time averaging is found numerically (see below). We also include spatial averaging, motivated by recent experiments on the Bcdhb signaling pathway [9]. By averaging over a number of different nuclei/cells, N spat , the effects of internal fluctuations can be further reduced by a factor ∼ N −1/2 spat (up to a constant of O(1), which we have verified numerically does not alter our conclusions). From [9], N spat ≈ CD 0 τ where C is a constant which depends on the particular arrangement of nuclei/cells.
Width due to internal fluctuations. For the three models we find that the widths ω due to internal fluctuations are: In Fig. 1C, we compare (7) with simulation results for a range of parameter values. We find good agreement between the two approaches, where k 2 = 0.56 ± 0.06 is a fitting parameter. A similar procedure for the other two models yields k 0 = 0.53 ± 0.07 and k 1 = 0.60 ± 0.05. Optimizing kinetic parameters. Importantly, the underlying kinetic parameters in the linear and quadratic decay models can be optimised to minimise ω lin and ω quad at the threshold position x T . The existence of optimal decay rates is a general feature of morphogen gradients and occurs for all decay exponents n > 0 [22]. Maximum precision is achieved when λ = x T [8] and x 0 = 3x T . For the non-linear decay model we verified numerically the existence of such an optimal decay rate, see Fig. 1D. In general, in the experimentally relevant range x T ∼ L/2, x 0 ∼ λ, we find ω s−s < ω lin < ω quad . Although both the morphogen density and slope affect ω, the value of the slope is the more important: the sourcesink model has the steepest, and the quadratic model the least steep, profile, thereby generating the above ordering. Comparing this inequality with our earlier result W quad < W lin < W s−s , the ordering is reversed, so that, for example, the quadratic decay model performs best on external fluctuations but worst on internal fluctuations.
The robustness of our three models to the combined effects of internal and external fluctuations can now be compared. Since internal and external fluctuations are statistically independent, the total width is given by ǫ = √ ω 2 + W 2 . We find that the minimum in ǫ lin with respect to λ (and in ǫ quad with respect to x 0 ) becomes more pronounced than is the case with internal noise alone (see Fig. 2). As a result, the penalty in terms of reduced precision when using non-optimized parameters is now more severe than in the internal noise-only case. Note that the optimal values of λ and x 0 are reduced compared to their values when considering internal noise alone, see Fig. 2.
Optimizing morphogen profile shape. We investigate which model gives the most precise positional information when subjected to embryo-to-embryo fluctuations (parameterised by δJ/J) and internal fluctuations (parameterised by the averaging time, τ ). For each δJ/J and τ , we compute the optimal decay rates for the linear and quadratic decay models as above. We then build up an effective phase diagram in the τ − δJ/J plane, determining the most precise morphogen profile for particular levels of external and internal noise, see Fig. 3A. The parameters values used are similar to those found from experiments on hb expression in the Drosophila embryo [9,11,12,24]. J is kept constant between the models, representing fixed resource expenditure. The colourmap is defined as (ǫ lin − ǫ a )/ǫ lin × 100% where ǫ a is the smaller of ǫ quad and ǫ s−s . The percentage change is negative when the linear decay profile is most precise. For small times or small δJ/J, profiles from the sinksource model are preferred as this model is best able to buffer the dominant internal fluctuations. For large δJ/J or large τ quadratic decay profiles are selected, as now the quadratic model is now able to best buffer the dominant external fluctuations. At intermediate values of δJ/J and τ , when W and ω are similar, exponential profiles are preferred as they provide the best compromise between robustness to both internal and external noise. From such diagrams we can build up a qualitative understanding of why particular morphogen profiles are selected depending upon the dominant sources of noise. We emphasise that our methodology is general and not dependent on the specific parameter values used above. For example, if D = D 0 then quadratic decay would be favoured in a wider region of the phase space, but there would still be conditions when linear decay, or the sourcesink model, would deliver higher precision. If we choose a general decay exponent n, then, for a fixed δJ/J, we find a qualitatively similar picture with small n favoured at short averaging times, with larger n selected at longer times [22]. Multiple targets. In many developmental systems morphogens are directly involved in the expression of several target genes at different positions throughout the system (e.g. orthodenticle (x T ≈ 125µm) [25] and hb (x T ≈ 240µm) are both targets of Bcd). We consider optimising the morphogen profiles at one expression boundary and then compare how precisely the profiles determine a second boundary at a different position (keeping δJ/J fixed). Using the same optimised morphogen profiles as above we find that exponential profiles, unlike the other gradients, lead to the best precision at both short and long distances for short times (τ < 15mins), Fig. 3B. This flexibility is a potential explanation for the widespread use of exponential profiles in developmental systems.
We have applied our analysis to simple (though still experimentally motivated) models of gradient formation, yielding linear, exponential and algebraic profiles. Intriguingly, the best characterized morphogens Bcd, Dpp and Wingless in Drosophila all have exponential profiles and their decay lengths are significantly less than their respective values of x T [12,13]. This is consistent with our conclusion that decay lengths adapt to buffer the combined effects of internal and external fluctuations and that exponential profiles perform well when buffering combined internal-external fluctuations for a range of x T (since all three morphogens have multiple target genes). However, more complex processes may also be involved in the formation of these gradients, including, for example, transcytosis, pre-steady-state measurement or mRNA gradients [3,26,27]. In some circumstances, morphogen systems are able to scale with embryo length [9,12,28]. Nevertheless, once these additional effects are better characterized, our concepts can still be readily applied. Attaining maximal robustness to the combined effects of internal and external noise may be a powerful unifying principle in understanding the fundamental design of morphogen systems.
We thank Enrico Coen and Richard Morris for insightful comments. MH thanks The Royal Society for funding.