Modeling and computational analysis of EGF receptor-mediated cell communication in Drosophila oogenesis

Mechanisms of intercellular signaling are crucial in development (Freeman, 2000; Wolpert et al., 1998). By regulating cell differentiation, migration, growth, and death, cell communication guides the development of tissues and organs. In adult organisms, the same mechanisms are responsible for tissue repair and maintenance. Cell communication in development is traditionally studied using a combination of genetic, cellular and biochemical approaches. In a typical sequence of events, developmental genes are discovered, assigned a role in the cell and finally assembled into regulatory networks. It is the connectivity and the dynamics of these networks that determines the mechanisms responsible for the formation of organs as complex as lungs and brains. As a rule, the number of genes involved in any given developmental context is large. Hence, the identification of an essential core subset of these genes sufficient to account for extensive genetic and biochemical data is a difficult task. Furthermore, the quantitative characteristics of the regulatory network, such as the spatial ranges and the time scales of intercellular signals, are rarely known. In several well-studied developmental systems, experiments can be complemented by mathematical models and their computational analysis, in order to test the consistency and the predictions of the proposed regulatory mechanisms (Palsson et al., 1997; Reinitz et al., 1998; von Dassow et al., 2000; Yuh et al., 1998). In this paper we report a model-based analysis of a recently proposed mechanism in Drosophila oogenesis (Peri et al., 1999; Sapir et al., 1998; Wasserman and Freeman, 1998). Our model describes the spatial organization and the dynamics of autocrine and paracrine signaling through the Drosophila Epidermal Growth Factor Receptor (EGFR) (Casci and Freeman, 1999; Nilson and Schupbach, 1999; Perrimon and Perkins, 1997; Schweitzer and Shilo, 1997). We use the model to estimate the relative spatial ranges and time scales of the relevant feedback loops, to interpret the phenotypic transitions in eggshell morphology and to predict the effects of new genetic manipulations.

The mature egg of Drosophila melanogaster is characterized by the presence of two dorsal appendages, a pair organ that supplies the developing embryo with oxygen (Hinton, 1981; Spradling, 1993; Waring, 2000). Dorsal appendages are formed by two groups of cells in the follicular epithelium that overlays the dorsal side of the oocyte (Fig. 1A). These two groups of cells are induced to their appendage producing fate by the interaction between the oocyte and the follicular epithelium (Schupbach, 1987). At stage 8 of Drosophila oogenesis, the nucleus of the oocyte migrates to the dorsal-anterior part of the egg and, using a sophisticated system for mRNA and protein transport and processing, creates a localized source of TGFα-like protein Gurken (Fig. 1B) (Neuman-Silberberg and Schupbach, 1993; Nilson and Schupbach, 1999; Riechmann and Ephrussi, 2001). Gurken serves as a localized input to the Epidermal Growth Factor Receptors (EGFRs), a receptor tyrosine kinase uniformly distributed in the overlaying layer of follicle cells (Nilson and Schupbach, 1999; Sapir et al., 1998). Acting through the canonical Ras-MAPK pathway, EGFR activates several genes in the follicle cells directly stimulated by Gurken. As a result of inter- and intracellular communication within the follicular epithelium, the domain of activation induced by the primary signal is modulated in space and time to produce a biochemical blueprint for the formation of dorsal appendages. Specifically, according to the observations of Wasserman and Freeman, the domain of high MAPK activity is first expanded and then subdivided into two smaller regions; the number of these regions corresponds to the number of dorsal appendages in a mature egg (Wasserman and Freeman, 1998). The spatial modulation of the oocyte-derived Gurken signal determines several downstream events that establish the eggshell morphology (Dobens and Raftery, 2000; Waring, 2000).

We have analyzed the mechanism that converts a single-peaked Gurken input into a more complex pattern of MAPK signaling in the responding epithelial cells. Our interest in this particular event in Drosophila oogenesis mechanism is twofold. First, it provides an excellent example of versatility in signaling, demonstrating how simple stimuli define complex patterns in development (Hogan, 1999; Metzger and Krasnow, 1999; Van Buskirk and Schupbach, 1999). More importantly, this mechanism is arguably one of the best studied at the genetic and biochemical level. It therefore provides a good target for the development of modeling methodologies and experimental tests of modeling predictions.

The oocyte-derived signal is modulated by a network of feedback loops in the follicular epithelium. One positive feedback loop is established when the EGFR, which acts through the Ras/MAPK pathway, induces the expression of rhomboid (Hsu et al., 2001; Mantrova and Hsu, 1998; Ruohola-Baker et al., 1993). Rhomboid is an intracellular protease that, together with intracellular protein Star, is responsible for processing and secretion of Spitz, another TGFα-like ligand of the EGFR (Bang and Kintner, 2000; Hsiung et al., 2001; Lee et al., 2001; Rutledge et al., 1992; Tsruya et al., 2002). The secreted Spitz directly interacts with the EGFR, further stimulating the intracellular MAPK (Sapir et al., 1998; Schweitzer et al., 1995b; Wasserman and Freeman, 1998). Thus, Rhomboid acts as a positive regulator of EGFR signaling (Fig. 1C). This positive feedback, together with another one mediated by a secreted EGFR ligand Vein, both amplifies and spatially expands MAPK signaling induced by Gurken (Peri et al., 1999; Sapir et al., 1998; Wasserman and Freeman, 1998).

The amplified and expanded signal is then downregulated by a number of negative feedback loops in the Drosophila EGFR system (Perrimon and McMahon, 1999). Although there is a large number of negative regulators of EGFR signaling in Drosophila, the original mechanism invokes a single endogenous inhibitor, argos (Wasserman and Freeman, 1998). The expression of argos is induced at high levels of MAPK signaling. Argos is a secreted protein that directly interacts with the extracellular domain of the Drosophila EGFR, and inhibits intracellular MAPK signaling induced by Gurken, Vein and Spitz (Fig. 1D) (Golembo et al., 1996; Jin et al., 2000; Queenan et al., 1997; Schweitzer et al., 1995a; Vinos and Freeman, 2000; Zhao and Bownes, 1999). According to Wasserman and Freeman, the amplification of the Gurken signal by Spitz and Vein, and its inhibition by Argos first expands the domain of MAPK signaling and then splits it into two smaller domains, establishing in this way the two groups of appendage-producing cells. In this mechanism, the two peaks detected in the MAPK signaling pattern (Wasserman and Freeman, 1998) are closely linked and co-localized with the two stripes that have been repeatedly observed in the pattern or rhomboid expression (Queenan et al., 1997; Ruohola-Baker et al., 1993; Sapir et al., 1998).

The mechanism of Wasserman and Freeman links gurken, Egfr, spitz, rhomboid, vein and argos into a system of interconnected feedback loops that are jointly regulated by the EGFR and intracellular MAPK signaling (Wasserman and Freeman, 1998). Although the mechanism is supported by genetic and biochemical data, several important questions remain unanswered. First of all, is the proposed network actually capable of converting single-peaked inputs into persistent two-peaked outputs? An independent measurement of the MAPK dynamics in oogenesis reports patterns that are more complex, and cannot be as straightforwardly correlated with the number of appendages of the mature egg (Peri et al., 1999). Although the necessary role of argos is supported by genetic approaches (Hsu and Schulz, 2000; Wasserman and Freeman, 1998), the role of other known EGFR inhibitors remains to be clarified. Argos is the only secreted inhibitor of EGFR in Drosophila (Casci and Freeman, 1999; Perrimon and McMahon, 1999). If it is the only necessary inhibitor, does this mean that the other negative regulators merely serve to modulate the basic pattern established by the secreted Argos? This must be reconciled with the phenotypic transitions in eggshell morphology induced by changing the levels of other inhibitors (Pai et al., 2000; Reich et al., 1999). Several papers have indicated that the expression of argos is detected later than the relevant changes in the rhomboid expression (Peri et al., 1999; Queenan et al., 1997). Can this be used to argue against the mechanism with a single secreted inhibitor? To summarize, the sufficiency of the proposed mechanism with a single diffusing inhibitor and its consistency with the large body of genetic data remains to be established.

We present a mathematical model of EGFR signaling in Drosophila oogenesis. Based on the computational analysis of the model, we validate the mechanism with a single diffusing inhibitor. We find that the proposed regulatory network is indeed sufficient to convert a single-peaked extracellular input into a two-peaked pattern of intracellular signaling. On the basis of our analysis, we are able to rule out some of the previously proposed patterning mechanisms. Specifically, within the framework of our model, we can show that positive feedback is necessary for the existence of two-peaked patterns. We quantify the relative time and length scales of feedback loops that are necessary to produce a two-peaked pattern. Using our model to interpret the observed phenotypic transitions, we point to the deficiencies of the proposed mechanism. While most of the work on dorsal appendage morphogenesis focused on D. melanogaster, which has only two dorsal appendages, it is known that the wild-type appendage phenotype of related fly species can be more complex (Hinton, 1981). For example, the eggshell of Drosophila virilis has four dorsal appendages. Notably, the spatial distribution of the oocyte-derived Gurken in D. virilis is similar to the one observed in D. melanogaster, indicating that the emergence of a more complex phenotype is determined by interactions within the follicular epithelium (Peri et al., 1999). We find that the mechanism with a single secreted inhibitor predicts that eggshells with more than two dorsal appendages can be observed upon quantitative variation of the parameters of the underlying regulatory network. At least indirectly, this can be used to account for more complex eggshell morphologies observed in related fly species.

We have developed a mathematical model that describes early patterning events in the formation of dorsal appendages on the Drosophila egg. The validity of our model is limited to stages ∼9-10B of Drosophila oogenesis (Riechmann and Ephrussi, 2001; Spradling, 1993). The first modeling assumption is that the follicle cells overlaying the oocyte (Fig. 1A) are static during the events guiding the spatiotemporal pattern of MAPK activity. We model the intracellular signaling and cell communication in the follicular epithelium presented with a localized Gurken signal. We are assuming that the signaling response of a follicle cell is determined by the extracellular signals received by it, neglecting the potential dependence of this response on the developmental state of the cell. The patterning mechanism proposed by Wasserman and Freeman was formulated in one spatial dimension (Wasserman and Freeman, 1998). Hence, as a first modeling attempt, our model is one-dimensional, i.e. it tests the peak-splitting capability of the autocrine network along the dorsoventral axis.

Our model describes the dynamics of signaling patterns in terms of three variables. The first variable, R, represents Rhomboid in the follicular epithelium; the second and third variables, A and S, represent Argos and Spitz in the narrow gap between the follicular epithelium and the oocyte. The whole system is driven by a spatially localized input, G, which represents Gurken generated by the oocyte. Argos and Spitz diffuse with different diffusivities in the narrow gap between the oocyte and the follicular epithelium (Fig. 2B). The generation rates of R and A depend on the level of signaling that is positively stimulated by G and S, and repressed by A. This reflects the positive regulation of argos and rhomboid expression through MAPK that, in turn, is activated by EGFR. The generation of S is stimulated by R, reflecting the role of Rhomboid in processing of Spitz. This results in the presence of a positive feedback through Rhomboid and Spitz, and a negative feedback through Argos (Fig. 2A). Note that we do not explicitly account for the positive feedback by Vein. However, the positive feedback in our model describes localized production and diffusion of stimulatory ligand, effectively combining the two positive feedbacks by Spitz and Vein into one.

The dynamics of each of the three variables (R, S and A) in our model are specified by generation and decay terms. The decay terms are modeled as first-order processes. The generation of the active form of Spitz is positively regulated by Rhomboid; we model it by a linear term. The generation rates of Rhomboid and Argos are characterized by sigmoidal functions of the level of MAPK signaling, representing the threshold behavior observed in the expression of the corresponding genes. The level of MAPK activity itself is a strongly nonlinear (sigmoidal) function of the level of EGFR activation that is defined jointly by Gurken, Spitz and Argos. In the model, the EGFR activation level is positively correlated and linear in the concentration of Gurken and Spitz and negatively correlated and linear in Argos. Hence, the level of receptor activation is a function of the extracellular concentrations of its ligands, i.e. it instantaneously responds to the extracellular concentrations of Gurken, Spitz and Argos. This assumption is reasonable because (1) EGFR binding equilibrium is expected to be attained essentially instantaneously (within minutes) on the time scale of the patterning events in oogenesis (hours) and (2) EGFR signaling in oogenesis proceeds in the ligand-limiting regime (Jin et al., 2000; Nilson and Schupbach, 1999; Shvartsman et al., 2001). The steady level of MAPK activity is estimated to be attained within 10-15 minutes (Schweitzer et al., 1995a; Schweitzer et al., 1995b), so it is an algebraic function of G, S and A on the time scale of our problem. These feedbacks lead to compositions of two sigmoidal functions in the generation terms of Rhomboid and Argos. As a composition of two sigmoidal functions is again a sigmoidal function, we model each of these generation terms by a single sigmoidal function characterized by its offset, steepness and amplitude. Importantly, for fixed levels of Gurken and Argos, the homogeneous dynamics of the positive feedback by Rhomboid and Spitz is switch-like. Thus, our network contains three structural elements: the diffusing inhibitor (represented by Argos), the autocrine messenger (represented by Spitz), and the autocrine switch (represented by Rhomboid).

It is observed that the fully developed two-peaked patterns of rhomboid expression have the shape of long, nearly parallel stripes along the midline. Therefore, we will assume that the relevant dynamics can be captured by a one-dimensional model, with the spatial variable corresponding to the circumference of the egg. Furthermore, as the visualized distribution of the MAPK varies smoothly across the pattern, and the characteristic size and distance between the peaks in the gene expression pattern is on the order of several cells, the continuum description of this spatially distributed network is reasonable.

Our model is ‘lumped’ in the sense that it combines numerous biochemical and biophysical processes, such as receptor dynamics, ligand internalization, protein processing, etc., and does not explicitly distinguish between genes and their protein products. However, it is a mechanistic model: in the spirit of the quasi-steady state approximation, it identifies the slowest relevant processes and variables responsible for signaling patterns guiding the formation of the dorsal appendages. We believe that such a one-dimensional continuum model with phenomenological kinetics is a first step in the analysis of patterning mechanisms in Drosophila oogenesis. Further investigations will focus on the effects of additional molecular components, and more faithful description of geometry and transport in oogenesis.

Based on these assumptions, the reaction-diffusion equations for the dynamics of Spitz, Rhomboid and Argos have the following form:

Here, S, R, A and G denote the concentrations of Spitz, Rhomboid, Argos and Gurken, respectively, k’s represent the rate constants for the linear decay terms, D’s represent the diffusion coefficients, g represents the generation rate constants, σ is a sigmoidal function, γ’s and δ’s represent the offsets and widths in the sigmoidal functions, and α and β are the coefficients that characterize the EGFR-mediated input to the MAPK cascade.

We chose the following form for the sigmoidal function describing the generation of Rhomboid and Argos:

The argument of this function is given by the linear combination of Gurken, Spitz, and Argos concentrations: i≡(S−αA+βG−γ_a,r)/δ_a,r, reflecting the strength of EGFR signaling that is determined by the concentrations of its extracellular ligands. Below a threshold determined by γ_r(γ_a) the rate of Rhomboid (Argos) production is identically zero. Above this threshold the generation term is described by a Hill function.

Let us introduce the following constants:

If we then measure lengths and times in the units of L and τ_r, respectively, and introduce the variables:

and dimensionless constants:

we arrive at the following set of dimensionless equations:

Thus, the model contains 11 dimensionless parameters that completely determine the behavior of the signaling patterns. The parameters b_a,r and c_a,r characterize the steepness and the offsets of the generation rates of Argos and Rhomboid, respectively; τ_a,s are the ratios of the time scales of Argos and Spitz, respectively, to the time scale of Rhomboid. The ratio of the length scales of Spitz and Argos is ε, and the relative strength of the negative feedback is λ. The spatial distribution of Gurken is characterized by the width x₀, and amplitude g₀, which is gradually turned on from zero to a steady level on the time scale τ_g. The oocyte-derived Gurken input g(x,t) is described by the following function: g(x,t)=g₀(1−exp(−t/τ_g))”exp”(−x²/x₀²). The width of the Gurken input, x₀, is constant; its amplitude is an increasing and saturating function of time, the time scale on which the amplitude attains its asymptotic value is τ_g. This corresponds to the establishment of the steady value of the oocyte-derived Gurken signal in oogenesis (Nilson and Schupbach, 1999; Peri et al., 1999).

The information available on the biochemical processes involved in patterning of the egg puts a number of restrictions on the choice of the parameters and allows for some quantitative estimates. In particular, the choice of the length and time scales is important. Measurements of patterning of the Drosophila eye, which also involves the EGFR/Spitz/Argos/Rhomboid system, indicate that Spitz diffuses on the length scale of 1-3 cells (each cell is about 4 μm), while Argos diffuses on the length scale of order 10-12 cells (Stevens, 1998), suggesting the choice of L≃40 μm and l≃4-12 μm. This, in turn, means that the ratio ε of these length scales, which is a crucial dimensionless parameter in our model, is ε≃0.1-0.3.

Furthermore, the time scale of the production of Rhomboid is determined by the slow process of gene expression, which takes around 30 minutes. The kinetics of Argos also involves the slow process of gene expression, so we use τ_a=1 (the ratio of the time scale of Argos to that of Rhomboid). However, the generation of Spitz involves release of its membrane-bound form, so it should be a much faster process. In the model, we use τ_s=0.1 for the ratio of the time scale of Spitz to that of Rhomboid.

Another question is the characteristic time and length scales of the Gurken input. The whole process is initiated by the translocation of the oocyte nucleus. In our model, we tried different values of τ_g, for example, in the transient simulation of Fig. 3, we have τ_g=40. Our analysis shows that the pronounced double-peaked solutions exist for x₀≃3, which in the dimensional units is on the order of 120 μm along the circumference of the egg. This is on the order of the size of the follicle cell patch stimulated by Gurken (Nilson and Schupbach, 1999).

Let us emphasize that we do not attempt to fit the parameters of the model to reproduce the experimental observations. The discussion above is only about the orders of magnitude and is intended to serve as a reference point.

Let us now discuss the appropriate choice of the dimensionless model parameters. The first requirement is that there exists a positive feedback with respect to Rhomboid in the absence of Argos and Gurken. The analysis of the model shows that for this b_r has to be sufficiently small. In the model, we chose b_r=0.2. However, in the absence of Gurken the system must have only one stable homogeneous state (the quiescent state). This leads to the requirement that b_a is somewhat smaller than b_r, so the sigmoidal function in the production of Argos should be steeper than that of Rhomboid. In the model, we chose b_a=0.05. The choice of the offsets c_a,r is restricted by the requirements that the threshold of Argos production is higher than that of Rhomboid and that in the absence of Gurken there are no persistent large-amplitude inhomogeneous patterns (this corresponds to the eggshell phenotype with no appendages) (Nilson and Schupbach, 1999). This implies that these offsets should be sufficiently high. However, they cannot be too high, as no patterns will exist in that case; in our computations, we chose c_r=0.4 and c_a=0.5. The choices of the amplitude g₀ of the Gurken input and the strength λ of the negative feedback are also dictated by the ability of the model to produce localized single and double-peaked stationary solutions. This then requires (for our choice of other parameters) that g₀≃1 and λ≃1.6 for the ‘wild-type’ two-peaked pattern.

The nonlinear partial differential equations for the dynamics of Argos, Rhomboid and Spitz were discretized using standard finite difference methods (Iserles, 1997). The resulting large-scale dynamical system was integrated in time using several explicit, semi-implicit, and fully-implicit methods. A stationary pattern computed in the transient simulation was used to initialize a pseudoarclength continuation procedure that can track (both stable and unstable) stationary patterns as functions of the parameters. At each step of the parameter continuation, a Newton-based fixed-point algorithm with a sparse linear system solver was used to converge to a steady state. Upon convergence, the right-most eigenvalues of the linearized problem were computed with inverse power iterations. A point of saddle-node bifurcation, which corresponded to the instability of the stationary pattern, was detected when the leading eigenvalue of the linearization transversely crossed into the right half-plane. The points of saddle-node bifurcations define the boundaries of existence of qualitatively different patterns; these boundaries were computed using AUTO continuation software (Doedel et al., 1991). The outcomes of the transitions following these instabilities were determined by time integration with adiabatically varying control parameter, i.e. the time scale of parameter variation was slow compared with the time scale of relaxation processes in the system. The codes for time-integration, steady-state continuation and stability analyses were implemented in MATLAB, FORTRAN and C. The codes are available from the authors upon request.

Measurements of both the mRNA and protein levels indicate that the spatial distribution of Gurken is steady during stages 9-10B of oogenesis, when patterning leading to the formation of the dorsal appendages takes place (Neuman-Silberberg and Schupbach, 1994; Neuman-Silberberg and Schupbach, 1993; Neuman-Silberberg and Schupbach, 1996; Peri et al., 1999; Peri and Roth, 2000; Queenan et al., 1999). This observation is supported by a simple linear model that accounts for localized ligand release, its extracellular diffusion, binding to cell surface receptors and receptor-mediated endocytosis. Model analysis indicates that the distribution of extracellular Gurken, and hence the input to the follicular epithelium, is equilibrated within 1 hour. This order-of-magnitude estimate is based on binding/”internalization rate constants measured for the EGFR in Drosophila and mammals, and on our estimates of the geometric parameters characterizing the interaction between the oocyte and the follicle cells (Jin et al., 2000; Shvartsman et al., 2001; Spradling, 1993; Waring, 2000).

In our computational analysis, we first determine whether the model network, starting from the state of no stimulation (i.e. zero level of Gurken), can evolve to the steady state with two distinct peaks when the input is switched on. To do so, we determine the dynamic response of our model to the increasing levels of Gurken. For the set of parameters used in the computations presented in Fig. 3A, the monotonic increase in the amplitude of the single-peaked Gurken distribution induces complex spatiotemporal dynamics. In the absence of stimulus, the response of the network is negligible and is, essentially, uniform in space. In the model, this is the result of the thresholds in the generation functions for the production of Rhomboid and Argos. As time goes on, a critical level of the stimulus necessary for the expression of rhomboid is reached. The activation of this positive feedback rapidly increases local Rhomboid concentration and sharpens the response into a narrow region of high Rhomboid/Spitz concentration (Fig. 3A). At the same time, as the expression of rhomboid occurs because of high MAPK activity, this results in the onset of the localized expression of argos in that region as well. However, because of its long-range character, Argos diffuses away from this region and builds up to stabilize a single-peaked pattern. Owing to further increase in the amplitude of the Gurken input, the formation of this localized zone of high Rhomboid/Spitz concentration is followed by its lateral expansion, accompanied by the increase in the level of Argos. The level of Argos is highest in the center of the pattern. As a result, at some point in time the local value of the output in the center is downregulated, and the pattern is split into two smaller regions; the two narrow domains are then stabilized. The final pattern, established for a particular amplitude of the steady single-peaked Gurken input, is composed of two distinct peaks (Fig. 3A). Each peak is characterized by high levels of Rhomboid, Spitz and Argos production. This particular computation confirms that a network with a single diffusing inhibitor is sufficient to convert a single-peaked input into a stable two-peaked output, supporting the mechanism with a single diffusing inhibitor (Wasserman and Freeman, 1998).

The mechanism proposed by Wasserman and Freeman does not specify the eventual outcome of the splitting event. Hence it is not clear whether the experimentally observed two-peaked pattern of EGFR activity is quasi-stationary or dynamic. Our major observation is that the two-peaked pattern is a stable steady state of the model network. At least indirectly, this is in line with experimental observations: a two-peaked pattern of rhomboid expression persists for hours, which is much longer than the time scale of binding and signaling in the EGFR system (Schweitzer et al., 1995a; Schweitzer et al., 1995b). In the following, we will describe the constraints on the quantitative parameters of the network in order for it to be able to establish the stable two-peaked pattern from single-peaked inputs.

The results of the transient simulation are presented in a different format in Fig. 3B. The value of Spitz in the center of the domain is plotted as a function of the maximum concentration of the Gurken input during the transient; qualitatively similar diagrams are produced by plotting the values of Rhomboid and Argos (not shown). Note the presence of two abrupt transitions in this parametric plot. The transition from a uniform distribution to a single peaked-pattern and the following transition between a single-peaked and a two-peaked pattern both occur within narrow ranges of the level of the Gurken stimulus. To characterize the nature of the transitions between the qualitatively different patterns, we have analyzed the steady states of the underlying model.

We have computed the spatially nonuniform steady state distributions of the network components as functions of the input amplitude, parameter g₀ in the model (see Fig. 5A). Again, the value of Spitz in the middle of the domain is used to characterize the pattern. Different branches of this one-parameter diagram correspond to patterns in Fig. 4. Upon the variation of the input amplitude, the low-amplitude steady state shown in Fig. 4A undergoes a discontinuous transition to the steady state characterized by a single, narrow, large amplitude peak (Fig. 4B). This peak is narrower than the inducing Gurken signal. In this case, a broad input is both amplified and sharpened. The transition between the low-amplitude and a single-peaked pattern is accompanied by a hysteresis, i.e. there is a region of inputs where these two types of patterns co-exist. Upon further increase of the input amplitude, a single-peaked pattern undergoes a discontinuous transition to the branch of the steady states with two distinct peaks, shown in Fig. 4C. A circle in Fig. 5A denotes the steady state previously attained in a transient simulation presented in Fig. 3A. Once again, this transition is hysteretic: single- and two-peaked patterns co-exist for a wide range of inputs. Another hysteresis is predicted to occur at a higher value of the inputs between the two-peaked pattern and patterns characterized by a single broad peak (Fig. 4D).

Several conclusions can be drawn from the diagram in Fig. 5A. First, in a large region of input values the patterned steady states in our model are localized: the spatially distributed Gurken input induces outputs that are confined to the middle of the domain. This is a consequence of our choice of the parameter characterizing the threshold of rhomboid expression. This choice is dictated by the observation that the spatiotemporal dynamics guiding the formation of respiratory appendages is constrained to the dorsal-anterior part of the follicular epithelium (Nilson and Schupbach, 1999). Second, patterns obtained upon continuous variation of the input naturally break up into a number of discrete classes; the three classes of interest are characterized by the presence of zero, one and two peaks (Fig. 4A-D). Third, transitions between patterns belonging to different classes are abrupt and discontinuous. In the model, the boundary of existence of a given class corresponds to an instability of a particular steady state (see the Materials and Methods section for the definition of the instability and the techniques for its detection). Overstepping this stability boundary, we obtain a pattern that is qualitatively different from the one that is observed at a parameter value ‘just before’ the instability. Fourth, we believe that the presence of discrete classes of solutions in our model can be used to account for the presence of distinct appendage phenotypes in Drosophila oogenesis. In fact, it is known that the results of multiple perturbations of the genotype, e.g. variations of the amplitude of the Gurken signal, generate phenotypes with zero, one or two dorsal appendages (Nilson and Schupbach, 1999). Within the framework of our model, eggs with no dorsal appendages are identified with the pattern with no large-amplitude nonuniformities (Fig. 4A). It has been established that eggs produced by females with hypomorphic EGFR or reduced levels of Gurken signal have one narrow dorsal appendage (Nilson and Schupbach, 1999; Wasserman and Freeman, 1998). We can identify this phenotype with the pattern characterized by a single narrow peak (Fig. 4B). The wild phenotype is identified with the two-peaked pattern (Fig. 4C). Most recently the phenotype with a single broad appendage was reported for the genotype with increased level of Gurken (Ghiglione et al., 2002). We identified this phenotype with the pattern with a single broad peak (Fig. 4D). The hysteretic transitions between different classes predict that different phenotypes can co-exist for a given range of the Gurken (input) level. In particular, according to the model, decreasing the level of Gurken is predicted to generate a mixture of phenotypes with either zero and one, or one and two appendages. Notably, these changes in the eggshell morphology are indeed observed upon decrease of the stimulus, both at ligand and receptor level (Ghabrial et al., 1998; Schupbach, 1987). Fifth, there is a wide range of inputs for which a slow increase of the stimulatory signal invariably establishes a two-peaked pattern. This can be used to explain the fully penetrant nature of the wild-type phenotype with two appendages.

The mechanism predicts that, once formed, a stable two-peaked pattern will be robust with respect to decreases in the level of Gurken input. This is manifested by a large domain of co-existence of two-peaked and narrow single-peaked patterns. This feature of the mechanism is physiologically significant. It is known that the two-peaked pattern in the rhomboid expression persists beyond stage 10B of oogenesis, which corresponds to the time when the follicle cells start secreting the vitelline membrane that gradually shuts off the oocyte-derived Gurken input (Spradling, 1993; Waring, 2000). This agrees with the predictions of the model: there exists a large hysteresis between the two-peaked and narrow one-peaked patterns as the Gurken amplitude is decreased (Fig. 5A). Hence, the mechanism can explain why the two-peaked patterns are detected even past the onset of the formation of the physical barrier between the oocyte and the follicular epithelium.

Sharp transitions between qualitatively different steady states are a generic feature of systems with positive feedback loops (Ferrell and Xiong, 2001; Shvartsman et al., 2002). In oogenesis, this positive feedback is provided by ligand-induced ligand expression (vein) and processing (rhomboid/spitz). Fig. 5B shows the effect of the strength of the positive feedback on the steady states in our model. Experimentally, it can be realized by changes in the doses of spitz, rhomboid or vein. The qualitative effect of the positive feedback is the same as the effect of the input amplitude. Once again, we have a sequence of hysteretic transitions between the patterns with zero, one and two peaks. The transition to single peaked-patterns at lower levels of rhomboid is supported by the genetic data (Wasserman and Freeman, 1998). At high levels of (Gurken) input or high strengths of the positive feedback (rhomboid/spitz and vein) we predict a transition to the pattern with a single broad peak (Fig. 5A,B; Fig. 4D). The transition to a single broad peak is preceded by a slight separation (up to 20%) and increase in the width of the two peaks. We note that this is at variance with genetic and biochemical results reporting that the two appendages as well as the two stripes in the rhomboid expression pattern move apart when the dose of gurken or rhomboid is increased (Neuman-Silberberg and Schupbach, 1994; Wasserman and Freeman, 1998). It should be noted, however, that a small fraction of eggs laid by females with four extra copies of Gurken had a single broad appendage (Neuman-Silberberg and Schupbach, 1994). Furthermore, recent experiment that used a different technique to increase the Gurken dose reported a phenotype with a single broad appendage (Ghiglione et al., 2002).

The single diffusing inhibitor in the mechanism of Wasserman and Freeman has to be sufficiently strong in order to generate the two-peaked patterns. Decreasing the relative strength of the inhibitor (parameter λ in our model) leads to a pattern with a single broad peak, Fig. 5C. This is consistent with the result of genetic experiments, reporting that decreasing the level of argos, or its upstream activator pointed (an ETS transcription factor), leads to phenotypes with a single broad appendage (Morimoto et al., 1996; Wasserman and Freeman, 1998). Note that the results in Fig. 5C were obtained by decreasing argos level in a spatially uniform manner; this is different from the cited data that were obtained by introducing clones of cells deficient in argos and pointed. We have verified that our results hold in the case when the strength of inhibition is decreased nonuniformly, mimicking the argos and pointed clones.

The mechanism proposed by Wasserman and Freeman relies on the diffusion of both activating and inhibitory ligands of EGFR. However it is unclear whether the inhibitor has to be long-ranged in order to generate the two peaks in the pattern of EGFR activity (Stevens, 1998). Based on the genetic studies of Argos and Spitz interaction in Drosophila eye development (Freeman, 1997), it has been argued that in oogenesis the feedback loops mediated by Argos and Spitz might be long- and short-ranged respectively. Within the framework of our calculations, emphasizing the importance of stationary two-peaked patterns, we have corroborated this hypothesis. Indeed, the spatial range of a diffusing inhibitor must exceed the range of secreted stimulatory ligands by about an order of magnitude in order to produce stable two-peaked patterns, Fig. 5D. In fact, there is a critical value for the ratio of the spatial ranges of positive and negative feedback loops (parameter ε in our model), at which the two-peaked patterns disappear. Importantly, there has been a recent report confirming that the secreted Spitz is localized to a small neighborhood of the producing cells, further supporting that it is a short-range ligand (Bergmann et al., 2002).

Fig. 6 presents the patterns generated in our model upon variation of the input strength and width. For narrow and low amplitude Gurken inputs, the response pattern has low amplitude and is almost uniform in space, the purple area in Fig. 6. This can be interpreted as a requirement on the oocyte-derived signal: it has to be sufficiently wide and strong to pattern the egg. Strongly nonuniform patterns are found outside this area. The majority of patterns fall into a class with a single peak located in the center of the domain (blue area in Fig. 6). Patterns with one and two peaks co-exist in the green area; as we have described above, this co-existence is an immediate consequence of the positive-feedback loop in the mechanism. For inputs of moderate width, the left boundary of the red area in Fig. 6 corresponds to the instability and disappearance of the single-peaked patterns. Crossing this boundary leads to stable two-peaked patterns. Our major finding is the existence of a large area in which a pattern with two distinct peaks is the only output induced by a single-peaked input. In terms of our model, the signaling dynamics guiding the formation of dorsal appendages can be viewed as a sequence of transitions from the purple to the red area in Fig. 6. This can be accomplished by slow monotonic increase of the Gurken signal that is switched on when the oocyte nucleus migrates to the dorsoanterior part of the egg.

Our model also predicts patterns with more than two peaks. In particular, for wider inputs, the ‘wild-type’ two-peaked patterns co-exist with patterns with three peaks. For inputs of fixed width in this area, a monotonic increase of the input amplitude would result in a transition from a one-peaked to a three-peaked pattern, in which the extra two peaks appear at the sides of the original peak. Note that the transition from a one-peaked to a three-peaked pattern exists only for ideally symmetric Gurken inputs. We found that a slight asymmetry in the Gurken profile leads to transitions to two-peaked instead of the three-peaked patterns. This may further corroborate the robustness of the two-peaked patterns.

In the model, the two-peaked patterns arise as a result of the instability of single-peaked patterns. We find that, for strong and broad inputs, two-peaked patterns themselves undergo discontinuous transitions to more complex patterns. Crossing the upper boundary of the red area leads to patterns with four peaks (Fig. 6). We found two different types of stable four-peaked solutions in our model; they correspond to two different instabilities of patterns with two peaks.

Transitions to more complex patterns with three and four peaks were also observed upon variation of other parameters in the model. In particular, a large area of four-peaked solutions has been observed in the two-parameter diagram computed as a function of the input amplitude and the strength of inhibition (not shown). The emergence of more complex patterns is a robust property of our model, and, hence, of the mechanism with a diffusing inhibitor. If the number of peaks in the spatial distribution of receptor activity defines the number of dorsal appendages, then this finding indicates that the mechanism can, at least partially, account for more complex eggshell phenotypes observed in the related fly species, and for some of the mutants of D. melanogaster (Hinton, 1981; Peri et al., 1999; Reich et al., 1999).

We have developed a mathematical model of cell communication in Drosophila oogenesis. Based on the computational analysis of our model, we have validated the recently proposed mechanism of pattern formation by peak splitting (Wasserman and Freeman, 1998). We have found that the proposed network of positive- and negative-feedback loops in the EGFR system can indeed convert a quasistatic, single-peaked input, into stable, two-peaked outputs. Analysis of the parametric dependence of the stationary patterns in the model indicates that the underlying mechanism accounts for a number of experimentally observed phenotypic transitions.

Stationary patterns predicted by the model fall into several qualitatively different classes characterized by the number of peaks in the signaling profiles. The transitions between the classes are discontinuous; this might explain why numerous experiments can be classified in terms of a small number of phenotypes (Fig. 4). Stable two-peaked patterns not only exist in our model, but are also kinetically accessible from the state of zero stimulation and are realized through inputs with a single maximum (Fig. 3). These stable two-peaked patterns are robust for a wide range of network inputs and strengths of the positive and negative feedback loops (Fig. 5). The existence of these stationary patterns requires intermediate input amplitudes and widths, intermediate strengths of the positive feedback, as well as a sufficiently strong and long-ranged inhibition (Figs 5, 6). The fact that large parameter changes induce transitions to qualitatively different patterns is consistent with a large body of genetic data (Nilson and Schupbach, 1999).

Several of the predictions of the original mechanism are at variance with experiments. The first quantitative disagreement is observed in the behavior of the two-peaked patterns upon increases in the doses of the stimulatory signal or the strength of the positive feedback. In the model, these changes produce a discontinuous transition to the pattern with a single broad peak, Fig. 5(A,B); this collapse of the two peaks is preceded by their slight separation. In experiments, eggs laid by the females with extra copies of gurken, heat-shock activated rhomboid, or deficient in Cbl (a negative EGFR regulator) have significantly increased inter-appendage distance (Neuman-Silberberg and Schupbach, 1994; Pai et al., 2000; Wasserman and Freeman, 1998). However, a broad appendage phenotype has been reported in a recent experiment that had used tissue-specific gene expression to increase the level of oocyte-derived Gurken (Ghiglione et al., 2002).

We also note that the initial stage of the transient produced by our model is different from the one measured by (Peri et al., 1999). This suggests a different mechanism at an early stage: EGFR and MAPK are first activated in the large subset of the follicle cells. Even before splitting, this large domain then decreases under the action of non-diffusing EGFR inhibitors with low thresholds (Perrimon and McMahon, 1999).

A more serious problem is related to the role of argos. First, it is unclear why the onset of argos expression is detected much later than the major changes in the patterns of rhomboid and MAPK activity (Peri et al., 1999; Queenan et al., 1997). The second question, based on the analysis of our model, is related to the relative position in the maxima of rhomboid and argos gene expression patterns. According to our model, the maxima in the expression of the two genes should be co-localized; this is an immediate consequence of the fact that both genes require receptor activity for their production (Queenan et al., 1997). At the same time, several independent measurements find that for a while argos is expressed between the two regions of high rhomboid expression. If the resulting two-peaked signaling patterns are quasi-stationary, as suggested by our analysis of the original mechanism, then it is unclear what maintains argos expression in the region of decreased MAPK activity (Queenan et al., 1997; Wasserman and Freeman, 1998). One possibility is that the rate of argos expression is much slower than that of rhomboid, and that the experimentally observed patterns should be interpreted as a transient in the model. We have computationally explored and rejected this mechanism; we have found that making the generation of Argos much slower than that of Rhomboid generates pathologic oscillatory instabilities of the two-peaked patterns. Another possible explanation for the observed relative location of the maxima in the gene expression patterns might involve an additional feedback loop. In this extended mechanism, a yet undiscovered positive feedback regulates the expression of argos, making it insensitive to decreases in EGFR signaling after peak splitting.

We have found that in order for the peak-splitting mechanism to work, the differences in the thresholds of Argos and Rhomboid production must not be too large. In our simulations, the difference in these thresholds is 25%. We have found that in the case when Argos generation is characterized by a significantly higher threshold than that of Rhomboid, only the one-peaked patterns are realized and peak splitting does not occur. In fact, peak splitting requires a rather delicate balance between the spatially distributed stimulation by Gurken and inhibition by Argos.

Our analysis supports the original hypothesis about the differences in the spatial ranges of Argos and Spitz. The cause of this separation of length scales is still unclear. Argos and Spitz are secreted into the gap between the oocyte and the follicle cells, where their transport is accompanied by binding to EGFRs in the follicular epithelium. The strength of ligand/receptor binding can regulate the range of the secreted signal in this situation. If this is the case, then the binding constant characterizing the Argos/EGFR interaction should be less than that of Spitz. Surprisingly, it was found that Argos has a higher affinity for the EGFR (Jin et al., 2000) [but, at the same time, an alternative Argos-like EGF mutant has a lower EGFR-binding affinity (van de Poll et al., 1997)]. Another mechanism regulating the range of secreted ligands relies on their receptor-mediated endocytosis (Narayanan and Ramaswami, 2001; Teleman et al., 2001). Based on the fact that all ligands of mammalian EGFRs are rapidly internalized (Wiley and Burke, 2001), we believe that this mechanism can control spatial ranges of EGFR ligands in Drosophila oogenesis. Furthermore, it is known that Argos interaction with EGFR prevents receptor dimerization and phosphorylation (Jin et al., 2000). As these processes are required to initiate receptor-mediated endocytosis of receptor tyrosine kinases, this further supports the mechanism in which the ranges of Argos and Spitz are controlled by the rates of their internalization.

We must emphasize that our model does not rely on the difference between the diffusivities of Spitz and Argos. This is in contrast to the classical activator-inhibitor mechanism for morphogenesis proposed by Turing (Turing, 1952). Instead, it relies on differences of ranges of these molecules, which depend on a combination of their diffusivities and rates of degradation. This is also true for the Turing mechanism; however, in the latter, the difference between the time constants of the activator and the inhibitor leads to oscillatory instabilities. Moreover, our model is different from the activator-inhibitor models: while Argos plays the role of a long-range inhibitor, the positive feedback is operated by an autocrine switch with a non-diffusing component (Rhomboid) and a short-ranged diffusing messenger (Spitz). Here, it is the time constant of Spitz degradation that determines the range of the positive feedback; however, the existence of the oscillatory instability is determined by the ratio of the time scales of Rhomboid and Argos. Furthermore, it is not difficult to show that with the given relationship between the thresholds of Rhomboid and Argos production (that is, when ca>cr) the Turing-like instability (as a result of the homogeneous increase of the level of Gurken input) is not realized at all in our model. Another difference between our model and that of Turing is that in our model, pattern formation occurs as a result of the instabilities leading to the abrupt formation and transitions between large-amplitude localized patterns. Large amplitude localized patterns, often referred to as autosolitons, are frequently encountered in different nonlinear systems (Kerner and Osipov, 1994; Lee et al., 1994; Meinhardt and Gierer, 2000; Muratov and Osipov, 1996).

In addition to patterns with one and two peaks, the mechanism supports more complex patterns. If the number of peaks in the profile of the receptor activity determines the number of respiratory appendages, then this finding predicts that more complex eggshell phenotypes are to be expected upon quantitative variation of the parameters of the EGFR network. At this point there is a single published observation of eggs with four appendages in mutants of D. melanogaster (Reich et al., 1999). Occasionally, eggs with three appendages are laid by the mutants with defects in Gurken accumulation (T. Schupbach, personal communication). Notably, eggs with four and even six appendages represent wild-type phenotypes of several related fly species: the eggs of D. virilis have four appendages, while the eggs laid by the flies of subgenus Pholadoris have six appendages (Hinton, 1981). According to our model, eggs with multiple appendages can be generated by increases in the width and amplitude of the stimulatory signal (Fig. 6). Further experimental, modeling and computational studies are required to check whether these more complex phenotypes can be realized in oogenesis or whether they manifest a pathological feature of the mechanism with a single diffusing inhibitor. Eggshells with more than two dorsal appendages were reported in the experiments with binuclear Drosophila oocytes (Roth et al., 1999). The two-peaked distribution of Gurken mRNA detected in binuclear oocytes may generate a wide distribution of Gurken protein. As discussed above, even a single-peaked broad input can lead to the formation of extra dorsal appendages.

Within the framework of our model, we have also analyzed a previously proposed mechanism that does not involve positive feedback loops (reviewed by Hsu and Schulz, 2000; Nilson and Schupbach, 1999). In this mechanism, a high concentration in the peak of a broadly distributed Gurken input induces localized expression of argos; Argos is then ‘subtracted’ from the EGFR response induced by Gurken, establishing in this way a response with two peaks. In our model, this mechanism would correspond to the absence of the positive feedback by rhomboid and spitz. This simplifies the model to a single equation that describes the induction of argos by the spatially distributed Gurken. For a wide class of input profiles, we were able to prove that neither the steady state nor the transients admitted by this equation can produce two peaks in the output (in this case, the difference of Gurken and Argos profiles, which is amplified in the MAPK signaling profile). Based on this, we rule out the mechanisms without the positive feedback.

We have tested the proposed peak-splitting mechanism in one spatial dimension. As the patterning events in the follicular epithelium are intrinsically two-dimensional, the next round of modeling should account for an additional spatial dimension. Our preliminary work indicates that, in two-dimensions, an oval-shaped Gurken input leads to a pattern in which the region of high signaling (the distribution of Rhomboid) has the shape of a hollow ellipse. A one-dimensional cross-section of this pattern produces a two-peaked distribution of signaling activity. Additional work is needed to determine the mechanisms that convert this pattern into a two-striped pattern detected in oogenesis. Note that we have analyzed the patterning capabilities of the Rhomboid/Spitz/Argos network in the presence of a single gradient (Gurken). At the same time, other morphogenes might provide additional inputs specifying the formation of dorsal appendages (Dobens and Raftery, 2000; Riechmann and Ephrussi, 2001). In particular, the anterior posterior gradient of Dpp is crucial for positioning of dorsal appendages (Dobens and Raftery, 2000; Peri and Roth, 2000).

The highly conserved nature of EGFR signaling makes it a good testing ground for the modeling and computational analysis of ‘subroutines in development’ (Hogan, 1999). The models should aim to directly interpret and predict the phenotypes arising from complex interactions in the underlying regulatory networks. If successful in complementing and extending the information obtained using genetics, cell biology and biochemistry, our model of EGFR signaling in Drosophila oogenesis can be used as a starting point for the quantitative analysis of the role of EGFR in other stages of Drosophila development (Buff et al., 1998; Halfon et al., 2000; Martin-Blanco et al., 1999; Wessells et al., 1999) and in local tissue regulation in higher organisms (Doraiswamy et al., 2000; Kashimata et al., 2000; Umeda et al., 2001; Wessells et al., 1999; Wiesen et al., 1999).

The authors are indebted to Trudi Schupbach (Princeton) for many helpful discussions during the formulation of the model and for her critical reading of the manuscript. We are very grateful to Eusebius Doedel (Concordia University) for his help with setting up the continuation software. S. Y. S. acknowledges helpful discussions with Alan Michelson (BWH), Laurel Raftery (MGH), Matthew Freeman (MRC) and Uri Abdu (Princeton). We are grateful to Allan Spradling (Carnegie Institute) for helpful discussion and the reference to the book by Hinton. During the initial stages of this work S. Y. S. was supported through the F32-GM20847 NIGMS postdoctoral fellowship in Quantitative Biology at the Department of Chemical Engineering at MIT. C. B. M. was partially supported by the NJIT SBR program.