Similar to other organisms, *Drosophila* uses its Epidermal Growth Factor Receptor (EGFR) multiple times throughout development. One crucial EGFR-dependent event is patterning of the follicular epithelium during oogenesis. In addition to providing inductive cues necessary for body axes specification, patterning of the follicle cells initiates the formation of two respiratory eggshell appendages. Each appendage is derived from a primordium comprising a patch of cells expressing *broad* (*br*) and an adjacent stripe of cells expressing *rhomboid* (*rho*). Several mechanisms of eggshell patterning have been proposed in the past, but none of them can explain the highly coordinated expression of *br* and *rho*. To address some of the outstanding issues in this system, we synthesized the existing information into a revised mathematical model of follicle cell patterning. Based on the computational model analysis, we propose that dorsal appendage primordia are established by sequential action of feed-forward loops and juxtacrine signals activated by the gradient of EGFR signaling. The model describes pattern formation in a large number of mutants and points to several unanswered questions related to the dynamic interaction of the EGFR and Notch pathways.

## INTRODUCTION

Epidermal Growth Factor Receptor (EGFR) regulates a multitude of processes in development. Eyelid morphogenesis in mice and feather array specification in chickens are just two examples from a long list of developmental events that depend on EGFR signaling (Atit et al., 2003; Mine et al., 2005). Some of the most advanced studies of EGFR signaling have been carried out in *Drosophila*, which, similar to other organisms, uses its EGFR pathway throughout development (Shilo, 2003). *Drosophila* EGFR signaling is involved in patterning of essentially all tissue types and compartments in the embryo (Golembo et al., 1996; Szuts et al., 1998; Yagi et al., 1998; Parker, 2006). In larval and pupal stages, EGFR controls growth, patterning and morphogenesis of multiple organs and structures, including wing veins and photoreceptor arrays (de Celis, 2003; Roignant and Treisman, 2009).

During oogenesis, a gradient of EGFR activation patterns the follicle cells, which form an epithelium over the developing oocyte (Fig. 1A). This gradient is induced by Gurken (Grk), a TGFα-like ligand secreted from the oocyte (Schupbach, 1987; Neuman-Silberberg and Schupbach, 1993; Neuman-Silberberg and Schupbach, 1994; Pai et al., 2000; Goentoro et al., 2006; Chang et al., 2008). Patterning of the follicle cells provides inductive cues for body axes specification in the early embryo and initiates the formation of two respiratory eggshell appendages (Van Buskirk and Schupbach, 1999; Roth, 2003). Each of the appendages is derived from a primordium comprising a patch of ~55 cells expressing the transcription factor *broad* (*br*) and an adjacent stripe of ~15 cells expressing *rhomboid* (*rho*), a protease in the EGFR pathway (Berg, 2005; Ward and Berg, 2005). The *br*- and *rho*-expressing cells form the upper and the lower parts of eggshell appendages, respectively (Fig. 1B-F).

Several conceptual and mathematical models of eggshell patterning have been proposed in the past (Wasserman and Freeman, 1998; Shvartsman et al., 2002; Boyle and Berg, 2009; Lembong et al., 2009; Zartman et al., 2011). Although some of these models can explain the origin of the two-dimensional pattern of *br*, none of them fully accounts for the highly coordinated expression of *br* and *rho*. In particular, the connection between the long-range gradient of EGFR activation by Grk and a single-cell wide stripe of *rho* expression is unclear. To address this issue, we synthesized the existing information into a revised model of follicle cells patterning. Based on the computational analysis of this model, we propose that *rho* expression is generated by sequential action of feed-forward loops and juxtacrine signals set in motion by the Grk gradient.

## MATERIALS AND METHODS

### Computational screen

In the computational screen for parameter values (Fig. 2D), a set of ODEs (Eqns 1, 2, 3, 4 and 5) was integrated for 10 dimensionless time units with zero initial conditions for all variables. Juxtacrine signal arriving in a given cell is given by arithmetic mean of stimuli (G_{5} values) in its neighbors in the lattice. Grk distribution in the two-dimensional model was based on a previously published model, adapted to Cartesian coordinates (Goentoro et al., 2006; Zartman et al., 2011). To evaluate the inductive inputs distribution (Fig. 3B), a reaction-diffusion problem given by Eqn 6 was solved at steady state with appropriate boundary conditions (zero flux, except for the anterior boundary for I_{2}, where a unit flux is applied). A standard finite elements method was used. The set of ODEs given by Eqns 1, 2, 3, 4 and 5 (this time with two-dimensional inputs I_{1} and I_{2}, as explained in the Results) was integrated for 10 dimensionless time units with zero initial conditions for all variables. Note that a separate set of equations is solved for each cell on a hexagonal lattice, as each cell has its own pair of values for I_{1} and I_{2}.

### In situ hybridization and immunohistochemistry

Description of the fluorescence in situ hybridization protocol and dpERK staining are provided elsewhere (Fuchs et al., 2011; Zartman et al., 2009). Digoxigenin- and biotin-labeled antisense probes were transcribed in vitro, and purified with RNeasy mini kit (Qiagen). Primary antibodies used were sheep anti-digoxigenin (3:500, Roche Applied Science), mouse anti-biotin (1:100, Jackson ImmunoResearch) and rabbit anti-dpERK (1:100, Cell Signaling), followed by Alexa Fluor-conjugated secondary antibodies (1:500, Molecular Probes). Confocal images were acquired using a Leica SP5 confocal microscope, and processed with ImageJ (Rasband, US National Institutes of Health, Bethesda, Maryland, USA, 1997-2011). Eggshell images were obtained with a Hitachi TM-1000 table top scanning electron microscope.

## RESULTS

### Description of the core mechanism

In this section we describe the network that can convert a graded signal into an ordered arrangement of the expression domains of three genes, similar to that observed during the patterning of the follicle cells by EGFR (Fig. 2A,B). We start by considering a row of cells, which can be thought of as a result of a cut along the follicular epithelium (arrow, Fig. 2A). The lattice is presented with a graded input, I, which models the Grk gradient (Fig. 2C). The state of each cell in the lattice is characterized by expression levels of five genes, G_{1}-G_{5}. G_{1} models a high threshold-target of Grk, such as *sprouty* (*sty*) or *pointed* (*pnt*). G_{2} corresponds to *rho*. G_{3}, which corresponds to *br*, is activated by intermediate levels of Grk. In addition to G_{1}, G_{2} and G_{3}, two more components are needed to generate a configuration in which a G_{2}-expressing cell is positioned between broader expression domains of G_{1} and G_{3} (Fig. 2C).

We attempted to derive the simplest network that may generate the specific pattern of *rho* expression (Fig. 2A). G_{1}, G_{2} and G_{3} are essential components, which are well established by experiments. G_{4} is required to reproduce the initial phases of *rho* expression dynamics, of initially wide domain that later reduced to the ‘T’-shaped domain. G_{5} is essential to generate the single cell-wide rows of *rho* expression, as explained below.

The proposed regulatory network comprises several incoherent feed-forward loops, network motifs in which an input activates both the target gene and its repressor (Alon, 2006). First, an inductive signal (I) activates G_{3} and its repressor, G_{1}. Second, G_{3} represses G_{4}, which is activated by I. Third, G_{4} represses G_{2} and activates G_{5} that in turn activates G_{2}. Finally, the inductive signal I activates another repressor of G_{3}, G_{5}. Regulatory interactions in this network can be either cell-autonomous or juxtacrine, as denoted by solid and broken lines in Fig. 2B, correspondingly. For cell-autonomous interactions, the product of a given gene affects the state expression levels of other genes in the same cell. For juxtacrine interactions, the effect is on the nearest neighboring cells.

Some of the network interactions are deduced from experiments. The feed-forward loop comprising I, G_{1} and G_{3} is based on the fact that high levels of EGFR signaling activate both *br* and its repressors (Morimoto et al., 1996; Deng and Bownes, 1997; Boisclair Lachance et al., 2009; Zartman et al., 2009). As a result, *br* is expressed only at intermediate levels of EGFR activation. Second, Br (the protein product of G_{3} in the model) represses genes such as *Fas3*, which, in the wild type, is expressed only at high levels of EGFR signaling (Ward and Berg, 2005). This supports the feed-forward loop formed by I, G_{4} and G_{3} that generates the initially wide domain of G_{4} expression, which is later reduced to the ‘T’-shaped domain.

The incoherent feed-forward loop composed of G_{4}, G_{2} and G_{5} mimics the specific dynamics of *rho* expression, i.e. of initially wide domain that is later refined to the ‘L’-shaped single cell-wide rows. The product of G_{4} induces G_{5}, which represses G_{3} and activates G_{2} in the neighboring cells. This juxtacrine mechanism is based on the fact that the Notch pathway, which commonly acts in a juxtacrine mode, is active at this stage of oogenesis in the dorsal midline (Ward et al., 2006). Furthermore, Notch signaling was shown to repress *br* and activate *rho* (Ward et al., 2006). Based on this, cell-autonomous loss of G_{3} in our model will lead to ectopic expression of G_{4}, which in turn will lead to ectopic expression of G_{2}. This is consistent with the results of experiments that demonstrated that loss of *br* can lead to ectopic expression of *rho* (Ward et al., 2006). Thus, the proposed network is largely consistent with previous observations. To summarize, G_{1} represents a high threshold EGFR target that represses *br*, G_{2} represents *rho*, G_{3} represents *br*, G_{4} is a hypothetical component that is repressed by *br* and induces G_{5}, which models the Notch pathway.

We cannot formally prove that the proposed mechanism is the simplest mechanism leading to the ordered arrangement of *br*- and *rho*-expressing cells. At the same time, we could not come up with a model that would have a smaller number of components or one that would not involve sequential action of feed-forward and juxtacrine loops.

### Mathematical model

In the one-dimensional model, the graded input follows an exponential distribution: I=exp(−ϕζ), where 0<ζ<1 is a dimensionless coordinate which spans 20 cells (Fig. 2C). The expression levels of five genes in a given cell *k* are denoted by . For each of these variables, the dynamics is governed by the sum of a linear degradation and a nonlinear generation function describing the combined effects of regulatory signals.

In these equations, τ denotes dimensionless time. For simplicity, we assumed that all species have the same decay rate. Following previous models of gene regulation networks (Plahte, 2001; Bolouri and Davidson, 2002; Albert and Othmer, 2003; Thieffry and Sanchez, 2003; Longabaugh et al., 2005; Perkins et al., 2006), gene expression is approximated by Heaviside functions. The Heaviside function is equal to one when its argument is greater than zero and is equal to zero otherwise. Thus, the rates of gene expression in our model assume only two values, zero and maximal. Switches in the rate of expression of a gene in a given cell happens at particular thresholds denoted by θ* _{ij}*. In this notation,

*i*and

*j*denote a target gene and its regulator, respectively (index

*0*corresponds to the inductive input I).

For example, G_{4} is activated by the inductive input *I* and repressed by G_{3}. To describe this dual regulation, the rate of G_{4} production is given by , where is a Heaviside function of the level of the inductive signal (*I*) minus the threshold (θ_{43}), while is the Heaviside function of the threshold (θ_{43}) minus the level of G_{3}. The Heaviside function with superscript + is equal to one when its argument is positive and zero otherwise, whereas the opposite is correct for the Heaviside function with superscript −. Similar to our earlier work on epithelial patterning (Pribyl et al., 2003; Reeves et al., 2005), the effect of a graded input on a specific cell is given by the integral of the input over this cell; *I ^{k}* denotes integral for a cell

*k*.

### Selection of parameter values

To analyze the model, we needed to specify the values of its ten parameters. We took advantage of previous studies of the Grk gradient and published gene and protein expression images that provide information about the spatial distribution of several model components. First, based on the earlier estimates of the spatial distribution of the Grk gradient, the length scale of the inductive signal, denoted by ϕ in the model, is set to 3 (Goentoro et al., 2006). This leaves nine parameters, corresponding to the thresholds of gene regulation functions. As the input levels and values of all model variables vary between zero and one, these thresholds also vary in the same range. Importantly, most of the thresholds can be chosen based on the experimentally observed gene expression patterns.

Based on the patterns of high threshold EGFR-target genes (Morimoto et al., 1996; Yakoby et al., 2008a), we chose the G1 activation threshold in such a way that G_{1} is activated in three cells corresponding to the high levels of the input. Given the distribution of *I*, this leads to θ_{10}=0.65. Second, it is known that both *rho* and *br* can be potentially expressed in the same cells, which correspond to the intermediate levels of EGFR activation (Peri et al., 1999). In our one-dimensional lattice, this corresponds to 11 cells (Fig. 2C). Based on this, the activation thresholds of G_{3} (which models *br*) and G_{4} (which activates G_{2}, which in turn models *rho*) are given by θ_{30}=θ_{40}=0.2.

The threshold of repression of G_{3} by G_{1} (*br* repression by Pnt) is set to a low value, θ_{31}=0.1, ensuring that *br* is repressed without a considerable delay following the induction of its repressor (Yakoby et al., 2008b; Boisclair Lachance et al., 2009; Fuchs et al., 2011). Similarly, the activation thresholds of G_{5} by G_{4}, and G_{2} by G_{5} are chosen to be low: θ_{54}=θ_{25}=0.1. This ensures that G_{2} is initially expressed in the same cells as G_{3} and G_{4} (Peri et al., 1999). In addition, repression of G_{2} by G_{4} is set to have a high threshold, θ_{24}=0.9, providing a delay of G_{2} repression in cells exposed to high levels of inductive signal. This is based on the fact that *rho* is initially expressed in a large region, overlapping the domains of its repressors (Peri et al., 1999).

Selecting the values for the remaining two thresholds, θ_{35} and θ_{43}, corresponding to repression of G_{4} by G_{3} and of G_{3} by G_{5} is less straightforward. However, given the simplicity of the model and low computational cost of solving it numerically, the values for remaining two thresholds can be constrained computationally. To do this, we solved the model for multiple pairs of θ_{35} and θ_{43}, chosen on a uniform grid (Fig. 2D), with all other parameters selected as described above (Fig. 2B). For each of the resulting solutions, we checked whether G_{2} is initially expressed in a domain that overlaps the expression domains of G_{1} and G_{3}, which is later refined to a single cell located between the expression domains of G_{1} and G_{3} (Fig. 2C). Following this procedure, we could identify the values of θ_{35} and θ_{43} that satisfy these criteria (Fig. 2D).

### Qualitative dynamics in one dimension

For each of the points inside the identified parameter region, the time-dependent solution is qualitatively similar to the representative examples shown in Fig. 2C,E. Initially, G_{1}, which is a high-threshold target of the inductive signal, is expressed in a narrow region, while G_{2} and G_{3} are expressed in a wider domain. The G_{2} expression domain is then refined by the repressive effects of G_{3} and G_{4}, while G_{4} is repressed by G_{3}. Note that G_{4} is both a repressor and an indirect activator of G_{2} (Fig. 2B).

In our model, the activating and repressive effects of G_{4} on G_{2} are separated in time, by setting low and high thresholds of activation and repression, respectively. In cells where G_{4}, which is essential for G_{2} activation (via G_{5}), is repressed by G_{3}, G_{2} is also repressed. This corresponds to the fact that late *rho* expression is abolished in *br*-expressing cells (Peri et al., 1999). Thus, *br* is repressed at the posterior border of cells expressing G5, which models the effect of the Notch pathway. The juxtacrine repression of G_{3} by G_{5} ensures that G_{3} is removed from cells expressing G_{2}. This agrees with the fact that *br* is not expressed in cells expressing *rho* at late stages in the patterning process (Dorman et al., 2004; Ward et al., 2006). In order to obtain the qualitative dynamics represented by Fig. 2C,E, the repression of G_{3} by G_{5} should be delayed, to ensure that G_{3} represses G_{4}. This is reflected in Fig. 2D, where a wider range of θ_{43} is obtained for higher values of θ_{35}.

Note that G_{2} (*rho*) is induced next to cells inducing G_{5}, which is downstream of G_{4}, which is in turn repressed by G_{3} (*br*). Thus, G_{2} is not induced in the G_{3}-expressing cells. The formation of the single-cell stripe of *rho* expression is a direct consequence of the action of the incoherent feed-forward loop composed of G_{4}, G_{5} and G_{2}. Initially, G_{2} is activated in the domain overlapping the expression domain of G_{4} (note that G_{2} activation threshold is low). G_{2} is then repressed by G_{4}, but remains in the neighboring cell (where G_{4} is not expressed), owing to the juxtacrine character of the G_{2} activation by G_{5}.

Importantly, the ordered arrangement of cell fates predicted by the model, with a single G_{2}-expressing cell, bounded by the wider areas of cells expressing G_{1} and G_{3} (Fig. 2C,E), persists over a time period that exceeds the decay constants of individual components by at least an order of magnitude. We established that such solutions are robust with respect to small perturbations of model parameters, as long as the remaining seven thresholds satisfy certain inequalities. For example, the activation threshold for G_{3} should be significantly lower than that for G_{1}. Thus, our analysis identified multiple parameter sets that generate the desired arrangement of gene expression domains (Fig. 2C).

The ordered arrangement of gene expression domains in the presented model is generated in a stepwise fashion, by sequentially acting feed-forward loops and juxtacrine signaling. This mechanism may control *rho* expression in the follicle cells. Indeed, similar to G_{2} in our model, *rho* is initially expressed in a wide domain induced by a long-range Grk gradient. Later, *rho* is downregulated in the midline and lateral follicle cells, as well as along the anterior boundary, with expression remaining only in two stripes at the dorsal-anterior borders of *br* expression patches (Peri et al., 1999). Importantly, according to our model, the repression of G_{3} (*br*) by G_{5} (Notch) is not necessary to generate the single cell-wide row of *rho*, but is needed only to abolish *br* from the *rho*-expressing cells.

### Pattern formation in two dimensions

With only a few modifications our model can be used to analyze patterns in two dimensions, on a hexagonal lattice of cells (Fig. 3). First, we modified the spatial distribution of the inductive signal (now denoted by I_{1}, Fig. 3C), making it two dimensional and consistent with the earlier studies of the Grk gradient (Goentoro et al., 2006). We approximate the dorsal-anterior part of the main body follicular epithelium, where the relevant patterning events take place, with the lateral surface of a cylinder, as it is schematically shown in Fig. 3A. The I_{1} distribution is obtained then by solving a steady-state reaction-diffusion equation (Eqn 6, *i*=1.2 corresponds to I_{1} and I_{2}) with zero-flux boundary conditions, with ϕ_{1}=3 and *f*_{1} representing a source function, which equals 1 in a restricted cusp-like dorsal-anterior domain and 0 everywhere else (Zartman et al., 2011). The resulting distribution is shown in Fig. 3B.

In this equation, ϕ_{1} is the ratio of the size of the domain and the characteristic length scale of morphogen degradation.

Second, we added another inductive signal (denoted by I_{2}), which models anterior-to-posterior gradient of signaling by the Dpp pathway (Fig. 3B). This gradient is short ranged and is generated by anteriorly secreted ligand. The distribution is obtained by solving Eqn 6 with ϕ_{2}=14, *f*_{2}=0 and a unit flux along the left (anterior) boundary of the computation domain. Earlier studies established that Dpp represses *br* (Yakoby et al., 2008b). In our model, this corresponds to repression of G_{3} by I_{2} (Fig. 3C). Finally, based on our recent experiments, we restricted the patterning effects of the dorsoventral inductive signal (I_{1}) to the anterior region of the system (Zartman et al., 2011). The nature of this prepatterning effect is still unclear, but it is possible that it results from the earlier phase of EGFR signaling, before the anterior migration of the oocyte nucleus (González-Reyes and St Johnston, 1998; Nilson and Schupbach, 1999).

Adding these effects to the mathematical model requires specifying the threshold of G_{3} repression by I_{2} and the size of the anterior domain within which I_{1} can induce its target genes. Based on the wild-type expression patterns of *br* and *rho*, we selected these parameters in such a way that G_{3} is repressed in the anteriormost two rows of cells by I_{2} and the inductive effect of I_{1} is restricted to the anterior part of the computational domain.

With these modifications, our model was able to provide a quantitative description of wild-type dynamics of expression of *pnt* (G_{1}), *rho* (G_{2}) and *br* (G_{3}) in two dimensions (Fig. 3D,E, where G_{4} and G_{5} are also shown for clarity). The notable feature of the simulated patterns is the formation of two L-shaped single cell-wide patches of *rho* expression (Fig. 3D,E). Thus, parameter sets identified for a one-dimensional model predict the dynamics of two-dimensional pattern formation.

### Analysis of mutants

Given a ‘wild-type’ parameter set, one can readily explore pattern formation in mutant backgrounds. As an example, we provide computational predictions of *br* and *rho* patterns in mutants with altered distribution of Grk (Schupbach, 1987; Price et al., 1989). As was shown earlier, increasing the level of Grk increases the separation between the two patches of *br* (Deng and Bownes, 1997). At the same time, below a crucial value for the amplitude of the Grk gradient, the patches of *br* expression fuse into one (Deng and Bownes, 1997; Lembong et al., 2009). Furthermore, in a mutant with a delocalized Grk, dorsally located patches of *br* expand to the ventral follicle cells, generating an anterior ring of *br* expression (Ward and Berg, 2005). Importantly, in all of these backgrounds, the expression of *rho* is always located at the dorsal-anterior border of the *br* domain and maintains its single cell width (Neuman-Silberberg and Schupbach, 1994).

All of these effects are correctly predicted by our model and are direct consequence of the proposed mechanism in which *br* plays an active role in shaping the pattern of *rho* expression (Fig. 4Aa-d). These are just a few from a long list of predicted patterning defects that can be directly compared with a long list of experimental observations. For example, in the absence of the anterior prepatterning signal (I_{2}) expression of G_{3} expands anteriorly (Fig. 4Ae). This is accompanied by the loss of the anteriormost part of G_{2} expression domain, in agreement with the effect of loss of Dpp on the expression of *br* and *rho* (Yakoby et al., 2008b).

Our computational model can be readily used to predict the results of experiments with genetic mosaics, whereby a component of the network is perturbed in an arbitrarily chosen group of cells. Several examples of computational mosaics are shown in Fig. 4B, where we modeled the effects of inactivating the Notch pathway (motivated by experiments by Ward et al., 2006), ectopic high levels of EGFR signaling (motivated by experiments by Pai et al., 2000), as well as the effects of *pnt* and *br* loss-of-function clones (Zartman et al., 2009; Ward et al., 2006). Thus, we see that loss of Notch signaling in the model leads to loss of *rho* and ectopic expression of *br* in the same cells (Fig. 4Ba). The effect is confined to the anterior part of the Notch clone and leads to ectopic expression of *br* in a single row of cells, consistent with the fact that *br* is also repressead by G1, which is expressed in the midline and independently of Notch. At the same time, a patch of cells with higher level of EGFR signaling can split the dorsal appendage primordium in two smaller domains, with *rho*-expressing cells located at the border of *br*-expressing domains (Fig. 4Bb). This result is consistent with ectopic dorsal appendages that were observed in response to generation of clones of cells hypersensitive to EGFR signaling (Pai et al., 2000). In another example, a midline-located clone of cells lacking *pnt* leads to ectopic *br* expression and an associated *rho* border (Zartman et al., 2009; Boisclair-Lachance et al., 2009). Finally, loss of *br* in the model leads to ectopic *rho* expression, but only in a subset of cells that lack *br*. Within the framework of our model, this reflects the effects of G_{4}, which is de-repressed in the absence of *br*. These computational results demonstrate how the proposed model may be used to summarize the existing experimental observations and to predict spatiotemporal dynamics in a complex regulatory system. Testing these predictions requires future experimental studies of the dynamic interaction of the EGFR and Notch signaling pathways.

## DISCUSSION

Previous studies of *Drosophila* oogenesis characterized the gradient of EGFR activation by Grk, identified many of its transcriptional targets, and began to assemble the regulatory network converting EGFR signaling to complex gene expression patterns (Cheung et al., 2011). A key feature in the emerging eggshell patterning network is an incoherent feed-forward loop. This is the mechanism by which a single-peaked gradient of EGFR activation by Grk gives rise to the two-domain expression pattern of *br*. Cascades of cell-autonomous feed-forward loops can establish patterns of progressively increasing complexity, in which the expression domains of multiple genes are arranged in a specific order. However, cell-autonomous mechanisms cannot readily account for the emergence of fine-grained patterns, such as the pattern of *rho*.

Studies by Celeste Berg and colleagues suggested that, in addition to cell-autonomous mechanisms, *rho* is regulated by juxtacrine signaling, mediated by the Notch pathway (Ward et al., 2006). These studies suggest that the border of the Notch pathway activation represses *br*. Furthermore, it has been shown that loss of *br* leads to ectopic expression of *rho*. Here, we incorporate these observations into a feed-forward model of pattern formation by the Grk gradient and demonstrate that the combined mechanism can indeed account for the wild-type patterns of *br* and *rho* and for changes of these patterns in response to genetic perturbations.

Our model is based on the switch-like description of gene regulation, commonly used in computational studies of gene regulatory networks (e.g. Longabaugh et al., 2005; Perkins et al., 2006). The main parameters in this model are thresholds in gene regulation functions. We demonstrated that most of these thresholds can be estimated based on the qualitative analysis of experimental gene expression patterns, e.g. the size of the *br* domain and its relation to the length scale of the Grk gradient (Goentoro et al., 2006). The remaining parameters can be estimated based on the numerical solution of the model and qualitative comparison of the computationally predicted and experimentally observed patterns. We provided a detailed description of this process, which starts with the analysis of a simple one-dimensional model and ends with the analysis of two-dimensional patterns in the wild-type and mutant backgrounds.

Although our model is based on a number of established components and interactions, some of them remain to be established experimentally. Hypothesized components and interactions may be revealed by studies of the regulatory regions of *br* and *rho*, the outputs of the analyzed network (Cheung et al., 2011). Experiments aimed at the identification of these regulatory regions are currently under way in our group (Fuchs et al., 2012). Future studies of these regions and their control by signaling pathways will provide crucial tests of the model and are likely to shed light on the mechanisms responsible for diversification of *Drosophila* eggshell morphologies (Hinton, 1981; Kambysellis et al., 1995; James and Berg, 2003; Nakamura and Matsuno, 2003; Nakamura et al., 2007; Kagesawa et al., 2008).

## Acknowledgements

We thank Celeste Berg, Jitendra Kanodia, Cyrill Muratov, Laura Nilson, Miriam Osterfield, Trudi Schüpbach, Nir Yakoby and Jeremy Zartman for helpful discussions.

**Funding**

This work was supported by the Human Frontiers Science Program [RGP0052/2009-C].

## References

**Competing interests statement**

The authors declare no competing financial interests.