## Abstract

The patterning of many developing tissues is organized by morphogens. Genetic and environmental perturbations of gene expression, protein synthesis and ligand binding are among the sources of unreliability that limit the accuracy and precision of morphogen-mediated patterning. While it has been found that the robustness of morphogen gradients to the perturbation of morphogen synthesis can be enhanced by particular mechanisms, how such mechanisms affect robustness to other perturbations, such as to receptor synthesis for the same morphogen, has been little explored. Here, we investigate the interplay between the robustness of patterning to the changes in receptor synthesis and morphogen synthesis and to the effects of cell-to-cell variability. Our analysis elucidates the trade-offs and constraints that arise as a result of achieving these three performance objectives simultaneously in the context of simple, steady-state morphogen gradients formed by diffusion and receptor-mediated uptake. Analysis of the interdependence between length scales of patterning and these performance objectives reveals several potential mechanisms for mitigating such trade-offs and constraints. One involves downregulation of receptor synthesis in the morphogen source, while another involves the presence of non-signalling cell-surface morphogen-binding molecules. Both of these mechanisms occur in *Drosophila* wing discs during their patterning. We computationally elucidate how these mechanisms improve the robustness and precision of morphogen-mediated patterning.

## 1. Introduction

Embryonic development and patterning are commonly orchestrated by secreted signalling molecules known as morphogens. Morphogens diffuse away from their site of production, bind receptors on cell membranes and may interact with other diffusive or fixed molecules. They form stable gradients from which cells obtain positional information that often dictates cell fate. For example, in *Drosophila melanogaster* wing imaginal discs [1–7], Decapentaplegic (Dpp), a bone morphogenetic protein (BMP), functions as a long-range morphogen to specify pattern and control growth. Dpp is secreted from a strip of cells next to the anterior–posterior compartment boundary and diffuses into both compartments to form gradients. Other morphogens are present in the same system, such as Wingless (Wg), which regulates dorsal–ventral patterning [8,9], and Hedgehog (Hh), which patterns the central part of the disc and specifies the location where Dpp is produced [10].

Robustness of a morphogen signal, which is often referred to as the capability of resisting or buffering fluctuations in gene dosage or environmental conditions, has been observed in many morphogen systems [11–16]. One common approach for studying genetic and environmental perturbation is through heterozygous mutations in genes that encode morphogens. Observations of where development tolerates heterozygosity for most of these genes suggest that morphogen systems are highly robust.

The robustness to perturbation in morphogen synthesis is one of the major focuses in most studies on the subject [14,17–21]. One generic strategy to make gradients robust to perturbation in morphogen synthesis is through ‘self-enhanced ligand degradation’ [14]. Morphogen receptor synthesis is regulated by morphogen signalling in such a way that morphogen degradation increases with the strength of signalling. Several morphogen systems, including Hh and Wg gradients in patterning of *Drosophila* wing disc [14] and retinoic acid in hindbrain patterning of zebrafish [22], exhibit features consistent with this strategy. In this strategy, the robustness is enhanced by decreasing the length scale of the signal gradient near the source [17–20]. The improvement in the robustness can be shown through reducing the sensitivity coefficient of the positional information to morphogen synthesis, which is measured through the fold change in patterning position with respect to any fold change in morphogen synthesis.

Stochastic events, such as morphogen–receptor binding processes and randomness in morphogen movement, particularly at low morphogen levels, produce uncertainty in morphogen signals and affect the precision of patterning and gene expression borders [17,23–25]. The ability of morphogens to specify the tissue boundary locations of patterns depends on morphogen levels. The transition width, which is used to estimate the size of the region of ‘salt-and-pepper’ cell responses, provides a measurement for the sharpness of the borders separating cells of different fates driven by morphogens [23]. Studies of the transition width show that the precision of the sharpness of the cell response borders decreases from the source, and the positional information may reach a maximal precision within a few cells near the morphogen source [17,24]. To study how the kinetic parameters affect the maximal precision, three relevant classes of morphogen models (linear, exponential and algebraic) were examined to investigate which of them is most precise in producing a pattern when subject to both external and intrinsic stochastic events [26].

When the sensitivity coefficient of positional information to morphogen synthesis, which decreases with the distance from the source of the morphogen, is considered together with the transition width, which increases with distance, a trade-off is clearly seen between these two performance objectives. Such a trade-off can result in a limitation in the size of the spatial region that can display both a desirable sensitivity coefficient and a short transition width [17]. The relative width of such a patterning region is a function of the length scale near the source and becomes zero when the length scale and the maximum pattern width are longer than certain levels. Because the length scale of a morphogen gradient determines the range of morphogen spreading, the limitation in the length scale of the morphogen due to such a trade-off restricts the ability to have long-range patterning of morphogens.

The trade-off between the two performance objectives potentially becomes worse if more performance objectives are considered. While we know that receptors that bind morphogens to transmit the morphogen signals to the individual cells clearly play an important role in patterning, the sensitivity coefficient of positional information to receptor synthesis is much less explored. Experimental evidence in the *Drosophila* wing disc reveals complex roles of the receptor in morphogen-mediated patterning, in which Dfz2 (Wg receptor) broadens the range of Wg action by stabilizing it from degradation [8] and overexpressing Tkv (Dpp receptor) shrinks the Dpp gradient [27]. Both experimental [28] and theoretical [29] studies suggest that downregulation and overexpression of Tkv have crucial effects on the shape and robustness of morphogen activity.

Here, we investigate multiple performance objectives simultaneously and obtain the constraints that may arise from achieving them collectively. We consider a system of morphogens that are synthesized in a local region of tissue and that diffuse and bind with non-diffusible receptors to initiate signal transductions. By analysing the robustness to both morphogen and receptor syntheses, along with the noise-induced cell-to-cell variability and its impact on patterning borders, we find trade-offs arising among these three performance objectives. In particular, the three performance objectives severely restrict the spatial range over which accurate and precise patterning may be achieved by morphogen gradients without complex regulations.

We first study a model consisting of only morphogens and their receptors. We derive an explicit analytical formula for the sensitivity coefficient of positional information to any given parameters in the model. Through a further analytical study on the transition width, we then estimate the range of the length scale for a robust and precise patterning region that is constrained by the three performance objectives. These analytical estimates, along with direct simulations, show that it is impossible for a morphogen system without any regulation to form patterns that are robust both to morphogen and receptor syntheses and to a ‘reasonable size’ of the transition width. The derived explicit relationships for the constraints suggest several mechanisms for relaxing the trade-offs and improving the balances among the three performance objectives.

The constraints can be partially alleviated through mechanisms such as receptor downregulation in morphogen synthesis regions [20,30,31] and the presence of non-signalling binding sites [20,32–34], as observed for Dpp gradients in the *Drosophila* wing disc. Finally, through both analytical estimates and computational exploration of the model, we investigate how each mechanism relaxes the trade-offs and widens the range of the length scale of the signal gradient, leading to formation of a long-range robust and precise patterning region. The approach applied here may be useful for explaining some of the complex regulatory architectures of morphogen systems.

## 2. Results

### 2.1. Mathematical model

We model morphogens, receptors, their bindings and other regulations through reaction–diffusion equations. This approach has been generally applied for studying the formation and robustness of morphogen gradients [17–20,26]. We use the experimental observations of the wing disc in *Drosophila* [28,32,35,36] to constrain and calibrate our model, which may be applied to other morphogen-mediated patterning systems.

In the model, morphogens are synthesized in a local region of tissue and diffuse and bind with non-diffusible receptors to form a gradient of morphogen–receptor complexes, which generate intracellular signals. These complexes govern patterning by placing boundaries between cell fates at particular thresholds of signal concentration. Here, we simplify the problem by representing the patterned tissue as a one-dimensional domain (figure 1). The midpoint of the morphogen production region is denoted *x* = −*x*_{p}; the boundary between the production region and the patterning region is *x* = 0 and the edge of the pattering region is *x* = *x*_{max}. The variables in the model are the concentrations of morphogens ([*L*]), receptors ([*R*]), and morphogen–receptor complexes ([*LR*]), time (*t*) and distance (*x*) from the morphogen source
2.1
2.2
2.3

In this scheme, morphogens bind to receptors, and their binding is assumed to be governed by a second-order binding rate constant *α*_{R}. Morphogen–receptor complexes dissociate according to a first-order rate constant *β*_{R}. The constants *d*_{R} and *d*_{LR} are first-order degradation rate constants for the receptors and morphogen–receptor complexes, respectively. *D*_{L} is the diffusion coefficient of the free morphogen. *V*_{L}(*x*) and *v*_{R} are the synthesis rates of the morphogen and receptor, respectively. *V*_{L} has a constant value *v*_{L} in the production region and is zero elsewhere. *v*_{R} is taken to be a constant in the entire region (although we will relax this assumption later). We assume that the morphogen gradient is symmetric with respect to *x* = *−x*_{p} and that the concentration of morphogens is effectively zero at *x* = *x*_{max}, i.e. that *x* = *x*_{max} is chosen to be sufficiently far from the morphogen source such that few morphogens are present. Thus, the boundary conditions are
2.4

For analytical treatments, we typically choose *x*_{max} = *∞*. This minimal morphogen gradient model fits many features of the Dpp and Wg gradients of the *Drosophila* wing disc [37].

The formation of morphogen gradients typically requires only a few hours, which is relatively short compared with the timescale of tissue development, which is usually longer than a day. In other words, the morphogen gradients quickly ‘relax’ to their steady states whenever the tissue changes its size (on a much slower timescale). Most of the previous studies on developmental patterning used the steady-state assumption [6,13,19]. Hence, we assume that morphogen gradients are realized by cells at the steady state. We thus take the intracellular signal, Sig, produced by the morphogen gradient to be equal to the steady-state concentration of the morphogen–receptor complex: 2.5For the readers' convenience, the definitions of all the parameters and variables used in this paper are listed in table 1.

### 2.2. Sensitivity coefficients of the positional information to synthesis rates

The robustness of systems to external or internal variations is often quantified by the sensitivity coefficient, which corresponds to the fold change in an output of interest in response to a given fold change in a particular input [38]. In general, for any input *a* and output *b*, we define the sensitivity coefficient as
2.6

To express the sensitivity coefficient of the output of a morphogen gradient to any parameter *p*, we introduce a function *x*(Sig, *p*) for the spatial location at which the signal is at level Sig, given the parameter value *p*. Then, the sensitivity coefficient of the positional information to *p* is *S _{x}*

_{,p}.

Most of the existing experimental studies have focused on the sensitivity to morphogen or receptor synthesis independently [8,14,27]. The majority of the experimental results are for Dpp/BMP systems. The sensitivity to other processes has hardly been explored in experiments or in theory. Hence, we focus on the case where *p* = *v*_{L}, the morphogen synthesis rate, or *p* = *v*_{R}, the receptor synthesis rate.

We also introduce the notion of the length scale *λ* of a morphogen gradient, which quantifies the spatial range of action of morphogens, as determined by the balance between morphogen diffusion and clearance. As these processes are not necessarily spatially uniform, the length scale is defined in a local manner as the ratio of the signal to the slope of the signal gradient [14,17]
2.7

The negative sign is used because Sig is a decreasing function and we consider *λ* to be a positive value.

Given the definitions of the sensitivity coefficient (2.6) and the length scale (2.7), we can write the sensitivity coefficient of positional information in terms of any parameter *p* at a signal level Sig, i.e. *S _{x}*

_{,p}(Sig), in terms of

*x*and

*λ*. First, by the definition of

*λ*, we have 2.8where Sig

_{0}(

*p*) is the level of signal adjacent to the source of the morphogen, Sig(0,

*p*). The derivative of

*x*with respect to

*p*is 2.9which leads to 2.10Then, the sensitivity coefficient

*S*

_{x}_{,p}(Sig

_{1}) can be written as 2.11

When *p* = *v*_{L}, for the case of the sensitivity coefficient to the morphogen synthesis rate (figure 2*a*,*b*), the first term of the right-hand side in (2.11) is zero for all *x* > 0 because there is no morphogen production for *x* > 0, so the length scale is independent of *v*_{L}, i.e. . Consequently, the sensitivity coefficient of positional information to the morphogen synthesis rate will be proportional to the length scale of the signal gradient evaluated at *x* = 0, a result consistent with previous studies [14,17].

When *p* = *v*_{R}, for the case of the sensitivity coefficient to the receptor synthesis rate, the first term of the right-hand side in (2.11) is negative because the length scale *λ* decreases as more receptors are produced; the second term of the right-hand side in (2.11) is positive because Sig_{0} increases with the expression of receptors. The sensitivity coefficient drops to zero near the source when the second term decreases to the threshold, which is equal to the magnitude of the first term. In the region far from the source, the second term continues to decrease to zero. The first term becomes dominant in the right-hand side, and the sensitivity finally approaches a constant value (approx. 0.5, shown in figure 2*d*). This non-monotonic behaviour of the sensitivity coefficient to the receptor synthesis rate (figure 2*c*,*d*) is different from that of the morphogen synthesis rate, for which the sensitivity coefficient monotonically decreases away from the source (figure 2*a*,*b*).

This result agrees with experimental findings in *Drosophila* wing discs that overexpression of the Dpp receptor Tkv reduces the spatial range of Dpp action, while heterozygous mutations of Tkv lead to increases in the range [27,32,39].

### 2.3. Trade-offs among sensitivity coefficients to morphogen and receptor synthesis rates and the transition width

The transition width measures the width of a ‘salt-and-pepper’ transition zone in a signal due to cell-to-cell variability and stochastic variation in morphogen concentration [17,24]. The transition width, denoted by *W*(*x*), is estimated as
2.12where CV(*x*) is the coefficient of variation (standard deviation divided by the mean) of the amount of signal molecules at position *x*. Similar to previous studies [17,26], we assume that the signal in a cell is a Poisson-distributed random variable [24]. Then, *W*(*x*) can be estimated by
2.13where *J* (with unit μM^{−1}) is defined as the number of signalling molecules per micromolar concentration to express the concentration in terms of the number of molecules per cell. The signal is estimated through the steady-state solution of the deterministic equations (2.1)–(2.5). Compared with direct stochastic simulations, this approximation can save us computational cost, particularly when a large number of parameter sets are explored.

Intuitively, the length scale in (2.12) is used to transform the coefficient of variation of signal molecules to a range of positional fluctuation because the length scale is the ratio of the positional difference (positional fluctuation) to the corresponding relative change in the signal level (coefficient of variation). For a system to achieve good robustness and precision of patterning, the region for forming pattern borders should have a short transition width with low sensitivity coefficients to morphogen and receptor synthesis rates.

Deriving a set of explicit formulae for sensitivity coefficients and the transition width, we next study the trade-offs that arise among the three performance objectives: sensitivity coefficient to morphogen synthesis, sensitivity coefficient to receptor synthesis and transition width. We also discuss how these trade-offs influence the maximum range over which patterns can reliably form.

In all following analytical studies, we assume that receptor saturation has a low value. With this assumption and setting the left-hand sides of systems (2.1)–(2.3) to be zero (assume that all the variables are time independent), Sig(*x*) can be approximated by solving
2.14and
2.15where [*L*]_{ss} is the steady-state solution of the morphogen concentration and *K*_{R}, the Michaelis constant of binding, is defined as . The length scale of the signal gradient can be approximated by a constant
2.16Computational simulations, without any prior assumption of receptor saturation, are used to confirm the analytical results.

#### 2.3.1. Upper bounds of length scale and distance for forming robust and precise patterns

A robust patterning system with good precision should have a short transition width and low sensitivity coefficients to both morphogen and receptor synthesis rates. Here, we define a ‘useful patterning region’ as a spatial region satisfying the following three conditions:
2.17where *C*_{W}, *C*_{L} and *C*_{R} are the upper bounds for the three performance objective functions. The upper bounds can be determined by the sizes of the cells and the tissues for a specific biological system. For example, if the diameter of each cell is approximately 1 μm, as in the *Drosophila* wing disc, and we assume that the transition width at any location is less than four cells in the useful patterning region, *C*_{W} should be less than 4 μm.

Because the signal is considered to be a Poisson-distributed random variable, the effect of the noise in the signal becomes important in the region where the signal level is low. In such a region, the level of receptor saturation is likely to be very low, such that the transition width can be approximated as
2.18where *x* > *x*_{1} and *x*_{1} is large enough to satisfy the condition for low receptor saturation. By assuming low receptor saturation everywhere, we can also obtain
2.19Using this approximation, along with two of the three constraints, , we obtain two inequalities involving *x*

Hence, we estimate that the useful patterning region *x* is within
2.20and that the width of this region is approximately
2.21

The useful patterning region exists only if the width of the region is positive; thus, one of the necessary conditions for the existence of the useful patterning region is
2.22According to (2.22), the constant *λ*_{U} has an upper bound that restricts the distance that morphogens can reach from the source. Additionally, from (2.21), the relative width of the region (the width of the region over the upper bound of the region) can be approximated by , which is a decreasing function of *λ*_{U} if *λ*_{U} satisfies the condition (2.22). This shows that the relative width decreases when *λ*_{U} increases for a long-range signal that governs farther patterning from the source.

We directly simulate the original model (2.1)–(2.3) for the case without any assumption on low receptor saturation. Instead of selecting one set of parameters for the model, we study a large number of sets of parameters randomly chosen in the ranges based on experimental observations for Dpp gradients in the *Drosophila* wing imaginal disc (table 3). The simulations are consistent with the analytical estimate, demonstrating that the region for *W*(*x*) < *C*_{W} = 6 μm and exists when the length scale is smaller than a certain level (figure 3*a*,*b*). As the length scale increases, the relative width of the useful patterning region decreases, exhibiting a trade-off between long-range patterning and the relative width of the useful patterning region.

#### 2.3.2. Lower bounds of the length scale and distance for forming robust and precise patterns

The sensitivity coefficient of positional information to the receptor synthesis rate can be approximated by solving the system (2.14) and (2.15) and substituting the solution into (2.11). The solution of the system (2.14) and (2.15) is

By substituting the solution into (2.11), we obtain 2.23

If *C*_{R} is less than 0.5, the condition implies that the useful patterning region is within
2.24

if *C*_{R} is larger than or equal to 0.5, the condition leads to
2.25

In this paper, we assume that *C*_{R} < 0.5, which implies that the pattern border changes less than 50 ln(2)/100 = 34.7% when the receptor synthesis rate is doubled.

Here, the condition (2.24) with another constraint implies that the useful patterning region is within 2.26and the width of the useful patterning region should be 2.27The useful patterning region exists only when the width is positive, which results in a necessary condition for the existence of the useful patterning region 2.28

The left-hand side of the condition (2.28) increases with respect to *λ*_{U}, so the condition (2.28) provides a lower bound for the length scale *λ*_{U}. This condition is usually very difficult to achieve; for example, when *C*_{R} = *C*_{L} = 0.3 (, as used in [17]), the condition (2.28) can be rewritten as
2.29

The left-hand side has an upper bound −(2/3)*x*_{p}; thus, the condition (2.29) fails for any *C*_{W} and any set of parameters.

This result is consistent with computational studies in the system (2.1)–(2.3). By systematically searching in a large parameter space with *C*_{L} = *C*_{R} = 0.3, we do not find a region satisfying both and . Both the analytical estimate (2.28) and the simulations of the system suggest that the useful patterning regions can be found only when *C*_{R} > 0.35. This result implies that if there is no regulation for robustness, the pattern border changes at least 35 ln(2)/100 = 24% when the receptor synthesis rate is doubled. The results in figure 3*c*,*d* show that a region satisfies both and only when the length scale *λ*_{U} is larger than 10 μm.

### 2.4. Mechanisms to relax the constraints

In §2.5, we showed that if there is a patterning region with the sensitivity coefficient to the morphogen synthesis rate less than 30% () and the transition width less than 6 μm (*W*(*x*) < *C*_{W} = 6 μm), the length scale *λ*_{U} and the distance for patterning need to be less than 12 and 60 μm, respectively, from the source (figure 3*a*,*b*). Instead of the transition width, if the constraint for the sensitivity coefficient to the receptor synthesis rate is considered () along with the constraint for the sensitivity coefficient to the morphogen synthesis rate (), the length scale *λ*_{U} and the distance for pattering should be larger than 20 and 80 μm, respectively, from the source (figure 3*c*,*d*). These results indicate that it is impossible to have a patterning region when all three performance objectives, as shown in figure 4*a*, are considered simultaneously; this is also verified with a large number of sets of parameters (figure 4*c*).

Receptor production regulations and the presence of non-signalling binding sites have been suggested to play an important role in the robustness and precision of pattern formation [7,14,19,20,30,32,33,40–42]. Hence, we consider two additional mechanisms: (i) instead of a uniform receptor synthesis rate, lowering the receptor synthesis rate in the morphogen production region [31] and (ii) including non-signalling binding sites (non-receptors) [32–34] in the system (2.1)–(2.3). In particular, we use the analytical estimates of the trade-offs among the three performance objectives to investigate how the size of the useful patterning region may be increased through these two mechanisms.

#### 2.4.1. Lowering the receptor synthesis rate in the morphogen production region improves the lower bound of the length scale and the maximum distance over which robust patterning can be achieved

In *Drosophila* wing imaginal discs, the expression of Tkv (Dpp receptor) in the Dpp production region is lower than that in other areas due to the downregulation by Hh through master of thickveins (mtv) [31]. Such a mechanism can decrease the lower bound (obtained from (2.28)) of the length scale and the distance for a robust patterning region. We demonstrate this result using an analytical approach followed by direct simulations of the full system.

First, the non-uniform receptor synthesis rate is added to the system (2.1)–(2.3) by replacing *v*_{R} with a discontinuous function
2.30where *θ* ≤ 1 represents the ratio of the receptor synthesis rate in the morphogen production region to that in the other region. When the level of receptor saturation is low in the entire domain, the sensitivity coefficient of positional information to the receptor synthesis rate can be approximated as
2.31where
3.32If *C*_{R} is less than 0.5, the condition implies that the useful patterning region is within
2.33The function *F*(*θ*) is decreasing with respect to *θ*, showing that the maximal distance of the region (2.33) is increasing when the level of receptor synthesis reduction is increasing. The improvement in the size of the region (2.33) may be understood as the following. When the receptor synthesis rate increases, more morphogen molecules are captured by receptors, leading to a decrease in the length scale of the morphogen gradient. Hence, the level of morphogen decreases in the region far from the production region. Lowering the receptor synthesis rate in the morphogen production region allows more morphogens to diffuse out from the production region because, near the production region, the level of the signal increases more compared with the case of the uniform receptor synthesis rate. This balance improves the robustness in regions far from the production region.

Additionally, the sensitivity coefficient of positional information to the morphogen synthesis rate can be approximated as
2.34In addition to the constraint , the range (2.33) is refined to
2.35Interestingly, this implies that when , the original condition (2.28) will be replaced by
2.36which is independent of *λ*_{U}. In other words, a strong reduction of receptor synthesis in the morphogen production region can lead to the useful patterning region without any restriction on *λ*_{U}.

Although no restriction on *λ*_{U} resulted from (2.36), the sum of twice the threshold of the sensitivity coefficient to the receptor synthesis rate and the threshold of the sensitivity coefficient to the morphogen synthesis rate must be larger than one, suggesting that the condition is more sensitive to the threshold of the sensitivity coefficient to the receptor synthesis rate. In particular, even with a strong reduction of receptor synthesis in the morphogen production region, the useful patterning region still cannot be formed when *C*_{R} = 0.35 and *C*_{L} = 0.3, which violate the condition (2.36).

Figure 4*d* demonstrates the contribution of lowering the receptor synthesis rate in the source through direct simulations in a large parameter space for the system (2.1)–(2.3) with the modification (2.48). In the simulations, we consider a small *θ*, which is equal to 0.01. Compared with figure 4*c*, which corresponds to the model without any regulation, lowering the receptor synthesis rate in the source increases the number of cases in which the useful patterning region exists and where the relative width of the useful patterning region can be larger than 0.1 at a distance of 40 μm from the source.

#### 2.4.2. Non-signalling binding sites reduce the sensitivity coefficient to the receptor synthesis rate through buffering

Other than signalling receptors, morphogens also bind with non-signalling binding sites (termed ‘non-receptors'). Heparan sulfate proteoglycans, such as dally, syndecans and glypicans, are found to act as non-receptors for Dpp, Wg and Hh systems [32–34]. Modelling studies have suggested that non-receptors could play an important role in enhancing the robustness of pattern formation [19,20,43]. Here, we include non-receptors in our system to study how they may relax the constraints to ease the trade-offs.

With the addition of non-receptors to the system (2.1)–(2.3), the one-dimensional system becomes
2.37
2.38
2.39
2.40
2.41where [*N*] and [*LN*] denote the concentrations of non-receptors and the morphogen–non-receptor complexes, respectively. Other than receptors, non-receptors also bind with morphogens to form morphogen–non-receptor complexes with a binding rate constant *a*_{N} and a dissociation rate constant *β*_{N}. In the system (2.37)–(2.41), non-receptors are produced in the entire domain with a rate *v*_{N}, and non-receptor and morphogen–non-receptors degrade with rate coefficients *d*_{N} and *d*_{LN}, respectively.

With the assumption of low receptor and non-receptor saturations in the entire region, the steady-state equation can be approximated as
2.42where *K*_{N}, the Michaelis constant of morphogen–non-receptor binding, is defined as . The length scale of the signal gradient can be approximated as a constant
2.43The constants and represent the rate constants of degradative flux induced by receptors and non-receptors, respectively. When *ρ* is defined as the ratio of the rate constant of degradative flux by receptors to the sum of the rate constants of degradative flux, i.e.
2.44the sensitivity coefficient to the receptor synthesis rate can be approximated by
2.45The sensitivity coefficient to the receptor synthesis rate in the region far away from the source can be approximated by 1/2*ρ*, which decreases as *ρ* decreases.

When non-receptors reduce the sensitivity coefficient to the receptor synthesis rate through their actions as a buffer, the contribution of non-receptors is accompanied by an increase in the sensitivity coefficient to the non-receptor synthesis rate, which may introduce another constraint for robust patterning. For such cases, the sensitivity coefficient to the non-receptor synthesis rate is approximated by
2.46Equation (2.46) shows that the sensitivity coefficient to the non-receptor synthesis rate increases when the ratio of degradative fluxes *ρ* decreases. An extra constraint may restrict the maximum range of the useful patterning region even though non-receptors relax the condition imposed by the sensitivity coefficient to the receptor synthesis rate.

When we consider the useful patterning region with four performance objectives, , , and *W* < *C*_{W} = 6 μm, as shown in figure 4*b*, with different *ρ*, *ρ* = 0.1, 0.25, 0.75 and 0.9, we observe that the useful patterning region can exist only in the case that *ρ* is 0.75 (figure 4*e*).

The presence of non-receptors reduces the sensitivity coefficient to the receptor synthesis rate through buffering. When the ratio of non-receptor concentration is small, i.e. *ρ* is close to one, the sensitivity coefficient to the receptor synthesis rate is still too large. This explains why the useful patterning region does not exist when *ρ* = 0.9 in figure 4*e*.

However, for non-receptors to work appropriately to alleviate the trade-off, the degradative flux rate constant of non-receptors cannot be much larger than that of the receptors, i.e. *ρ* cannot be too small. Otherwise, most of the morphogens are captured by non-receptors, leading to a decrease in signal, which in turn causes an increase in the transition width. Another constraint, , also increases with small *ρ*, although non-receptors reduce the sensitivity coefficient to the receptor synthesis rate. Thus, a useful patterning region does not exist for small *ρ*, *ρ* = 0.1, 0.25, in figure 4*e*.

#### 2.4.3. A mixed mechanism of non-signalling binding sites (non-receptors) and lowering the receptor synthesis rate in the morphogen production region broadens the robust and precise patterning region

While including non-receptors or lowering the receptor synthesis rate in the source relaxes certain constraints, as shown above, a direct computation of the system shows that inclusion of one mechanism alone still does not yield more patterning regions. For example, when *C*_{L} = 0.3, *C*_{R} = 0.4, *C*_{R} = 0.3 and *C*_{W} = 6 μm, the distance of the useful patterning region is still less than 50 μm from the source (figure 4*d*); only a few cases can achieve the conditions for robust patterning (figure 4*e*). Here, we apply direct simulations (figure 4*f*) to show that a mixed mechanism of non-receptors and lowering receptor synthesis in the morphogen production region can increase the number of cases and expand the range of the useful patterning region.

In *Drosophila* wing imaginal discs, the level of dally (non-receptor) expression is high in regions where Dpp is synthesized, but by contrast, tkv (receptor) expression levels are low in that region; dally expression is downregulated in other regions where tkv expression levels are high [33], as shown in figure 4*f*. Now, we model this expression pattern by modifying the system (2.37)–(2.41) considering non-uniform receptor and non-receptor synthesis rates:
2.47and
2.48where *θ*_{1} ≤ 1 represents the ratio of the receptor synthesis rate in the morphogen production region to that in the other region; *θ*_{2} ≥ 1 represents the ratio of the non-receptor synthesis rate in the morphogen production region to that in the other region.

Figure 4*f* displays two sets of simulations with different combinations of *ρ*, and both cases show that a mixed mechanism of non-receptors and lowering receptor synthesis in the morphogen production region expands the useful patterning region with four performance objectives farther than 50 μm from the source. Although there is one more constraint for non-receptors, the non-uniform receptor and non-receptor synthesis rates can overcome the trade-offs to form farther and larger regions for robust patterning.

## 3. Discussion

Morphogens are a key factor for growth and patterning during development. Morphogens interact with many extracellular components and closely tie the details of patterning to the amount of receptor occupancy at different spatial locations. Thus, the ability of development to proceed normally in the face of genetic and non-genetic perturbation demands robustness of the morphogen patterning system.

As shown here, a morphogen system might not require a set of complicated mechanisms to be reasonably robust with respect to one particular genetic perturbation, for example, perturbation in morphogen synthesis. However, the system has different challenges in the face of multiple genetic perturbations, for example, fluctuations in both receptor and morphogen syntheses. With additional performance objectives, such as the transition width that measures the effect of local stochastic fluctuations, a trade-off due to various constraints usually arises.

We have delineated such trade-offs by examining constraints on the length scales of the signal gradient to achieve the prescribed robustness. Through a general analytical formula derived without specifying the values of the system parameters, we have demonstrated that a mixed mechanism combining both lowering the receptor synthesis rate in the morphogen production region and the presence of non-signalling binding sites enlarges the size of the useful patterning region.

The results described here can be applied to the Dpp gradient in *Drosophila* larval wing discs, as substantial evidence points to a role for non-signalling binding sites (heparan sulfate proteoglycans) in controlling the shape of the Dpp gradient, and the reduction of receptor Tkv synthesis in the Dpp synthesis region is well established. It is also known that in *Xenopus* embryos, removing the heparan sulfate binding domain from BMP4 (an orthologue of Dpp) greatly increases its range of action as a morphogen [44]. This work potentially provides a useful approach for analysing other regulation mechanisms, such as self-enhanced clearance regulation observed in Wg gradients [14,17] and negative feedback from the signal to the receptor in the Dpp gradient [32] in the *Drosophila* wing disc, in the context of patterning robustness. One interesting question would be how the two opposite regulation mechanisms affect potential trade-offs and the constraints observed in this work.

Achieving multiple performance objectives and overcoming constraints in a morphogen system require multiple strategies, complex regulations and other temporal and spatial information in morphogen systems. One natural strategy for improving the robustness of a morphogen gradient is using the slope information of the signal gradient, which has been observed to play an important role in regulating cell proliferation [45]. The temporal information of the morphogen gradient on patterning could also contribute to mediating trade-offs and to relaxing constraints [46,47]. Pre-steady-state decoding of morphogens [48,49] and intracellular signalling cascades, which result in the activation of transcriptional effectors [35], may also be critically important in the robustness of developmental patterning. Cell-to-cell variability [24] and embryo-to-embryo variability [50] will likely lead to additional complexity for attaining robustness and dealing with stochastic dynamics in morphogen systems. Interestingly, noise in downstream gene regulations has been found to enable sharpening of borders between different gene expression domains specified by morphogens [51]. How do these factors affect the system's capability in achieving multiple performance objectives simultaneously? The presented theoretical and computational approaches will be useful in addressing this challenging question.

## 4. Material and methods

Newton's method is used to calculate the solution of steady-state systems in which the diffusion is approximated by a central difference scheme with ▵*x* = 0.1 For Newton's method, the tolerance for convergence is 10^{−6}. Finer resolutions with smaller ▵*x* and smaller tolerance are used to check the convergence of the steady-state solution. The ranges of parameters used in the simulations are displayed in tables 2 and 3. It is worth noting that we use wide ranges of the binding rates *α*_{R} and *α*_{N} containing all possible values observed in the experiments. Therefore, some values in the ranges may be very close to the limits for the diffusion for bimolecular interactions. To study how the range of parameters may impact the solution, we also consider another range for binding rates, 10^{−}^{3}*–*10 μM^{−}^{1} s^{−1}. We find no change in the qualitative results observed in figure 4, and the range of the useful patterning region remains the same.

All the cases considered here must satisfy , which guarantees that the convexity of signal gradients is consistent with experimental observations. Over 80% of cases satisfy this condition in the parameters searched.

The calculations are performed using FORTRAN 77, and plots and data analyses are produced using matlab7 (Mathworks, Natick, MA, USA).

## Funding statement

This work was partially supported by NIH grant nos. R01GM67247, P50GM76516 and R01GM107264 and NSF grant nos. DMS0917492 and DMS1161621. W.L. has been supported in part by the Mathematical Biosciences Institute through NSF grant DMS 0931642.

- Received September 17, 2014.
- Accepted October 28, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.