MORPHOGEN PRODUCTION AND MORPHOGEN GRADIENT FORMATION

Partial differential equations and auxiliary conditions governing the activities of the morphogen Dpp in Drosophila wing imaginal discs were formulated and analyzed as Systems B, R, and C in [7] [9] [10]. All had morphogens produced at the border of anterior and posterior chamber of the wing disc idealized as a point, line, or plane in a one-, two-, or three-dimensional model. In reality, morphogens are synthesized in a narrow region of finite width (possibly of only a few cells) between the two chambers in which diffusion and reversible binding with degradable receptors may also take place. The present investigation revisits the extracellular System R, now allowing for a finite production region of Dpp between the two chambers. It will be shown that this more refined model of the wing disc, designated as System F, leads to some qualitatively different morphogen gradient features. One significant difference between the two models is that System F impose no restriction on the morphogen production rate for the existence of a unique stable steady state concentration of the Dpp-receptor complexes. Analytical and numerical solutions will be obtained for special cases of System F. Some applications of the results for explaining available experimental data (to appear elsewhere) are briefly indicated. It will also be shown how the effects of the distributed source of System F may be aggregated to give an approximating point source model (designated as the aggregated source model or System A for short) that includes System R as a special case. System A will be analyzed in considerable detail in [6], and the limitation of System R as an approximation of System F will also be delineated there.

1. Introduction.Morphogens are molecular substances that bind to cell surface receptors and other kinds of (nonreceptor and co-receptor) molecules.The gradients of morphogen-receptor complex concentrations are known to be responsible for the patterning of biological tissues during the developmental phase of the biological host.For a number of morphogen families (including Dpp in the wing imaginal 240 A. D. LANDER, Q. NIE AND F. Y. M. WAN disc of Drosophila fruit flies), it is well established that the concentration gradients are formed by morphogens transported from a localized production site and bind to surface receptors of cells near and away from the production site (see [2], [3], [4], [18], and other references cited in [7]).Recently, the mechanism of morphogen transport has been re-examined by theoreticians and experimentalists, resulting in considerable uncertainty about the role of diffusion in transporting morphogens and other transport mechanisms suggested as replacements.Observations perceived to be inconsistent with diffusive transport were summarized and addressed in [7] on the basis of results from a quantitative analysis of appropriate mathematical models in the form of a system of partial differential equations and auxiliary conditions (defining an initial-boundary value problem, abbreviated as IBVP) describing the known morphogen activities in the wing imaginal disc of Drosophila fruit flies.The first group of results from this quantitative study reported in [7] shows that diffusive models of morphogen transport (such as Systems B and C in [7]) can account for much of the data obtained on biological systems, including those that have been used to argue against diffusive transport.When observations and data are correctly interpreted, they not only fail to rule out diffusive transport, they favor it.
The mathematical underpinning of the case for diffusive transport of morphogens made in [7] has been provided in [10] and elsewhere.Specific results obtained include a necessary and sufficient condition for the existence of nonnegative and monotone decreasing steady state concentration gradients for free and bound morphogens, characterization of the shape of these gradients and their stability with respect to small perturbations.One remarkable outcome of our analysis shows that the governing boundary value problem (BVP) for the steady state behavior of the more complete model that allows for endocytosis and for receptor synthesis and degradation (designated as System C in [7]) may be reduced to the same BVP for the corresponding model without internalization and with a fixed receptor concentration (System B in [7]), but now with its amplitude and shape parameters, β and ψ, replaced by the relevant "effective" amplitude and shape parameters.A similar reduction to the same BVP was carried out for the intermediate System R, which allows for receptor renewal but not for endocytosis [10].The mathematical equivalence of the steady state problem for these systems justifies the use of extracellular models (without internalization) to simplify the analysis of morphogen gradient formation.
In the point source systems B, C, and R, morphogens are synthesized at the border between the anterior and posterior chamber of the wing disc, with the border idealized as a point, a line, or a plane in a spatially one-, two-, or three-dimensional model.In reality, morphogens are produced in a narrow region of finite width between the two chambers in which diffusion and reversible binding with renewable receptors also take place.The present investigation revisits the point source System R investigated in [10], now allowing for a Dpp production zone of finite width between the two chambers.The theoretical foundation for the more realistic model, designated as System F, will be established and an analytical solution for special ranges of parameter values will be obtained.Preliminary results from System F and related distributed source models to investigate robustness of the developmental process in Drosophila wing discs have been reported in [13] with more detail to be published in the near future.A report on applications of System F type distributed source models to investigate the effect of diffusible nonreceptors such as Sog can be found in [12].
It will also be shown that the extracellular model with a spatially distributed morphogen production, System F, leads to some qualitatively different morphogen gradient features.One significant difference from System R is that there is no longer any restriction on the morphogen production rate for the existence of a steady state concentration of Dpp-receptor complexes.We will also discuss briefly how the effects of the distributed source of System F may be aggregated to give an approximating point source model (designated as the aggregated source model or System A for short) that includes System R as a special case.System A will be analyzed in considerable detail in [6], and the limitation of System R as an approximation of System F will also be delineated there.
2. An extracellular formulation with receptor synthesis.As in [7], we simplify the development of the wing imaginal disc of a Drosophila fly as a onedimensional phenomenon, ignoring variations in the ventral-dorsal direction, and the apical-basal direction, because extensions of the one-dimensional model to account for developments in these other directions are straightforward (see, for example, [9]).To investigate the consequence of spatially distributed morphogen production, we will work with an extracellular formulation similar to System R in [10], where we have shown that the results for such a model may be re-interpreted as the corresponding results for a model where morphogen-receptor complexes internalize (through endocytosis) before degradation.
Let L(X, T ) be the concentration of the diffusing morphogen Dpp at time T and location X in the span from the midpoint of the morphogen production region X = −X min to the edge of the posterior chamber of the wing disc at X = X max , with morphogens produced only in −X min < X < 0. Let R(X, T ) and [LR](X, T ) be the concentration of unoccupied receptors and morphogen-occupied receptors, respectively.For the underlying biological processes of the development described in [7], we add to Fick's second law for diffusive transport of Dpp (∂L/∂T = D∂ 2 L/∂X 2 , D being the diffusion coefficient) terms that incorporate the rate of receptor binding, −k on LR, and dissociation, k of f [LR], with k on and k of f being the binding rate constant and dissociation rate constant, respectively.In living tissues, molecules that bind receptors do not simply stay bound or dissociate; they also (endocytose and) degrade [18].In accounting for the time rate of change of the Dpp-receptor complexes, we allow for constitutive degradation of [LR] by introducing a degradation rate term with a rate constant k deg .There is also a separate accounting of the time rate of change of the concentration of unoccupied receptors as they are being synthesized and degrade continuously in time (with a degradation rate constant k´d eg as in [10]).In this way, we obtain the following reaction-diffusion system designated as System R for the evolution of three concentrations L, [LR], and R: for −X min < X < X max and T > 0, where V L (X, T ) and V R (X, T ) are the rates at which Dpp and receptors are synthesized, respectively.In [10], we were interested only in the portion of the wing disc corresponding to X > 0, where there is no morphogen production (so that V L (X, T ) = 0 for X > 0) with Dpp introduced into the region 0 < X < X max through a point source at the end X = 0. We will discuss the relation between this point source version of System R and the present model, which considers explicitly the activities in the region −X min < X < 0, where morphogens are produced.
With −X min being the midpoint of the Dpp production region, we have by symmetry The far end of the wing disc (i.e., the edge of the posterior chamber) is taken to be a sink so that At T = 0, we have the initial conditions To reduce the number of parameters in the problem, we introduce a reference unoccupied receptor concentration level R0 (to be specified later) and the normalized quantities In terms of these new quantities, we write the IBVP in the following normalized form, for t > 0, and 2.1.Time independent steady state behavior.For the existence of a timeindependent steady state behavior, the morphogen production rate and the receptor synthesis rate must be uniform in time: For the present investigation, we ignore possible feedback effects and, unless indicated otherwise, consider V L to be a step function with V L (X, T ) = VL H(−X) for some constant VL to simplify our discussion.Correspondingly, we have We will also take the nonnegative receptor synthesis rate V R to be with 0 ≤ ρ 2 = Vn / Vp ≤ 1 unless indicated otherwise.In that case, we have where With the initial receptor concentration taken to be the steady state receptor distribution prior to the onset of morphogen production, R 0 (x) = V R (X)/k deg , we take R0 = Vp We are interested in a time-independent steady state solution ā(x), b(x), and r(x) for the system (9)- (12).For such a solution, we may set all time derivatives in these equations to zero to get 20) where a prime indicates differentiation with respect to x .The nonlinear system of ODE ( 19)-( 20) is augmented by the boundary conditions With v L (x) and v R (x) both being piecewise constants, the form of ( 19)-(20) requires that ā(x) and its first derivative be continuous at x = 0.As in the steady state problem for the point source model of System R [10], the two equations in (20) may be solved for b and r in terms of ā to obtain r = α 0 r 0 (x) where The expressions in (22) are now used to eliminate r and b from (19) to get a second order ODE for ā alone: This second order ODE is supplemented by the two boundary conditions (21), keeping in mind also the continuity conditions on ā and ā´at x = 0.For our choice of synthesis rates V L and V R , we have v L = 0 and r 0 (x) = 1 for the range 0 < x < 1 so that In the complementary range −x m < x < 0, we have v L = vL and r 0 (x) = ρ 2 so that for some known value of ρ 2 in the range 0 ≤ ρ 2 ≤ 1.The governing ODE for ā(x) for the present extracellular model is identical in form to the corresponding ODE for System B (with a fixed receptor concentration for the point source case) investigated in [7] and [10].This observation effectively allows us to use simpler systems with a fixed receptor concentration (without receptor degradation and synthesis) to investigate the effects of other morphogen activities such as feedback mechanism and morphogens binding to non-receptors (e.g., proteoglygans) in the future.
2.2.Existence, uniqueness, and monotonicity.For the present model with a finite morphogen production region, ā(0) is determined by the ligand activities within the production region and is therefore not known a priori.The coupling between the morphogen activities in the two regions −x m < x < 0 and 0 < x < 1 (with ā(x) and ā´(x) continuous at x = 0) makes it necessary to consider a single BVP for the entire solution domain −x m < x < 0, which is structurally different from the corresponding BVP for the point source case.As such, we need to consider anew the issues of existence, uniqueness, monotonicity, and stability of the steady state concentration gradients for the new problem.We begin with the following existence and uniqueness theorem: THEOREM 1.There exists a unique set of nonnegative steady state concentration gradients ā(x), b(x), and r(x) for System R characterized by the two-point boundary value problem ( 21)-( 24) and the continuity conditions on ā(x) and ā´(x) at x = 0.
Proof.As in the cases of Systems B, R, and C with a point source, existence of a nonnegative solution of the boundary value problem is proved by producing a nonnegative upper solution and a nonnegative lower solution for the problem so that we can apply a theorem of Sattinger established in [15] (see also [1,16]).Evidently, a (x) ≡ 0 is a lower solution, since For an upper solution, consider . For this a u (x), we have and hence also so that a u (x) is an upper solution for the BVP for ā(x).The monotone method of [15] assures us that there exists a solution ā(x) of the BVP (21)-( 24), with Since a u (x) is already known to be positive for −x m ≤ x < 1, ā(x) must be nonnegative in the whole solution domain.
To prove uniqueness, let a 1 (x) and a 2 (x) be two (nonnegative) solutions and a(x) = a 1 (x) − a 2 (x).Then as a consequence of the differential equation ( 24) for a 1 (x) and a 2 (x), the difference a(x) satisfies the following differential equation: Upon integration by parts, observing continuity of ā(x) and ā´(x), and application of the boundary conditions in (21), the relation above may be transformed into Both integrands are non-negative and are not identically zero; therefore we must have a(x) ≡ 0, and uniqueness is proved.
REMARK 2. Note that unlike the point source case, there is no restriction on the relative magnitude of the (dimensionless) morphogen production rate v0 and the (dimensionless) degradation rate g 0 for the existence of steady state concentration gradients.
COROLLARY 3.For r 0 (x) > 0 in (0, 1), the steady state concentration ā(x) does not attain a maximum or minimum at an interior point of [−x m , 1] and hence is monotone decreasing in that interval.
Proof.First, it is easy to see that ā(x) does not have an interior maximum in the interval 0 < x < 1.If it does at some interior point x 0 , then we must have ā (x 0 ) ≤ 0. But since ā(x) ≥ 0 and hence The continuity requirements imply ā(0) = ā´(0) = 0.But it is impossible for the solution of the ODE (24) to satisfy both of these conditions as well as the boundary condition ā´(−x m ) = 0.
For the range −x m < x < 0, consider first the case ρ 2 = 0 (so that V R = 0 in that region).In that case, we have from the ODE (24) and the boundary condition ā Hence, there can be no interior maximum or minimum for the case ρ 2 = 0.For ρ 2 > 0, the ODE (24) for ā(x) and the reflecting boundary condition at where β = vL /g 0 , or where ām is the unknown value ā(x m ) to be determined by the solution process.If ζβ > ρ 2 , only (26) applies, since we must have ā(x) ≥ 0. Since the solution of the BVP exists, the right hand side of (26) must be nonnegative, and ā (x) must be < 0 to satisfy ā(1) = 0 (given that ā(x) cannot be a constant since (25) does not apply).It follows that ā(x) has no interior (local) maximum or minimum in the entire interval −x m < x < 0. If ρ 2 − ζβ is positive so that (25) applies, then we have keeping in mind the continuity conditions on ā and ā at x = 0.But we must have ā (0) < 0 for the correct solution in the range 0 < x < 1; hence, (25) cannot be a solution for the BVP in −x m < x < 0. This, together with the result for ζβ > ρ 2 , rules out an interior extremum in −x m < x < 0. This completes the proof.
3. Approximate solutions.We see in the proof of Corollary 3 that an exact solution of the BVP for ā(x) is possible.However, such an exact solution is not particularly informative regarding the dependence of solution behavior on the biological parameters of the problem.Accurate numerical solutions can be obtained by a number of methods for the relevant two-point boundary value problem, (21) and (24), or the original initial-boundary value problem, (10)- (13).In this section, some approximate analytical solutions will be obtained both to serve as useful benchmarks for numerical solutions and to provide a more explicit delineation of the effects of the biological parameters.

3.1.
No receptor synthesis in the morphogen production region.For the extreme case where there is no receptor synthesis in the morphogen production region −x m < x < 0 so that ρ 2 = Vn / Vp = 0 (and thereby no concentration of either occupied or unoccupied receptors in that interval), the exact solution of the ODE in that region and the reflecting end condition a The constant of integration c 0 is to be determined through the continuity conditions at x = 0.It turns out that the coupling between the solutions in the two different parts of the solution domain (that of the morphogen production region −x m < x < 0 and that of the nonproduction region 0 < x < 1) is such that we can determine the solution in the region x > 0 without knowing c 0 and then calculate c 0 from the solution obtained.This is because we have which is a known quantity.Because of the continuity of ā´at the junction x = 0, the condition (28) serves as the second boundary condition (in addition to ā(1) = 0) for the ODE This two-point BVP determines ā(x) in 0 < x < 1 ; it can be solved by anyone with the available numerical software.The continuity of ā at x = 0 then determines c 0 to be For ζ 1, an explicit solution for the problem is In particular, we have When ζ is not small compared to 1, we generally have to solve the BVP for ā(x) in the range 0 < x < 1 numerically to find ā(0), which in turn determines the morphogen gradient shape in that range.This has been done, and the results for a 0 and a m for a typical set of parameter values are given in columns 2 and 7 in Table 1.We see in columns 4 and 5 of Table 1 that the corresponding results by the leading term perturbation solution are in excellent agreement for β = vL /g 0 ≤ 1 and still to within a 10% relative error for β = 5.Table 1.ā(0) and ā(−x m ) vs. β = vL /g 0 for ρ 2 = 0.
a m 0.25 0.00159 0.00159 0.00158 0.00166 0.00176 0.00183 0.50 0.00319 0.00320 0.00316 0.00332 0.00334 0.00369 1.00 0.00644 0.00647 0.00632 0.00663 0.00715 0.00744 5.00 0.03488 0.03512 0.03159 0.03317 0.03877 0.03985 10.0 0.07660 0.07738 0.06317 0.06634 0.08511 0.08655 25.0 0.24295 0.25076 0.15793 0.16586 0.27213 0.26783 A significant simplification of the problem for general ζ is possible, however, when X min is small compared to X max .In that case, we may, for a good first approximation, take X max to be infinite.(Correspondingly, we should use some other reference length X 0 , such as X min or the typical span of the posterior chamber, which is about 100 cells deep, instead of X max in ( 7)-( 9).)For this approximation, we should have ā(x) → 0 and ā´(x) → 0 in the limit as x → ∞.With these auxiliary conditions, we get the following first integral of the governing ODE in the range x > 0: or, with With the end condition (28), evaluation of (31) at x = 0 leads to a relation between The right hand side of (32) is a monotone increasing function; hence, given β m µ 2 > 0, there is a unique positive solution of (32), denoted by z p (β m µ 2 ) with ā(0) = a 0 = g r z p (β m µ 2 )/µ 2 .Having z 0 = z p (β m µ 2 ), the ODE (31) can be solved to give z as a function of ξ = µx with β m µ 2 as the only parameter.We summarize the results as follows: THEOREM 4. In steady state, the free Dpp gradient for x > 0 in the case ρ 2 = 0 is given by ā(x) = g r z(x)/µ 2 = g r Z(µx; z 0 )/µ 2 , where z 0 = z p (β m µ 2 ) is the unique positive solution of (32) and Z(µx; z p (β m µ 2 )) is the solution of (31) with the auxiliary condition z(0 The dependence of a 0 and a m on vL (with all other dimensionless parameters fixed) is shown by the results of a typical set of other dimensionless parameters in column 3 and 6 of Table 1.There is an excellent agreement with the exact numerical solution for the range of values of β given in that table.The results of the last two rows of the table are significant, as they correspond not only to β > 1 but also to ζβ > 1; neither is allowed by the corresponding ad hoc point source models of Systems B, C, and R previously studied in [7], [9], and [10].It is also remarkable that for the limiting case of X max = ∞, we have succeeded in determining the complete structure of the steady state gradients without solving any BVP and that the solution of the relevant BVP for the actual concentration gradients depends on only a single parameter β m µ 2 .

Perturbation solution for ζ
1 (ρ 2 > 0).For ρ 2 > 0 , we no longer have an explicit determination of the solution for the range −x m < x < 0, as in the previous subsection.Nevertheless, a perturbation solution in ζ is still appropriate for a small Dpp synthesis rate: For sufficiently small values ζ so that ζā α 0 , the leading term ā0 , determined by the BVP is an adequate approximation of the exact solution.Here, we have, in terms of the Heaviside unit step function The ODE now must be solved for x < 0 and x > 0 separately, and the constants of integration are to be determined by the two end conditions and the continuity conditions on ā0 (x) and ā 0 (x) at x = 0.The exact solution for this linear BVP is where Results for ρ 2 = 1 obtained from these formulas for the same typical set of parameter values used for Table 1 are given in columns 4 and 5 of Table 2 along with accurate numerical solutions of the original BVP ( 21) and (24) for comparison.We see there is excellent agreement for β < 5. Higher-order correction terms of the perturbation series can also be obtained to improve on results for moderate values of ζβ = vL /g r = O(1).
Table 2: ā(0) and ā(−x m ) vs. β = vL /g 0 for ρ 2 = 1.Two other limiting situations are of interest.In the absence of receptor synthesis in the morphogen production region, the perturbation solution obtained in this subsection reduces to those of the last subsection.In particular, we have, as in (29), The other extreme case of interest occurs when X max X min so that x m 1.For this case, we consider the limiting behavior as X max → ∞ with ā0 (x) and its first derivative both → 0 (and with the normalizing length scale taken to be some typical length scale, possibly X min , instead of X max ).In that case, we have For the limiting cases of ρ 2 = 0 and ρ 2 = 1, the approximate expression for ā(0) becomes respectively.
To give some indication on the type of applications of interest to biologists, we consider here the mid-level location x h of the ligand-reception concentration b(x) defined by b(x h ) = 1  2 b(0).It serves as a rough measure of the span, steepness, and convexity of the concentration gradient.With b(x) proportional to ā(x), we have from the expressions (33) and ( 34) COROLLARY 5.At low morphogen synthesis rate, the location of mid-level ligandreceptor concentrations, x h , is given by ( 38) to a first approximation.It is independent of the morphogen synthesis rate and moves toward the origin as µ ω → ∞.
Another application of the approximate solution (33) to determine indirectly the effects of a diffusive non-receptor such as Sog on the gradient shape can be found in [12].A direct determination of the same effects has been found in [5] and [11] to be much more difficult (by an order of magnitude at least).

3.3.
Approximate solutions for large X max (ρ 2 > 0).Given the simplification realized for ρ 2 = 0 , we consider for ρ 2 > 0 also the limiting case of X max = ∞ as an approximation for the case of a very large X max , say X max X min .As before, we take for this limiting case ā(x) → 0 and ā´(x) → 0 as x → ∞.A first integral of the ODE (24) satisfying the reflecting end conditions ā (−x m ) = ā(1) = 0 can be written in the form with {z, z 0 , z m } = ζ{ā, a 0 , a m }/α 0 .Continuity of ā (x) and ā(x) at x = 0 requires which gives one condition for the two unknown parameters z 0 and z m .A second condition may be obtained from solving the ODE (39) in the range where {z, z 0 , z m } = ζ{ā, a 0 , a m }/α 0 as before.The two conditions (40) and (42) are to be solved simultaneously for the two unknown parameters z m = ζa m /α 0 and z 0 = ζa 0 /α 0 .Even without an explicit solution for z(x), we see from the structure of (39) that ā depends on x through ξ = √ 2µx with ā(x) = α 0 Z(ξ; z 0 , z m )/ζ.Once we have a 0 and a m , we can obtain ā(x) from (41) for x < 0 and its counterpart for x > 0 (or by directly integrating (39) for both x > 0 and x < 0 using z(0) = z 0 as the auxiliary condition for the first order ODE).
For the extreme case of ρ 2 = 1, the problem simplifies significantly, as (40) reduces to We now use (43) to express z 0 in (42) in terms of z m , so that the resulting relation is a condition for z m alone.The sequential determination of z m and z 0 effectively completes the solution process, except for the solution of an IVP in the form with ā(x) = α 0 z(x)/ζ.Some typical results of this limiting case of X max = ∞ for the same selected set of parameter values as used previously are given in columns 3 and 6 of Table 2. Again there is excellent agreement between the results for this limiting case and the corresponding exact numerical solutions for the entire range of β calculated, including those with ζβ > 1.It appears that for x m 1, approximating X max by ∞ results in a considerable amount of simplification of the BVP and very little penalty.This simplification was exploited in [12] yielding considerable useful results.

3.4.
High morphogen production rate.For a very high morphogen production rate so that ζ vL α 0 and ζ vL g 0 , we re-scale the unknown by setting ā(x) = vL A(x).The BVP for ā may be written in terms of A(x) as with A and A continuous at x = 0.A leading term approximate solution A 0 (x) for ζ vL g 0 and ζ vL α 0 is determined by the linear BVP with A 0 and A 0 continuous at x = 0.The exact solution of this problem is Since ζ vL g 0 so that ρ 2 ≤ 1 ζβ, the solution may be simplified to If we also have ζβx m 1, we may take A more refined analysis would have ( 46) and (47) applicable only up to some location x c when ā(x) no longer dominates the denominator of (44) and an exponential decay sets in.However, comparison with numerical simulations for the range of parameters of interest showed that our approximate expression (47) is accurate to within 10% for ζβ > 10 and that ( 46) is accurate to a few percents of relative error for ζβ 10.
Consider again the mid-level (half-peak) location x h of the sharp gradient front of the receptor-bound morphogen concentration so that b( COROLLARY 6.For β ω σ 0 /µ 2 ω , the location of the sharp gradient front of the receptor-bound ligand complexes characterized by the location of mid-level (or half-peak) concentration is given to a first approximation by Unlike the low-ligand synthesis rate case, the mid-level location now depends on the magnitude of the synthesis rate with as VL → ∞, as it should be for this case.The biological implications of these results are discussed in [12].

4.1.
A nonlinear eigenvalue problem.In addition to the existence of unique steady state concentrations ā(x), b(x), and r(x), it is important for these concentrations to be asymptotically stable.To investigate the stability of the steady state solution known to exist from Theorem 1, we consider small perturbations from the steady state solution in the form After linearization, the differential equations ( 10)- (11) become The relations (51) and ( 52) are then solved for b and r in terms of â making use of b = g r ā/[g 0 (ā + α r )] to get The expressions (54) and ( 53) are used to eliminate b and r from (50 where where we have set ā(x) = α 0 βm A(x), (57) with A(−x m ) = 1 so that ā(−x m ) = α 0 βm .Note that βm is known from the solution of the steady state problem of the previous section to be positive.Let Together, ( 55) and ( 59) define an eigenvalue problem with λ as the eigenvalue parameter.Though the ODE is linear, the eigenvalue problem is nonlinear, since λ appears nonlinearly in q r (x; λ) so that ( 55) and ( 59) is not a Sturm-Liouville problem.Given that r(x) and b(x) may be discontinuous at x = 0 (because a possible discontinuity in r 0 (x) there), we also have the continuity conditions on â(x)and â (x) at x = 0.In the next subsection, we will show that the eigenvalues of the homogeneous boundary value problem defined by the differential equation ( 55) and the homogeneous boundary conditions (59) must be positive.It follows then that the steady state gradients are asymptotically stable according to linear stability theory.

4.2.
Positive eigenvalues and asymptotic stability.Although it is not necessary to do so, we will work with the smooth function for ε sufficiently small, instead of r 0 (x) to simplify the analysis.In that case, r(x) and b(x) are no longer discontinuous at x = 0 (and we no longer need to be concerned about the continuity of â(x) and â (x) at x = 0).Analogous proofs for the same results working with a more general (possibly discontinuous) r 0 (x) are similar but more tedious.
LEMMA 8.All the eigenvalues of the nonlinear eigenvalue problem (55) and ( 59) are real.
Proof.Suppose λ is a complex eigenvalue and a λ (x) an associated nontrivial eigenfunction; then λ * is also an eigenvalue with eigenfunction a * λ (x), where ( )* is the complex conjugate of ( ).The bilinear relation (which can be established by integration by parts and applications of the boundary conditions in (59)) requires It is straightforward to verify q r (x; λ) − q r (x; λ * ) = −(λ − λ * )Φ(x; λλ * ), where , is a positive real value function of λ.In that case, the condition (61) becomes Since the integral is positive for any nontrivial function a λ (x; λ), we must have λ − λ * = 0. Hence, λ has no imaginary part.
THEOREM 9.All eigenvalues of the nonlinear eigenvalue problem ( 50)-( 52) and ( 59) are positive, and the steady state concentrations ā(x), b(x), and r(x) are asymptotically stable by a linear stability analysis.
Proof.Suppose λ ≤ 0 .Let âλ (x) be a nontrivial eigenfunction of the homogeneous BVP (55) and ( 59) for the non-positive eigenvalue λ.Multiply (55) by âλ and integrate over the solution domain to get After integration by parts and applications of the homogeneous boundary conditions (59), we obtain at least for 0 < x < 1.For any nontrivial solution of the eigenvalue problem under the assumption λ ≤ 0, the right-hand side of (63) is positive, which contradicts the assumption λ = − |λ| ≤ 0 (which makes the left hand side of (63) non-positive).
Hence the eigenvalues of the eigenvalue problem (55) and (59) must be positive and the theorem is proved.

4.3.
Approximate decay rate.Although knowing the eigenvalues to be positive is sufficient to ensure the (linear) asymptotic stability of the steady state morphogen concentration gradients, we want to know the actual magnitude of the smallest eigenvalue to give us some idea of how quickly the system returns to steady state after small perturbations.Since parametric studies require that we repeatedly compute the time evolution of the concentration of both free and bound morphogens from their initial conditions, the value of the smallest eigenvalue will also give us some idea of the decay rate of the transient behavior and thereby the time to reach steady state.
We can also get a (possibly) sharper lower bound for λ s by examining the dependence of λs (η, κ) on η. and κ but will not pursue this analysis here.4.6.Estimation of the decay rate of transients.To illustrate the effectiveness of the various bounds obtained in this section, we compare them with the corresponding decay rate for our distributed source model, System F. This decay rate is related to the time required for the solution of the initial boundary value problem defined by ( 10)-( 13) to reach steady state, denoted by t s (in seconds).An estimate of t s is done by requiring the time for b(0, t s ) − b(0) / b(0) to be ≤ .We then calculate λ s using the formula: for the given relative error .
The above estimation of λ s has been carried out for different sets of parameter values for our problem.In all cases computed for comparison with the corresponding bounds, we fix ρ 2 = 1, k deg = 2 × 10 −4 sec −1 , D = 10 −7 cm 2 sec −1 , X max = 0.01 cm, X min /X max = 0.1, and k on R 0 = 0.01 sec −1 so that g 0 = 0.2 and h 0 = 10.We will vary the remaining three parameters, vL /R 0 = 5 × 10 −5 sec −1 to 5 × 10 −4 sec −1 , k of f = 10 −6 sec −1 to 5 × 10 −4 sec −1 , and k g = 10 −4 sec −1 to 10 −3 sec −1 to examine the usefulness and accuracy of the bounds for the distributed source model.Computation and comparison have also been done for other sets of parameter values.However, unlike the case with System R, we do not know b(0) and have used instead b(0, T ) for a sufficiently large T to achieve a significantly lower relative error than .Since the accuracy in the value of λ s estimated in this way is expected to vary with , the actual λ s reported in Table 3 is an average of the estimates for a range of spanning (10 −4 , 10 −2 ).
In Table 3, comparisons are shown for two values of k g in the biologically meaning range min{g 0 , g r } < π 2 /4 so that λ ≡ min{g 0 , g r } and λ u as defined in Corollary 16 (with λ < λ s < λ c < λ u ) .Our comparison shows a good agreement between the theory and numerical results in all cases shown, including the leading term perturbations when applicable.The significance of the last five cases is that both β = vL /g 0 and β r = ζβ = vL /g r are larger than unity, confirming the decay rate of the transients not significantly different from cases with β ≤ 1 and/or β r ≤ 1 (so that the steady state behavior already known to exist is meaningful).The theoretical results for the present more realistic case of distributed finite Dpp synthesis over a finite region −x m < x < 0 provide us with the assurance that we can meaningfully compute the steady state gradients of interest.However, the presence of the two distinct regions of −x m < x < 0 and 0 < x < 1 with different morphogen activities poses unwelcome tedium to the solution process for the steady state concentration gradients (as well as for the decay rate of transients) of the problem in different ranges of the parameter space.It is therefore desirable to find an appropriate simplification of this more realistic model.One possible approach in this direction would be to reduce the problem to one for a single solution domain with morphogens produced and infused at an end point.Ad hoc models of this type have been developed and analyzed as Systems B, R, and C in [7] and [10].We will indicate briefly in this section how these models may be related to the present distributed source model System F. This is done in [6] by aggregating the activities in the region of the distributed source and some reasonable approximations to get the following boundary condition at x = 0 on ā(x) alone: The condition (83) and the end condition ā(1) = 0 provide the two boundary conditions for the ODE (19) for the determination of ā(x) in the range 0 ≤ x ≤ 1.
The two other concentration gradients b(x) and r(x) for 0 ≤ x ≤ 1 are obtained from once we know ā(x).The corresponding end condition for the time-dependent concentrations is [6]: For the dynamic problem, the concentrations a(x), b(x), and r(x) are to be determined concurrently.For X min X max , we have σ 0 1 so that the flux term appears to dominate the left hand side of (83), and we have as x m → 0 the limiting condition ā (0) = 0 as the boundary condition at the source end.With the morphogen synthesis rate VL held fixed and X min tending to zero, the total concentration of morphogens produced over the entire interval −X min < X < 0 tends to zero, which would constitute an upper bound for the morphogen concentration at X = 0.But a more appropriate comparison with the point source case would be to keep VL X min fixed (with VL X min = V 0 ) so that ā (0) = −v 0 for a fixed v 0 = −v L x m .
Unlike the previous ad hoc formulations, the aggregated source model (System A developed in [6] and) summarized above is an appropriate consequence of the distributed source System F where all approximations involved in deriving the new point source model are known explicitly from the development in [6].System A differs from System R by an additional flux term in the boundary condition at the source end.Such a flux term was neglected in System R (and other previously studied ad hoc point source models [7,10] because its contribution to the final solution was not expected to be significant when the shape parameter ψ is sufficiently large; that is, when the effective binding rate is high.In the aggregate source formulation, the coefficient of this flux term is obtained in terms of known biological parameters of the problem, and the magnitude of ψ needed for the flux term to be insignificant can now be determined quantitatively.System A will be shown in [6] to replicate the main features of the distributed source model System F and delimit the range of applicability of System R.

Conclusion.
A system of partial differential equations and auxiliary conditions is formulated to model the extracellular activities (diffusion, reversible binding, and degradation with receptor synthesis and degradation) of Dpp in Drosophila wing imaginal discs, allowing for distributed morphogen production in a finite region between the anterior and posterior chamber of the wing disc.It constitutes a more realistic counterpart of System R investigated in [10], where morphogens are introduced into the wing disc from an end source at the edge of the posterior chamber (corresponding to the border between the posterior and anterior chamber).Justified by the mathematical equivalence among Systems B, C, and R established in [10], [8], and [9], System F does not deal with internalization of free or bound receptors, but the results may be re-interpreted for applications to development with endocytosis.Existence, uniqueness, monotonicity, and (linear) asymptotic stability of steady state morphogen gradients have been established for the new System F. Approximate analytical solutions have been derived by asymptotic and perturbation methods; their adequacy for the relevant parameter range has been confirmed by accurate numerical solutions.
A surprising and significant finding is a qualitative difference between System F and the point source models analyzed previously: there is no longer any restriction on the rate of morphogen production for the existence of steady state gradients.For the ad hoc point source models, Systems B, C, and R, previously considered in [7,10], to reach a time-independent steady state, the morphogen synthesis rate V L /R 0 must be less than the relevant effective degradation rate.The inclusion of continuous synthesis of receptors and endocytosis does not eliminate this restriction.The requirement on the Dpp synthesis rate is a consequence of requiring the various morphogen concentrations at the boundary point to reach a steady state; no such requirement is imposed on the time evolution of the concentrations.The removal of the biologically unrealistic limitation (leaving no restrictions whatsoever) on the ligand synthesis rate is a welcome improvement, for there appears to be no such a restriction in actual Drosophila wing discs.
To investigate the cause for this important qualitative difference between Systems F and R, we aggregated the effects of morphogen activities in the morphogensynthesizing region and applied the aggregate effects at the border between the posterior and anterior chamber.In this way, we reduce the effects of distributed morphogen source and related activities to a point source at the inner edge of the posterior chamber.An essential difference between the derived point source end condition and the ad hoc end condition of System R (as well as Systems B and C) is the presence of a flux term in the end condition at X = 0 with a coefficient determined in terms of the biological parameters that appear in the model.Having System A allows us to assess the significance of the previously omitted flux term in [6].Among the results obtained for the aggregated source model in [6] is the fact that like System F, System A does not impose any restriction on the Dpp synthesis rate for the existence of a time-independent steady state.
A third qualitatively significant result is the approximate determination of the decay rate of transients for the new System F. The eigenvalues associated with a linear stability analysis of System F are solutions of the conventional eigenvalue problem defined by the ODE (55) and the two conventional homogeneous boundary conditions in (59), one being the Dirichlet type and the other the Neumann type.Unlike Systems B, C, and R, there are no other sources of eigenvalues for System F. Thus, we have the rather tedious task of solving the nonlinear eigenvalue problem (55)-( 59) to obtain the decay rate of the transients.To assist us with the solution process, we establish a number of bounds on the rate to provide a good estimate for the rate in some cases.For Dpp activities in the wing imaginal disc of Drosophila with g 0 < g 0 + f 0 < g r < π 2 /4, we have g 0 < λ s < g 0 + f 0 , which gives a very good estimate for the decay rate without solving the nonlinear eigenvalue problem numerically, since g 0 f 0 in Dpp in wing discs.Numerical solutions of the IBVP confirmed the accuracy and usefulness of the sharp upper and lower bounds for the lowest eigenvalue.For the typical case of min{g 0 .gr } < Λ s π 2 /4, it is remarkable that the slowest possible decay rate, given by our lower bound does not depend on the ligand or receptor synthesis rates but only on the two normalized degradation rate constants.
Preliminary results from applications of System F and related distributed source models to investigate robustness of the developmental process in Drosophila wing discs have been reported in [13], with more detail to be published in the near future.Other applications of System F type and related models can be found in [12].
) then β m = b(−x m ) is positive.(In contrast to Systems B, C, and R, there is no restriction on β m or the rate of morphogen synthesis in the present model, System F.) The ODE for â(x) is supplemented by the boundary conditions â (−x m ) = 0, â(1) = 0. (59)