The Role of Feedback in the Formation of Morphogen Territories

In this paper, we consider a mathematical model for the formation of spatial morphogen territories of two key morphogens: Wingless (Wg) and Decapentaplegic (DPP), involved in leg development of Drosophila. We define a gene regulatory network (GRN) that utilizes auto-activation and cross inhibition (modeled by Hill equations) to establish and maintain stable boundaries of gene expression. By computational analysis we find that in the presence of a general activator, neither auto-activation, nor cross inhibition alone are sufficient to maintain stable sharp boundaries of morphogen production in the leg disc. The minimal requirements for a self-organizing system are a coupled system of two morphogens in-which the auto-activation and cross-inhibition have Hill coefficients strictly greater than one. In addition, the GRN modeled here describes the regenerative responses to genetic manipulations of positional identity in the leg disc. of bound morphogen. In the system, there are no spatially dependent parameters.

The Drosophila leg imaginal disc can be used to illustrate pattern formation, which is the acquisition of positional information by cells. Drosophila has six leg imaginal discs that become patterned during the larval stage and then undergo morphogenesis in the pupal stage to form the six legs of the y. Each leg imaginal disc is a single layered field of cells that integrate positional information from three morphogens, WNT/WINGLESS (WG), DECAPENTAPLEGIC (DPP)/BMP, and HEDGEHOG (HH). Dorsal/Ventral (D/V) positional information is specified by antagonizing inputs from the dorsally expressed morphogen, DPP, and the ventrally expressed morphogen, WG [22,1,9,10,19,20,23,6]. The distal tip of the leg, where the dorsal and ventral territories abut each other is specified by the integration of Wg and Dpp signaling. Abnormal expression of these genes causes patterning defects in flies and developmental defects and cancer in humans.
Since restricted expression of these morphogens is critical for normal patterning, the domain of production of these proteins must be maintained as the disc grows. Here we model the GRN that maintains the mutually exclusive DPP/BMP and WG expression domains during development and reestablishes these territories in response to injury or genetic manipulations. The secreted signals that form this GRN are HH which is expressed in the posterior compartment, WG which is expressed in a ventral wedge in the anterior compartment and DPP/ BMP which is expressed in a dorsal stripe in the anterior compartment (Fig 1). Activation of wg and dpp gene expression, i.e. the production of WG and DPP/BMP, requires HH signaling. WG signaling inhibits ventral production of DPP [22,23,21,6,1,10], and DPP signaling inhibits dorsal production of WG [1,9,20,23], thus setting up a dorsal DPP expressing territory and a ventral WG expressing territory. WG and DPP have also been shown to autoactivate their own production [2,6,30,1,8,19,11,24]. We have generated a computational model for this GRN and used it to test the requirements for auto-activation and cross inhibition in maintaining stable domains of expression.
In this model, the domains of production of WG and DPP have sharp boundaries. The free morphogen can diffuse away from this domain and form a gradient of morphogen signaling. The domain of Wg production (Wp) is a function of WG signaling (WR) and DPP/BMP signaling (BR), where WR represents WG bound to its receptor and BR represents DPP/BMP bound to its receptor. Similarly, the domain of DPP/BMP expression (Bp) is a function of WR and BR. In this system, DPP/BMP acts as a dimer, which results in similar diffusivity coefficients for WG and DPP. This differs from a Turing-Meinhardt model, which predicts that the diffusion constants for two interacting morphogens are different, often by several orders of magnitude [25,17,12]. The objective of the modeling is to determine the elements of a GRN required to establish stable sharp mutually exclusive boundaries of expression of two morphogens from an initial shallow gradient of morphogen concentration. The initial shallow gradient is formed by processes discussed in [14,15]. Once this shallow gradient is formed, the combination of auto-activation and cross-inhibition will cause the gradients to steepen until two distinct territories are formed.
We now compare and contrast some of the recent mathematical models of GRNs in Drosophila. In [13], it is shown that a single diffusible morphogen with localized production and interaction with its receptor can lead stable spatial gradients of morphogen concentration. In [16,18] it is shown that the presence of a competing ligand can increase the strength of the spatial morphogen concentration gradient. In [15], a well defined production region is shown to result in the formation of shallow spatial gradients. In [26] positive feedback acting through a cellsurface bound BMP-binding protein can result in sharp gradients. In this paper, we will consider two morphogens interacting with their respective receptors. Morphogen production rate will be regulated by only auto-activation and cross-inhibition, and both the feedbacks are assumed to be functions of bound morphogen. In the system, there are no spatially dependent parameters.

A model and equations
The Drosophila leg initially develops as a flattened disc with dorsal DPP/BMP and ventral WG production domains, which respectively define dorsal and ventral territories. The disc is also divided into an anterior and posterior compartment defined by the absence or presence respectively of engrailed expression. The center of the disc will form the distal tip of the leg and the outer edge will form the proximal leg structures. At this stage the disc can be treated as a 2D structure. In Figure 1, we show the geometry of a leg imaginal disc, and the domains of gene expression/protein production of wg and dpp/bmp. To simplify the analysis, we consider a one-dimensional system in the dorsal-ventral direction as illustrated in (Figure 1 (b)).
In this system, the following assumptions are made: (a) DPP/BMP and WG are produced in the entire domain. The bio-chemical reactions for receptor ligand interaction are as follows, Applying the law of mass action, the number of receptors available to bind with ligand is RW -Wr and RB -Br for WG and BMP respectively. We will model the unbound morphogens as freely diffusing chemicals and the bound morphogens as spatially fixed. Thus the dynamics of the model are given by the following system of differential equations.
We have assumed no-flux boundary conditions for WG and DPP, because far from the interface, the concentration of the morphogens is roughly constant. The activation function (Act ) smoothly connects a minimum production value to a maximum and is strictly increasing, while the inhibition function (Inh) connects a maximum production value to a minimum production value and is strictly decreasing. These functions will be further defined later (Section 3). Representing the feedback as the product of the activation and inhibition function implies that the auto-activation and cross-inhibition mechanisms are linked.

Scaling
In sections §4.1 we will use asymptotic methods to construct a solution and analyze its stability. In order to perform this analysis, we will need to know the relative orders of each of the terms in (1). We will now scale the variables to normalize the solution and the length of the domain. We thus set, (2) (3) (4) (5) (6) We now define the new paremeters as, Since the diffusivities D ͞ D and D ͞ W are small and equal, we will let ϵ 2 = D ͞ D = D ͞ W and use ϵ as an asymptotically small parameter. With this choice, we expect the thickness of the transition layer to be O(ϵ). The new equations may new be written as, The remainder of the paper will proceed as follows. In §2 we will consider spatially homogeneous equilibrium solutions of (7). In §3, we consider explicit forms of the feedback functions. In §4, we will consider a simplification which reduces the model to two uncoupled second order differential equations. We will show that it is possible to construct a territoried solution with a single morphogen, but such a solution will not be sufficiently robust. We will also show that even if such a solution exists, it will be unstable. In §5 we will give some numerical results and compare with experimental results. Finally, we will discuss the relevance of the results for the simplified system to the full coupled systems and consider further avenues of study.

Construction of Solution
We are interested in steady state solutions with two distinct regions. In one of the regions, corresponding to the dorsal region, b will be a relatively large constant and w a smaller constant value. In the other region, corresponding to the ventral region, w will be the large constant and b will have a smaller constant value. These regions will be connected by a thin layer in which both w and b concentrations will have steep gradients. As a first step in the construction of a solution, we look at the spatially constant steady-states.
Before we can proceed, we must first provide some details for the functions Act and Inh. We let, In the most general case, each instance of Act(u) and Inh(u) may have distinct limits, but as this generalization will only make the notation more cumber-some, we will not consider it here.
We now find the values of the steady state solution in the two regions (dorsal and ventral). In the dorsal region, the value of b r will be relatively high and w r relatively low. Thus, in this region, the terms . In these regions, the solution is constant, so the diffusion terms will all be zero. The steady state equations in the dorsal region are then, (8) (9) (10) (11) The solution to the above system is given by, We may repeat the process in the ventral region to find the ventral steady state values: Note, we will require for a positive solution.
For a steady state solution connecting these two equilibria to exist, each equilibrium must itself be stable in the absence of diffusion. Thus, we set ϵ to zero and linearize about (12) and (13). The eigenvalues of this linearization are given by, For the equilibria to be stable we need the real parts of the above eigenvalues to be negative.
The two restrictions imply that all the steady states must be stable.

Feedback
We now consider an explicit form for the functions Act and Inh. We require that the activation function goes from K min to K max monotonically and the inhibition function goes from K max to K min monotonically, where K max > K min > 0. We will characterize these functions with four parameters. The two parameters K min and K max , have already been discussed. The two remaining parameters will control the sharpness of the transition (from K max to K min ) and the point at which the transition occurs.
The most critical parameter value is the point at which the activation and inhibition switches off/on. If this value is too high or low, the activation/inhibition will have no effect. We can find an appropriate value by considering the values of b r and w r . Since 0 ≤ w r , b r ≤ 1, we may set the switching value to .
A natural choice for the activation and inhibition functions are, (14) ( 15) Here m controls the steepness of the transition and is referred to as the Hill coefficient. u 0 controls the value at which the transition occurs, for our application, we set . This form of activation/inhibition was first considered by Archibald Hill in [7]. The basis for this form of activation/inhibition is considering the simultaneous binding of m ligands to an enzyme to produce a product which will either initiate or inhibit protein production. The fact that simultaneous binding is required makes choices of m large unreasonable [29]. Using a negative Hill coefficient for inhibition is considered in [28].
We will now consider the conditions necessary for the formation of morphogen territories with this feed-back model. A morphogen territoried solution corresponds to a heteroclinic connection between two equilibria of the system. To find the conditions which will allow for such a connection, we consider a simplification. When we consider steady-state solutions, (7) reduces to two second order differential equations which are coupled by the inhibition functions. To simplify the system, we take the inhibition function to be a constant and thus decouple the equations.

Heteroclinic Connection: The simplified system
In this section, we will construct a heteroclinic connection for a system in which the Inh function is a constant. This will reduce the steady state problem for (7) to a single second order ordinary differential equation. We find levels of inhibition which will result in a heteroclinic connection between two states. The existence of a heteroclinic orbit in the situation of a constant inhibition field will require a specific isolated value of inhibition, and thus the heteroclinic will not be robust. Further, we show that any such connection will be unstable.
To simplify our situation, we will consider the case Inh(x) = I where I is some constant. We now look for steady state solutions with a heteroclinic connection. Since w r does not diffuse, we may solve for w r in terms of w and eliminate it from the equation. (16) where . We plug into the steady-state equation for w to get, (17) where, (18) I will leave w r (w) as a function in the equations to simplify the analysis. We note that the function w r (w) : ℝ + → [0, 1) is one-to-one and onto.
We now consider some general results for the existence of heteroclinic connections. For an equation of the form, (19) to have a heteroclinic connection (a front solution) joining the values u = u − to u = u + , the following restrictions on Q(u) must be met: 1. Q(u) must have 3 consecutive roots u − , u 0 , and u + with u − < u 0 < u + .

3.
The constructed heteroclinic orbit will connect the two states u = u − and u = u + as y → ±∞.
The system we are considering is posed on a finite domain, and is of the form so the ϵ 2 w yy + Q(w; I) = 0. First we apply the coordinate change to magnify the region about the interface.
The equation then will be of the form . If we assume that we can satisfy the conditions necessary for the existence of a heteroclinic orbit, the constructed orbit will fail to satisfy the boundary conditions, but only by exponentially small terms. We may thus expect a solution to exist which is exponentially close to the constructed heteroclinic. In §4.1 we will carefully consider the effects of the finite boundary on the stability of the constructed solution. Now we will examine the polynomial in the numerator of Q(w; I), (20) to determine the conditions which will ensure the existence of a heteroclinic orbit. For Q(w; I) to have 3 positive roots, (20)(as a function of w r ) must have 3 roots in (0, 1). This immediately implies that m ≥ 2. This result agrees with numerical observations (see Figure 5). We will assume now m ≥ 2. Descartes' rule of signs implies that (20) has exactly one or three positive real roots (counting multiplicity). If m is even, we have no negative roots and if m is odd we will have exactly one negative root.
If we assume that we have 3 distinct roots in (0, 1), the second constraint will be satisfied due to the sign of the highest power of w r and the fact that the denominator of Q is positive and increasing. We will label the three roots w r− , w r0 and w r+ where w r− < w r0 < w r+ . We denote the corresponding values of w as, w − , w 0 and w + .
We find necessary and sufficient condition for (20) to have 3 distinct roots in (0, 1). We will not give this condition in an explicit form as the resulting expression is cumbersome and provides no illumination. We let w rmax be the value at which P attains its local maximum. (21) We now let Ī be the value of I for which P(w rmax (Ī)) = 0. So for all three roots to be in (0, 1), we require that w rmax (Ī) < 1. If this condition is met then for I > Ī and I − Ī sufficiently small, all three roots of P(w r ) must lie in (0, 1).
We must satisfy one final condition to ensure the existence of the heteroclinic orbit. We need to satisfy the integral condition . We would like to have a large range of I values for which all three roots will be in (0, 1) as this will make it easier to solve the integral condition. The closer w rmax is to at I = Ī, the easier it will be to satisfy the integral condition (this will result in a larger interval of I values for which we will have three positive roots in (0, 1). Finding an explicit restriction on parameter values which are necessary and sufficient for the satisfaction of this condition does not seem possible. The difficulty lies in the fact that changing the value of I to satisfy the integral condition may cause w r+ to move past 1.
We demonstrate that it is quite simple to find a value of I which satisfies the integral condition and for which the three roots of (20) are in (0, 1). We will use the same rate and diffusion constants before but with m = 2, but we can pick K max and K min . To check the integral condition, we make the change of variables . Then the integral condition is, (22) Using numerical integration is sufficient to determine if there are values of I for which the above integral is negative and values for which it is positive. We set K min = 0.001 and K max = 0.051 1 . Using I = 0.01, we have the following zeros of P(w r ), (23) and (24) For I = 0.0097 we get the three zeros of P(w r ) to be, (25) and (26) 1 Although it appears that and thus there shouldn't even be two positive homogeneous solutions, this is not really the case here. The value we need to consider in the inequalities are not Hence there must be a value of I for which the integral condition is satisfied exactly and maintains the root condition. We can solve for the heteroclinic orbit implicitly. We can write (17) as a first order system, This system is Hamiltonian with Hamiltonian function . With this choice of Hamiltonian function, the heteroclinic orbit is given by H(w, u) = 0. Thus, (29) We can solve this equation implicitly to find the following implicit expression for the heteroclinic orbit: (30) The ± determines which heteroclinic we find with the plus sign corresponding to a connection from w − to w + . We can solve this integral numerically to get the profile for the heteroclinic orbit (see Figure 3).
The existence of this solution requires a specific value of I. This implies that the front will be structurally unstable. Any small change in inhibition will mean that the connecting steady-state solution will cease to exist.

Stability of Heteroclinic
We now examine the stability of the heteroclinic orbit constructed in the previous section. We constructed the heteroclinic orbit on an unbounded domain. However, it will fail to satisfy the boundary conditions on any finite domain by exponentially small terms. We can expect a solution to exist on a bounded domain that is exponentially close to the constructed heteroclinic. To study the stability of this orbit, we will construct an eigenvalue problem by linearizing about the constructed solution. We will show the operator associated with the linearization will have an exponentially small eigenvalue. This eigenvalue is related to the translation invariance of the interface when the problem is posed on an unbounded domain. We estimate the eigenvalue in the limit ϵ → 0 and show it is positive, but exponentially small. So although the orbit is unstable, it can persist for an extremely long time.
First we list some asymptotic estimates for the behaviour of the heteroclinic solution. We assume we have a heteroclinic connection given by . Then we have that, Now we construct the associated eigenvalue problem by linearizing about the heteroclinic orbit. (34) and substitute (34) into (17) to get the following eigenvalue problem If we differentiate (17) with respect to y we find . So the operator L posed on an unbounded domain has a eigenvalue λ = 0 with eigenfunction . Since w H is a monotonic function connecting the two equilibria, must be of one sign. Thus it must be the principal eigenvalue. Since fails to satisfy the boundary conditions of (35) by only exponentially, small terms we expect there to be an exponentially small eigenvalue with eigenfunction ϕ 0 exponentially close to . We must determine the sign of the perturbed eigenvalue in order to determine the stability of the heteroclinic orbit.
We now construct this eigenfunction using boundary layer correction terms. In the interior of the domain we expect the eigenfunction to be very close to , the correction terms should be localized to the area near the boundaries. We thus write, We define boundary layer coordinates η − = ϵ −1 (y + 1) and η + = ϵ −1 (1 − y).
The boundary layer correction term will then satisfy, The solution to (38) are given by, Now we can start to estimate the small eigenvalue. First we define the dot product as . We then have the identity (43) We apply this identity to the function and ϕ 0 to get (44) Since in the interior of the domain, we have the following asymptotic estimate, Near x = ±1, we have that . Using all our asymptotic estimates we find that, As can be seen from (46), the principle eigenvalue is positive and thus the heteroclinic orbit is unstable. However the eigenvalue is exponentially small and we can thus expect that a heteroclinic solution may persist for long times.
In this section, we have shown that a heteroclinic solution is possible for a single morphogen with self regulation, but the requirements for the existence make such a system unlikely to exist in a natural setting. The restrictions on the auto inhibition function are far too severe to make this a viable alternative. Even if these restrictions are satisfied, the heteroclinic will be unstable and must eventually collapse.

Numerical Simulations
In the previous section, we found that a steady-state solution in which the two morphogens are confined to distinct regions in not a practical possibility with just auto-activation. We would like to show that the formation of such solutions is possible when in addition to auto-activation, the two morphogens mutually cross-inhibit each others production. The numerical method is based on the method of lines with a central difference on the diffusion terms and a secondthird order adaptive Runge-Kutta time integrator.
There are four pairs of biological parameter values in the model, excluding any for the activation and inhibition, for Wg and BMP respectively. In the numerical simulation, we use the same rate constants and diffusivity for both Wg and BMP. The size of the tissue is x max 0.02 cm ( [3]). The diffusion coefficients for Wg and BMP are chosen to be the same: D b = 1 × 10 −7 cm 2 / s based on the measurement [27]; The rate constants are f w = 1 × 10 −5 /s, hw = 0.12/(sµM), γ w = 5 × 10 −4 /s, and the total receptor concentration is R d = 1µM [14] [18]. For the initial conditions we set the concentration of bound receptor to be uniformly zero. We set the initial concentration of free morphogen to the following linear gradient: where w,ŵ,b̄ and b̂ are defined in (12) and (13).
The parameters we will vary are those dealing with activation/inhibition. Figure 4, Figure 5, Figure 7, Figure 9 are graphs of numerical solutions to (1). In each figure we provide 20 concentration profiles equally spaced in time. Initially we set no bound morphogen and a linearly decaying free morphogen concentration profile. The initial condition is displayed in green, the last time step is displayed in red and intermediate steps are displayed in blue.
In this first set of numerical simulations of (1), we demonstrate the robustness of the solution by repeating the simulation with twice the maximum production rates.
The result that we need m, n ≥ 2 is only proven for a single morphogen with no external inhibition, however numerical results suggest that the requirement m, n ≥ 2 is also true for the full system. Setting just one of the values to 1 results in an unstable front (see Figure 5).
In Figure 7 and Figure 9, we demonstrate simulations with one cross-inhibition and one autoactivation disabled. In both cases the fronts are unstable and move across the domain until the solution is spatially homogeneous.
In Figure 5, Figure 7 and Figure 9, the system has not reached a steady-state by the end of the run. This is done in order to display the initial transient in which the gradient vanishes. In Figure 6, Figure 8 and Figure 10 we repeated the respective simulations with much larger timesteps to demonstrate the convergence to equilibrium.

Modeling Regeneration
Based on this model, auto-activation is insufficient to maintain mutually exclusive territories of wg and dpp production. Stable sharp boundaries of wg and dpp expression are dependent on both auto-activation and cross inhibition. We tested whether this GRN could describe the phenotypes produced by genetically modifying wg and/or dpp production and/or signaling. We will now look at the models ability to capture the behavior of various experimental mutations observed in [6]. These mutations will be mimicked in the model by altering the value of either or both of Act and Inh on specified subintervals of the domain. In all cases, the experimental results are well represented by the appropriately modified version of (1). Each graph represents 20 equally spaced time steps. The initial condition is given in green. The final state is given in red. All intermediate states are displayed in blue. Where available, the results of the mutation are also displayed.
A mitotic clone of dsh blocks Wg signal transduction. In the model, this can be mimicked by absence of the WG receptor in a small region. In the region 0.013 < x < 0.018 we set wr and wrt to 0 ( Figure 11). This can also be modeled by defining WR as inactive or unable to signal ( Figure 12). In the second simulation, the wg receptors in a part of the ventral region can bind with wg, but there is no response. The mutation is in the region 0.01 < x < 0.014 ( Figure 12). When this occurs in a clone of cells in the fly leg, it results in pattern duplications that are predicted by the GRN modeled here. Similarly, the model predicts that production of WG at high levels in the dorsal dpp/bmp domain also produces pattern duplications. For this simulation (Figure 13), Wg production is fixed at a high rate independent of receptor binding for the region 0.012 < x < 0.016. The model predicts patterning duplications that mimic those seen in fly legs where a clone of dorsal cells is mutant for an antagonist of Wg signaling (GSK/ sgg). In these cells, Wg signaling is ectopically activated.

Discussions
As we have seen from the analysis, although it is possible for a single morphogen to form a territory solution based only on auto-regulation, such a situation is unstable and unlikely. However, numerical simulations of two interacting morphogens show that morphogen territory solutions are not only possible, but very robust. The only requirement is that the interaction must be strong. If we use a Hill equation to model the interaction, we require a Hill coefficient of at least two. Such a coefficient would be achieved by requiring two or more morphogen molecules bound to a receptor for activation/inhibition to be turned on.
The study of the full system is still open. Since the diffusion coefficients of the two morphogens are of the same order, it is not possible to separate (7) into two second order equations. It may be possible to force symmetries into the system with an appropriate parameter choices, however there is no numerical evidence that such restrictions are required.    Approximation of the heteroclinic connection for (17). For this simulation, I = 0.0099805992, m = 2, K max = 0.051, K min = 0.001. For these values, w − ~= 0.0000950097 and w + 0 .0075233436. Numerical simulation of full system with parameter values K min = 0.0001 m = n = 3. In both runs the final time is t = 150000 seconds. The plots of w signal and b signal are plots of the production rates of w and b respectively. A numerical simulation of system (1) with m = 1, n = 3, K min = 0.001 and K max = 0.01. The final time step is at t = 100 seconds. A continuation of the simulation given in Figure 5 with a much larger time-step in order to display the convergence to an equilibrium. The final time step is at t = 5000 seconds. A numerical simulation of system (1) for m = 3, n = 3 and cross-inhibition of B r on W shut off. The values of K max and K min are adjusted to produce maximum and minimum production rates of 0.01 2 and 0.001 2 respectively for all terms. Final time step is at t = 300 seconds. A continuation of the simulation given in Figure 7 with a much larger time-step in order to display the convergence to an equilibrium. Final time step is at t = 5000 seconds. A numerical simulation of system (1) for m = 3, n = 3 and auto-activation of W shut off. The values of K max and K min are adjusted to produce maximum and minimum production rates of 0.01 2 and 0.001 2 respectively for all terms. Final time step is at t = 300 seconds. A continuation of the simulation given in Figure 9 with a much larger time-step in order to display the convergence to an equilibrium. Final time step is at t = 5000 seconds. In this simulation there are no Wg receptors for 0.012 < x < 0.016 . The final time step is at t = 15000 seconds. The slight localized increase in W is due to lack of receptor to bind with. The localized increase in B is due to the local lack of inhibition. Simulation of region with defective Wg receptors in the region 0.013 < x < 0.018. Since Wg is still able to bind to its receptor there is no localized increase as seen in Figure 11. The final time step is at t = 15000 seconds. The experimental image is from [6]. Simulation of region with Wg production set to its maximum rate independent of receptor binding for 0.012 < x < 0.016 causing localized increases in W and decreases in B. The final time step is at t = 15000 seconds. The experimental image from [6]. Abbreviations Used W The concentration of free Wg.

W r
The concentration of bound Wg.

B
The concentration of free BMP.
B r The concentration of bound BMP.
R W The total number of Wg receptors.

R B
The total number of BMP receptors.

D W
The diffusivity of Wg.

D B
The diffusivity of BMP.