ROBUSTNESS OF SIGNALING GRADIENT IN DROSOPHILA WING IMAGINAL DISC.

Quasi-stable gradients of signaling protein molecules (known as morphogens or ligands) bound to cell receptors are known to be responsible for differential cell signaling and gene expressions. From these follow different stable cell fates and visually patterned tissues in biological development. Recent studies have shown that the relevant basic biological processes yield gradients that are sensitive to small changes in system characteristics (such as expression level of morphogens or receptors) or environmental conditions (such as temperature changes). Additional biological activities must play an important role in the high level of robustness observed in embryonic patterning for example. It is natural to attribute observed robustness to various type of feedback control mechanisms. However, our own simulation studies have shown that feedback control is neither necessary nor sufficient for robustness of the morphogen decapentaplegic (Dpp) gradient in wing imaginal disc of Drosophilas. Furthermore, robustness can be achieved by substantial binding of the signaling morphogen Dpp with nonsignaling cell surface bound molecules (such as heparan sulfate proteoglygans) and degrading the resulting complexes at a sufficiently rapid rate. The present work provides a theoretical basis for the results of our numerical simulation studies.

morphogens (aka ligands) are synthesized at a localized site. These morphogens disperse from their production site; some bind to cell receptors along the way, generally resulting in different receptor occupancies at different cell locations. The spatial concentration gradient of morphogen-receptor complexes (aka bound morphogens) induces spatially graded differences in cell signaling. The differential cell signaling in turn gives rise to different gene expressions from which follow different stable cell fates and visually patterned arrangements of tissues and organs during development.
In principle, the process of forming morphogen gradients leading ultimately to tissue patterning consists of syntheses of transportable morphogens and membrane bound cell receptors, their binding and dissociation, endocytosis and exocytosis of morphogenreceptors and their intracellular degradation. This collection of biological processes that explicitly include endocytosis and exocytosis has been modeled mathematically as System C in [18] from which we have deduced by analysis and computation how the shape of the signaling gradient depends on the system parameters such as synthesis rates of morphogens and receptors, binding and degradation rate constants, etc. We also see from the mathematical model that small changes of these system parameters may cause substantial changes in gradient shape [18]. In contrast, embryonic patterning is usually highly robust, resisting not only substantial changes in the expression level of individual genes, but also fluctuating environmental conditions (e.g., unseasonal heat waves). This suggests that additional biological processes must also be at work to ensure such robustness. Identifying the cause of robustness and ways of producing robust morphogen gradients have become a major research effort in recent years [6,7,8,13,14,15,23,24,31,32].
A reasonable supposition would be that robustness is the result of various types of feedback control mechanisms. For example, the amount of signal received by a cell may influence the amount of receptors it makes for the particular type of signaling morphogen. Another feedback mechanism may be the up-regulation of receptor-mediated degradation rate by cell signaling. The Drosophila wing imaginal disc is patterned by the gradient of the decapentaplegic (Dpp) morphogen, a member of the bone morphogenetic protein branch of the transforming growth factor-superfamily. Dpp signaling represses synthesis of its receptors, but enhances Dpp degradation [5,22]. Analytical and numerical simulations of a model system that includes feedback [15,20] showed that repression of receptor synthesis rate alone (without enhancing morphogen degradation) does not sustain the expected robustness. This is supported by the results in [7] showing that additional biological activities were needed to attain robustness for the morphogen gradient and prompted considerations for alternative paths to morphogen gradient robustness in [20]. A major finding by numerical simulations of various model problems in [20] is that feedback control is neither necessary nor sufficient for robustness for the Dpp gradient in wing imaginal disc of Drosophilas. In addition, the numerical results suggest that robustness of the Dpp gradient can be achieved by substantial binding of the signaling morphogen Dpp with nonsignaling cell surface bound molecules, such as cell surface heparan sulfate proteoglycans, called nonreceptor for brevity, and degrading them at sufficiently rapid rate. (The former is to be the consequence of high occupation of receptors and low occupation of non-receptors while the latter means a high degradative flux of non-receptor relative to that of the receptors.) No feedback is required throughout this slightly more complex process of gradient formation. That non-receptors can be solely responsible for robustness may explain, in part, existing and growing evidence that nonsignaling molecules are usually present and participating actively in morphogen gradient formation (through binding with the signaling morphogen) [25].
The robustness of a morphogen gradient is relevant only if the morphogen gradient is biologically capable of inducing differential cell signaling (multi-fate signaling). A multi-fate signaling should broadly distribute pattering information over the entire field of cells so that multiple types of cells can be developed. A quantitative measure of biologically realistic multi-fate signaling morphogen-receptor gradients was first introduced in [18] in terms of the magnitude, steepness and convexity of the gradient. The measure is further quantified numerically in [20,21] to give numerical yardsticks to geometrical features of "too steep/too narrow" "too wide" and "too convex". In [20], numerical simulations were carried out for 2 20 (or more than more than 10 6 ) random sets of parameter values for each of the model systems. After discarding the parameter sets that result in biologically unrealistic gradients that would result in most cells in the wing disc developing into the same cell type, the robustness of the remaining (biologically multi-fate signaling) cases were examined with respect to a range of discrete values of one of the five important normalized parameters for a very large number of random sets of the other parameters. These plots enabled us to see the existence of robust multi-fate gradients with the addition of non-receptors alone (without feedback) to System C and the nonexistence of robust multi-fate gradients with down regulating feedback alone (without non-receptors).
The findings in [20] positioned us to develop a theoretical foundation for the conclusions suggested by the numerical simulations. Specifically, we develop an existence proof of robust multi-fate gradients in a non-empty region of the parameter space of the biological system, the Dpp-receptor gradient in the imaginal disc of Drosophilas, with respect to a substantial (two-fold) change of ligand synthesis rate. We will use the same criterion for robustness introduced in [20,21] but will work with a more general set of criteria for a multi-fate gradient. Together, they will enable us to analytically locate a region of multi-fate gradients in the parameter space which are robust with respect to the ligand synthesis rate. As such, we will have a considerably more complete and explicit understanding of the dependence of multi-fate and robustness on the system parameters. It is expected that the same analytical method will also enable us to extract useful information on robustness with respect to other system parameters.

The Initial-Boundary Value Problem
The Dpp-receptor gradient system in the wing imaginal disc of Drosophilas analyzed in [18] involves concentrations of free ligands (Dpp) [L], cell membrane bound receptors (Tkv) in extracellular space [R] out and cell interior [R] in , and the extra-and intracellular morphogenreceptor complexes [LR] out and [LR] in . The morphogen-receptor complexes inside cells [LR] in provide the signal to activate the target genes potomotorblind ( omb) and spalt (sal); it is the gradient of [LR] in that determines the fate of cells. Synthesized locally over a few cells between the posterior and anterior compartment at a rate V uniformly in the directions orthogonal to the anterior-posterior axis, the ligands diffuse away from the source and bind to signaling cell receptors along the way with the ligand-receptor complexes transcytosing and degrading in the cell interior. Receptors are synthesized at a constant rate R to replenish losses through degradation of both [LR] in and [R] in . The evolution and interaction of ligands and receptors will be as characterized by the space-time model System C of [18]. Given the nature of the morphogen source, morphogen activities essentially vary only in the direction of anterior-posterior axis X, from the ligand source to the edge of the wing disc.
For the present robustness study, we add to System C concentrations of non-receptors [N] out and [N] in , morphogen-non-receptor complexes [LN] out and [LN] in ; they transcytose and degrade in a way similar to the morphogen-receptor complexes. For this extension of System C, we have the following spatially one-dimensional system of differential equations first introduced in [20]: (1) (2) (8) (9) for −d 0 < X < X max with V (X) given in terms of the Heaviside unit step function H(z): (10) where v 0 is a constant so that morphogens are synthesized uniformly in the region −d 0 ≤ X < 0 only. The synthesis rates of receptors and non-receptors are taken to be uniform in time but R and N may be piecewise constant in X. No feedback is considered in the present model.
We consider here the development of only the posterior compartment of the wing disc including half of the morphogen production region. At the border between the two compartments, X = −d 0 , we have the no flux conditions (11)  as a consequence of the condition of symmetry relative to the border. At the other end, there are very few free ligand molecules not bound to a receptor or non-receptor; hence we may treat the edge as a sink so that: (12) With V (X) discontinuous at X = 0, we stipulate also the continuity of [L] and ∂[L]/∂X at X = 0.
Before the onset of morphogen production at T = 0, we have no morphogen concentration of any kind so that (13) The receptors and non-receptors are expected to be in steady-state prior to the onset of ligand production so that (14) These conditions lead to the steady state values that constitute the remaining initial conditions: (15) for the concentrations of receptor and non-receptor.
The system above, designated as System CN, is formally a straightforward extension of System C and reduces to the latter in the absence of non-receptors. However, they differ substantively from System C in that we have included as in [19] the region of morphogen synthesis as a part of the solution domain. In acknowledging the presence of a region of ligand synthesis and the need to consider the molecular dynamics in that region in conjunction with the other ligand activities, we have made the problem more complicated and must deal with its consequences, including allowing the morphogen-production cells to have receptors and non-receptors that bind with some of the morphogens they produced [10,11].

Steady-State Behavior
After the onset of morphogen production, concentration gradients of the different concentrations form rapidly reaching some quasi-steady state configuration in a matter or hours or less. It is the robustness of the steady state gradients that is of current interest in development. Upon setting all time partial derivatives to zero, all but the first of the set of governing partial differential equations become algebraic equations that can be solved to obtain (16) where the intrinsic parameters are Upon setting (18) we obtain from 1, with ∂[L]/∂T = 0, the following dimensionless ODE for the normalized steady state free ligand concentration a(x): (22) with (23) where = (v 0 /R o )/k degobs is the ratio of the (normalized) ligand production rate to the (normalized) observed degradation rate of the ligand-receptor complexes first introduced in [18] to help characterize the steady state level of ligand concentration.
For the second order ODE with a discontinuous forcing term, we have the boundary conditions (24) and the continuity conditions of a and a' at x = d. The normalized concentration of the signaling morphogen-receptor complexes is given by (25) Hereafter, we will use uppercase letters (X, [L], [LR], etc.) to denote the original variables, and lowercases (x, a(x), b(x)) to denote the normalized/dimensionless variables; they are related by 18-21. Note that this convention does not apply to the system parameters.
Biologically, free and bound morphogens form gradients outside the production region rapidly. At steady state, the gradient of the signaling ligand-receptor concentration should be capable of inducing different cell fates at different cell locations. Moreover, the signaling gradient and the resulting tissue pattern should be highly robust notwithstanding substantial system parameter changes (e.g. a two-fold change in the expression of individual genes) resulting from fluctuation of environmental conditions (e.g. unseasonably high or low temperature). In this paper, we will be concerned with the role of non-receptors in the robustness of signaling gradients. In the next section, what constitutes an admissible signaling gradient for multi-fate development (or multi-fate gradient for brevity) will be defined quantitatively and an analytical measure of robustness will be introduced for signaling gradients. Together, they will provide us with the quantitative criteria for analyzing how non-receptor affects the robustness of multi-fate gradients.
In the boundary value problem (BVP) for the steady state free ligands defined by 22-24, there are five parameters d, , , p, and v. The parameter d is the relative width of the production region of the morphogen. In this paper, we will always take d to be a prescribed quantity (with d = 0.06 corresponding to the width of 12 m of the production region compare with 200 m, the width of the Drosophila wing imaginal disc). The parameter v is the only one involving the rate of morphogen infusion v 0 and may be taken as normalized ligand synthesis rate. The parameter = R k m / N j m , is seen to be the ratio of the saturation level of receptors to that of non-receptors found to be important factors for robustness in [20]. Similarly, the ratio is seen to be of the order of the ratio of degradative fluxes of the receptor, R = k degobs [LR] out , to that of non-receptor, N = j degobs [LN] out , previously introduced in [20] when the steady state ligand concentration is relatively low. (We may also consider the ratio to characterize the relative magnitude of the (normalized) synthesis rate of ligand receptor and that of nonreceptor since p/(1 − p) = O( R / N ).) In the degradative flux interpretation, is seen to be of the order of the sum of these fluxes. The numerical results in [20] suggested that for System CN, only a certain combination of these flux and saturation factors would induce robust signaling gradients capable of differential cell signaling. We will provide an analytical validation for the observations in [20] and quantify more precisely the conditions for their validity. For these and other results on robustness of signaling gradients, we will need some specific properties of the solution of the steady state BVP 22-24. These will be developed in next few subsections.

Monotonicity of the Unique Steady State Solution
For a given set of the non-negative values of the parameters d, , , p, and v, the following existence and uniqueness theorem for the BVP 22-24 can be proved by the same method as that used in [19] (see also Theorem A.1 of the Appendix of this paper): Theorem 2.1-A unique non-negative solution a(x) exists for the BVP 22-24 with 0 ≤ a(x) ≤ a u (x) where the upper solution a u (x) for the problem is given by (26) To study the properties of the steady state solution, we introduce the abbreviation It is easy to see that (27) (28) for all a > 0. The following monotonicity properties of a(x) are less straightforward: Proposition 1-Let a(x; , v) be the unique steady state solution of 22-24 (which of course depends also on the parameters p and as well). Then for all positive and v, we have (29) Proof: Upon differentiating 22-24 with respect to v and setting (x; , v) = ∂a/∂v, we see that is determined by the BVP (30) where H(z) is the Heaviside unit step function and (31) Apply Theorem A.1 to the BVP 30, we have (x) = ∂a/∂v > 0 for all x (0, 1).

Exact Steady State Solution
The second order ODE 22 is autonomous and can be integrated to give x as a function of a: The three unknown constants of integration a 0 a(0), a d a(d) and s 1 a'(1) are determined by a(1) = 0 and the continuity of a and a' at x = d with the last of these three conditions requiring (37) The derivation of this exact solution is similar to that for the case of no non-receptors in [19] and will not be given here.

Low Ligand Synthesis Rate (LLSR)
Though equations 33-35 give the exact solution of the BVP, insight to steady state behavior is not readily accessible from these expressions. We obtain in this subsection an explicit solution which is a leading term perturbation solution (and an adequate approximation) of the exact solution for low normalized morphogen synthesis rates (corresponding to low occupation for both receptors and non-receptors in [20]). As we shall see, it also provides a useful tool to decipher the implications of the exact solution for intermediate range of morphogen concentration.
For a sufficiently small normalized synthesis rate, we expect 0 ≤ a(x) 1 and 0 ≤ a(x) 1 (corresponding to [L] R k m and [L] N j m ). In that case, a leading term approximate solution of the steady state problem is determined by (38) Similar to the method in [19], we have for this low ligand synthesis rate (LLRS) case a L (x) The monotonicity properties of a(x) of Proposition 1 apply to K(x; ). Some of these can be seen directly from the explicit solution above. For example, we have K'(x; ) < 0 from (40) Consistent with the leading term LLSR approximation, we have (41) Note that in the LLSR range, b(x) b L (x) = vK(x; ) depends only on v and (and of course the synthesis region width d which is assumed to be fixed in this paper) and not on p and . Furthermore, the dependence on v is linear so that the magnitude of b(x) b L (x) varies in proportional to the ligand synthesis rate. It follows that development would be sensitive to a variation in v (caused by environmental changes for example) and therefore in some sense not "robust". On the other hand, the sensitivity with respect to the ligand synthesis rate would vary with the convexity of the gradient and hence with the value of the parameter . The significance of the actual variation from the combined effect would depend on how we quantify robustness. We will address the issue of an appropriate measure of robustness in the next section. For that purpose, the following properties of K(x; ) will be useful: The last two properties of k( , d) imply k( , d) < 0 for any d in (0, 1), and the third property is proved.

High Ligand Synthesis Rate (HLSR)
At the other end of the spectrum where the ligand synthesis rate is sufficiently high so that vd is large compared to max {1, 2 , 2 / }, we have a case of high occupation of receptors and non-receptors discussed in [20]. In that case, the leading term approximation a H (x) for the steady state solution is determined by the BVP (42) or (43) Correspondingly, we have (44) But unlike the LLSR case, the approximation of b(x) by b H (x) 1 is valid only for x away from a narrow interval adjacent to the x = 1 end. With a(1) = 0, a(x) is not large compared to unity near x = 1. Except in the boundary layer adjacent to the wing disc edge, b(x) is seen not to depend on any of the four parameters , , v and p to a first approximation for the HLSR range. As such, development is essentially insensitive (and therefore robust with respect) to system and environmental changes that may affect the system parameters. (A formal validation of this observation will be given after we have formulated a quantitative measure of robustness.) On the other hand, the concentration of morphogen-receptor complexes responsible for signaling and development is effectively uniform in nearly the entire solution domain and would not give rise to patterning. In other words, such a ligand-receptor gradient, though insensitive to changes, is not a multi-fate signaling gradient and would not be of interest to the study of real biological systems.

Normalized Root-Mean-Square Displacement
We saw from the explicit solution for the LLSR range that the signaling ligand-receptor gradient is generally sensitive to system parameter changes. Yet actual biological systems are generally robust to such changes. It is our goal to investigate factors that are responsible for such robustness. We do so by focusing on robustness with respect to a two-fold change in the ligand synthesis rate in our model problem as in [20]. The general methodology developed for this parameter change should be helpful to our study of robustness with respect to other parameter changes.
Let b(x) and be the normalized signaling ligand-receptor gradients for synthesis rate v and 2v, respectively and x 1 and x 2 the corresponding location where they attain the value , i.e., . With the change of ligand synthesis rate, x 2 is generally different from x 1 with x 2 − x 1 = x. The root-mean-square of x over the range of b(x) would be a meaningful measure of robustness. To minimize the effects of outliers, the range of b will be taken to be the interval (b 1/5 , b 4/5 ) with b 1/5 = b(d)/5 and b 4/5 = 4b(d)/5. The measure of robustness for our analysis, R v , is this root-mean-square deviation normalized by the interval x(b 1/5 ) − d: (45) In general, the displacement x depends on the normalized signaling ligand-receptor gradients for the two different ligand synthesis rates v and 2v. Since these gradients themselves depend also on the parameters p, and , we indicate these dependence by writing R v (p, , ). It is seen from 45 that the smaller the value of R v the more robust the system would be. As suggested in [20], the system is considered to be acceptably robust with respect to the ligand synthesis rate v if, for a two-fold increase in v, we have R < 0.2.
Numerical solutions obtained in [20] for the steady state behavior of System CN suggest that the corresponding system without non-receptors (System C) does not have any robust multifate gradients for any combination of parameter values. We will validate this observation in the next section after we quantify multi-fate gradients. Before doing this, we will show that with non-receptors, 1) the signaling gradient b(x) is generally robust for sufficiently high ligand synthesis rates but the gradient itself is not a multi-fate gradient, and 2) at low ligand synthesis rates, the signaling gradient b(x) is generally not robust with a value for R v (p, , ) significantly above 0.2 (tending to 0.43 in the limit).

The HLSR Case
As indicated above, we will focus herein on the robustness corresponding to a two-fold increase of the normalized synthesis rate of free ligands, i.e., when v is changed to 2v. For a sufficiently high ligand synthesis rate so that vd ⪢ max{1, 2 , 2 / }, we have from 43 and therewith (46) in the region where ligand is not produced, d < x < 1. When the synthesis rate is changed from v to 2v, the displacement x of gradient at concentration b is given by It follows that the relevant robustness measure R v (p, , ) is given by (47) where b d = b(d). It is easy to see from this expression that the robustness R v can be made smaller than any given positive number if v is large enough. The result summarized in the following proposition provides a mathematical justification of the intuitive expectation in subsection 3.3 using R v as a measure of robustness.: Proposition 3-For vd ⪢ max{1, 2 , 2 / }, the steady state behavior of the model biological system 1-15 is robust.
However, as indicated in subsection 3.3, the signaling gradient [LR] in for such a high ligand synthesis rate cannot form a "realistic" biological gradient for patterning since it is nearly uniform for the entire solution domain except in a narrow region adjacent to the edge x = 1 and at the same time too steep in that narrow layer. While these observations may be evident from a graph of 46, we need to have some quantitative measure for what constitutes a multifate signaling gradient before we can direct our attention to study factors responsible for robustness of such signaling gradients. We quantify in the next subsection what constitutes a multi-fate gradient and use the criteria developed and the robustness measure R v to investigate robustness of multi-fate gradients for the three special cases of high ligand synthesis rates, low ligand synthesis rates and systems without non-receptor. The role of non-receptors in promoting robustness for multi-fate signaling gradients with moderate ligand synthesis rates then be delineated.

Multi-fate Signaling Gradients
In order to induce spatially differential cell signaling, i.e., a multi-fate signaling gradient, a steady state b(x) should have the following characteristics: 1. The slope of the normalized signaling gradient b(x) = [LR] in / R R i that activates the target gene should not be too steep. From 25, we have b(0) < 1. Therefore, the average slope of b(x) over the interval (0, 1) is less than 1. The gradient is considered not too steep if the magnitude of the relative slope |b'(x)/b(x)| in the region of interesting is less than some threshold, i.e., for some > 0.

2.
The concentration of a patterning signal [LR] in should be higher than a certain threshold in the vicinity of the ligand production region. Before the onset of morphogen, the concentration of receptors inside the cell ([R] in ) equal to R i . Thus, the threshold can be assumed to be a fraction of R i . From 25, we see that [LR] in is generally less than R R i and approaches R R i from below when [L] is large enough so that we are near receptor saturation. Hence, for differential cell signaling, the concentration [LR] in threshold can only be a fraction of R R i . We let that threshold fraction be (< 1) so that mathematically the signal is activated if b(d) ≥ with 0 < < 1.

3.
The slope of a(x) at x = 1 should be substantially less than unity. Experimental results had shown that the Dpp form shallow gradient in the imaginal disc [9,29]. These results suggest that the free ligand decays quickly in the imaginal disc in steady state. Motivated by the corresponding relation for the LLSR case, we expect that the slope a'(1) should be considerably less steep than the average slope of free The observation above suggest that we quantify the characteristics of a multi-fate gradient by the following definition: For the purpose of obtaining specific results, we will take = 0.05, = 0.1(or a = 0.11), and = 0.002. While these choices of parameter values may seem somewhat arbitrary, we will see that the results are not particularly sensitive to our choices. (see also 47) so that the third condition for a multi-fate gradient above is not met. Proposition 4 is a negative result. To gain insight to robust development, we need to quantify the ranges of the four parameters p, , , and v for which the corresponding gradients are capable of differential cell signaling, i.e., for which they are multi-fate, by the requirements of Definition 4.1. The following theorem provides the quantification sought:

The LLSR Case
In this section, we examine how the quantitative requirements of a multi-fate gradient applies to the LLSR case and what the resulting expression R v (p, , ) tells us about the nature of steady state signaling gradient robustness when the normalized ligand synthesis rate is increased two-fold. For this case, we have from 41 that b L (x) = a L (x) = vK(x; ) with K(x; ) determined by the BVP 38. While K(x; ) is given explicitly by 39, it turns out to be more effective to work with x(a) as in Subsection 3.1 for the purpose of calculating x and R v (p, , ). By the method of that subsection, it is straightforward to obtain where (51) emphasizes the dependence of s 1 on the normalized ligand synthesis rate v (and of course on as well See [21] for the possibility of size-normalized robustness for the LLSR range. However, Proposition 4 tells us that such a signaling gradient would not be a multi-fate gradient; it would not be capable of inducing differential cell signaling for patterning. Together, they limit the ligand synthesis rate to a moderate range of v values. In the moderate v range however, the BVP for the steady state solution does not admit simplifications that would lead to an explicit solution or useful tool such as a steady state proportional to the ligand synthesis rate. Nevertheless, certain simplifications are still possible in the robustness calculations. In this section, we deduce some of these simplifications and use them to analyze the level of robustness possible for a system without non-receptors, i.e., for System C (instead of CN). Bounds on a(x; λ, v) and a'(x; λ, v) To simplify the expression for R v , we need to establish first some upper and lower bounds on the steady state free ligand concentration a(x; , v) and its derivative a'(x; , v). Let with the last inequality follows from the 54 and the fact that . The monotone increasing positive sequence bounded above by has a limit *; it is the solution of the equation (55) Note that 55 has only one solution since the right hand side is a decreasing function of *. Altogether, we have the following lemma:

5.1.
Lemma 5.1-The monotone increasing positive sequence defined by 54 is bounded above and therefore has a positive limit * ( , v) < which is the unique solution of 55.
The limit * ( , v) enables us to deduce an upper and a lower bound for both a(x; , v) ≥ 0 and a'(1; , v) ≥ 0. The relation 57 is deduced from the following two inequalities: and R v (p, γ, λ) for a multi-fate gradient Lemma 5.2 above and the requirements of a multi-fate gradient place a restriction on the magnitude of : if the second condition is met. It follows from the hypothesis on , the two conditions 60-61, and Lemma 3.1, Therefore, we must have by Lemma 3.1 and therewith .

Simplification of
We are mainly interested in the application of the lemma to multi-fate gradients for which 1. For example, we have for = 0.002.

R v (p, γ, λ) for System C
Let where minimization is taken over all acceptable pairs of { , v} that ensure a multi-fate signaling gradient. Note that if R(p, ) > 0.2, then is always larger than 0.2 for any admissible pair ( , v) that ensures a multi-fate gradient. In that case, it would not be possible to find a parameter pair ( , v) such that the steady-state is both multi-fate and robust. On the other hand, if R(p, ) < 0.2, there exist ( , v) parameter pairs such that the steady state signaling gradient is both multi-fate and robust. Consequently, the quantity R(p, ) provides a more succinct measure of robustness and is used subsequently whenever appropriate.
In the absence of non-receptors so that p = 1, numerical simulations carried out in [20] suggested that all multi-fate gradients are not robust with respect to a doubling of the ligand synthesis rate (see also [21] for a different kind of robustness for low synthesis rates). The theoretical lower bound of the robustness measure R(p, ) in the absence of non-receptors is given below to validate this observation:

Proof:
We have from Proposition 5 a d (2 v) > 2 a d (v). Then for any pair of ( , v), and hence R(p, ) > J(p, ). When p = 1, we have and J(1, ) can be determined numerically to be (0.354527… or) greater than 0.35.
As a direct consequence of Proposition 6, the robustness measure of System C can not be lower than 0.35 for any acceptable pair of ( , v). In other word, without nonreceptor, any multi-fate gradient of System CN (which, without non-receptor, is reduced to System C) cannot be robust. Quantitatively, the robustness of System C (or System CN without non-receptors) has a lower bound of 0.35 for for all parameter sets with a multi-fate signaling gradient. The simulation results of System CN show that this lower bound for can be lowered considerably to well below 0.2 with the addition of non-receptor for certain parameter sets. Experimental results also show that nonreceptor is essential in forming robust morphogen gradients of Dpp in the wing disc of Drosophila (see [1,2,3,4,10,12,14,16,17,28,30]). The theoretical results for System CN of the last few subsections mean that robustness can only be attained for relatively moderate values of the normalized ligand synthesis rate. They help to limit our search in the next section for a region (or regions) in the parameter space where robust multi-fate gradients can be found.

The Role of Non-receptors
From the results of the last section, we know that robust multi-fate gradients are not possible when there is no non-receptors in System CN (leaving us with just System C). From the earlier section, we also learned that multi-fate gradients are also not possible for low or high occupation of both receptors and non-receptors. If non-receptors should be responsible for robustness, the results of numerical simulations in [20] suggest that it would be occur at a level of high receptor occupancy (by ligand) and low non-receptor occupancy. We prove in the first subsection that low non-receptor occupation is necessary for robustness while sufficiency requires some additional consideration as we show in the next subsection.

Non-robustness in Parameter Space
Given by Proposition 6, signaling gradients cannot be robust for pairs (p, ) for which J(p, ) > 0.2 (for all acceptable ( , v)). The graph in Figure 1 shows that J(p, ) is an increasing function of . Thus, the (p, ) plane can be divided by the curve J(p, ) = 0.2 into two regions. Numerical computation shows that the curve J(p, ) = 0.2 can be approximated by (68) In other word, is always larger than 0.2 whenever or . To bring out the role of non-receptors in robustness more explicitly, we introduce a ratio of receptorto-non-receptor (effective) synthesis rate (69) and express the condition for robustness in terms of and (instead of and p up to now). In terms of and , the condition p > p* ( ) becomes > * ( ) with (70) and we have the following sufficiency result for non-robustness:
The two sufficient conditions for non-robustness, and > * ( ) (when ) may be rephrased as the following necessary condition for robustness : However, a multi-fate gradient is required when we write robustness in the form of 67 specified in Definition 5.4. As a consequence, the non-robustness range (with the parameter ) in Proposition 8 is insensitive to the choice of , and .

Region of Robust Multi-fate Gradients
As a direct consequence of Proposition 8, multi-fate gradients of System CN can be robust only if the two nonnegative parameters and are both small enough with 0 < < *( ).
With being a measure of the relative infusion of receptor to non-receptors, this means that non-receptors should play dominant role in forming the multi-fate signaling morphogen gradient. We show below that this condition is also sufficient for robustness. If we use another value R c instead of 0.2 as the upper bound for robustness, then 71 should be replaced by (73) possibly with a different function C b (p, , ) and the conclusion of Theorem 6.1 still holds.
When p and are small enough, we can, by Theorem 6.1, always find parameters ( , v) in G p, such that the system has robust signaling gradients. The region can be found numerically from Theorem 9 in the Appendix. The results for sample points on the boundary of are given at Table 1 (see also Figure 3(a)). From Theorem 6.1, for any pair , there exists a function C a (p, , ) (depends on the parameters d, , , as well) such that when ( , v) G p, and is bounded above and below by two curves: (74) then the pair ( , v) is acceptable to (p, ). For parameters studied in this paper, numerical computation shows that the curve v = C a (p, , ) is identical to the upper bound of the region G p, (See Figure 2).
From Table 1, the domain can be given approximately by (see also Figure 3(a)) (75) Furthermore, from the numerical results (not included herein), the function C b (p, , ) can be approximated by (76) where {C bi (p, )} are, respectively, ratios of second and first degree polynomials of p and . In particular, we have C b (p, , ) > 0. The function C b (p, , ) is found by minimizing the square error 2 = i |r i | 2 where r i are the difference between each original data point and its fitted value. In our sample study, we have used over 4300 data points giving only a 7% square error for C b (see Figure 4). The upper bound C a (p, , ) in 74 is consistent with the upper bound of G p, with analysis formulate given in Appendix A.3.2 (See Proposition 9).
In terms of , the relations 74-76 can be rewritten as (see also Figure 3 The two conditions (i) and (ii) in the part (2) of Theorem 6.2 may be given in terms of instead of p in which case 74, 75 and 76 would be replaced by 77, 78 and 79, respectively.

Concluding Remarks
In this paper, we examine the robustness of steady state morphogen gradients capable of differential cell signaling with respect to a two-fold change of morphogen production rate. Quantitative measures of multi-fate signaling gradients and robust of signaling gradients are specified and used to delineate the occurrence of robust multi-fate gradients in the parameter space. By mathematical analysis, we succeeded in validating the simulation results in [20]. The  not too high to saturate the available receptors so that the signaling ligand-receptor gradient remains multi-fate. As long as there are unoccupied receptors, high ligand synthesis rate would continue to produce more ligand to saturate them unless these additional ligands can be otherwise engaged and (proportionally) unavailable for binding with the unoccupied receptors. The presence of abundance of non-receptor with strong affinity for binding with ligand and for rapid degradation of the resulting non-signaling ligand-non-receptor compounds provides the mechanism to derail free ligands from association with signaling receptors. Numerical simulations in [20] support this scenario while the analysis of this paper delimit the region in the four dimensional parameter (p, , , v) space favorable to the occurrence of such robust multi-fate signaling gradients.
The presence of abundance of non-receptor with strong affinity for binding with ligand and for rapid degradation of the resulting non-signaling ligand-non-receptor compounds can be a mechanism to derail free ligands from association with receptors to result in robust development of other biological organisms. While the mathematical analysis leading to the delimitation of region in the parameter space favorable to the occurrence of such robust multi-fate signaling gradients may or may not be applicable to other gradient systems, the quantification of multi-fate signaling gradients and robust measures should remain central to robustness studies of the biological developments based on appropriate signaling morphogen gradients.
It is easy to verify that and are, respectively, upper and lower solution of 81. Existence of a solution follows from a theorem of Sattinger (Theorem 2.1 of [27]), and the solution satisfies .
To prove the uniqueness, assume that we have two solutions u 1 (x), u 2 (x). Let w(x) = u 1 (x) − u 2 (x); then w(x) satisfies where by q ,u (x, u) ≥ 0. Hence, we have

Integration by parts gives
With both integrands non-negative, we conclude w(x) 0, for all xin [0, 1].
The following result for the comparison of solutions of two different BVP in differential equation is a direct consequence of maximum principle (Theorem 4.1 of [26]).
For I', we note that the function is an increase function when u ≥ 0 (see the proof below). If ( , v) G p, so that ≤ (p, ), we have when a /(5 + 4a ) ≤ a(x) ≤ 4a /(5 + a ), and hence I' is also satisfied.
The fact that H(u) is an increase function is given below. We have where By h(0) = 0 and we have h(u) > 0 and thus H'(u) > 0.

C.2. Explicit Characterization of G p,γ
Let (83) and v 2 ( ) to be the solution of following equation for v (84) where * is the solution of 49 and is therefore a function of p, , and v. In particular 2 does not depend on v.

Proof
First, for ( , v) in the part of the parameter space originally defined in 50, we want to show that ( , v) is in G p, as specified by 85. In this case, it is evident that ≤ 2 and v ≥ v 1 ( ). We only need to show that v ≤ v 2 ( ). But from Lemmas 3.1 and B.1, we have (86) and, given the upper bound on v in 50, which implies v ≤ v 2 ( ).
Conversely, suppose ( , v) is as specified by 84, then the condition on and the first half of the condition on v in (50) are met. For the remaining upper bound on v, we have from (87) with the last inequality assured by 86 given v ≤ v 2 ( ). Thus, if the parameter pair ( , v) satisfies 85, it also satisfies 50.
To determine the set G p, for a given pair of (p, ) by Proposition 9, we only need to find v 1 ( ) and v 2 ( ) for any [0, 2 ] from the second relation in 83 and 84, respectively. Sample regions for several pairs of (p, ) are given in Figure 2.  The parameter ranges for multi-fate gradient and robust multi-fate gradient with given p and (the points with robust multi-fate gradient are marked by dots). ) (or ( , )) from the region and any ( , v) acceptable to (p, ), the system is not robust for the two fold of ligand synthesis (R v (p, , v) > 0.2). The region of 'R' means that for any (p, ) (or ( , )) from the range, there exist ( , v) that acceptable to (p, ) such that the system is robust for the two-fold increase of ligand synthesis (R v (p, , v) ≤ 0.2). In (b), the circles are original data from the simulations, and the dashed lines are fitted values through 70 and 78, respectively. The original data and fitted values of the function C b (p, , ) in 76. Table 1 Parameter range for good robustness. , v) such that the robustness is less than 0.2, we consider otherwise, . The values given here are maximum of those pairs for given p. Note that in the limit of p 0, tends to a finite limit mathematically but is not biologically realistic.
Discrete Continuous Dyn Syst Ser B. Author manuscript; available in PMC 2013 October 04.