A Hh-driven gene network controls specification, pattern and size of the Drosophila simple eyes

During development, gradients of intercellular signals (called morphogens) are read and modified dynamically by fields of target cells. As a result, spatiotemporal patterns of gene expression are generated. These patterns are then translated into cell function and into the development of functional body structures (Freeman and Gurdon, 2002; Davidson, 2006). Yet how the integration between intercellular signals and intracellular gene networks occurs is only beginning to be understood.

One of the best characterized family of morphogens is that of hedgehog (hh). Hh genes are evolutionarily conserved and participate in many key developmental processes. Not surprisingly, their malfunction has been associated with a number of developmental diseases and with cancer (Jiang and Hui, 2008; Varjosalo and Taipale, 2008; Ingham et al., 2011). Experimental work in Drosophila and in vertebrates indicates that the Hh signaling pathway is subject to extensive feedback among elements of the pathway, most notably that of the Hh receptor patched (ptc) (Chen and Struhl, 1996; Marigo and Tabin, 1996). Furthermore, the combination of experimental and modeling studies on developing fly wings and on vertebrate neural tube and limbs have uncovered new roles for this feedback in the dynamics of Hh gradient formation, in its robustness and in the generation of distinct patterns of target gene expression (Briscoe et al., 2001; Saha and Schaffer, 2006; Dessaud et al., 2007; Dessaud et al., 2008; González et al., 2008; Nahmad and Stathopoulos, 2009; Dessaud et al., 2010; Irons et al., 2010; Probst et al., 2011; Balaskas et al., 2012). Therefore, the iteration between mathematical modeling and experimentation is emerging as a productive way of illuminating the problem of Hh morphogen action during organ growth and patterning. Here, we investigate this issue in a particularly simple and genetically tractable model organ: the Drosophila ocellar complex.

The Drosophila ocelli are three simple eyes [one anterior (or medial) ocellus and two posterior (or lateral) ocelli] located at the vertices of a triangular patch of cuticle on the dorsal head. Together, the ocelli and the interocellar cuticle (plus its bristles) are referred to as the ‘ocellar complex’ (Fig. 1A). The development of the ocellar complex depends on hh. Flies homozygous for a hh temperature-sensitive mutation raised at the restricted temperature during larval development (Royet and Finkelstein, 1996) or expressing a dominant-negative Ptc receptor [PtcΔloop2 (Briscoe et al., 2001)] lack the ocellar complex (Fig. 1B). Therefore, hh signaling is required for the specification and pattern of two tissue types: ocellus and interocellar cuticle.

The ocellar complex forms by the fusion of the dorsal-anterior domains of the eye discs (Haynie and Bryant, 1986). Here, the ocellar field is specified by the action of, at least, two transcription factors: the pax6 gene twin of eyeless (toy) and the Otx family member orthodenticle (otd) (ocelliless, oc – FlyBase) (Finkelstein et al., 1990; Wieschaus et al., 1992; Royet and Finkelstein, 1995; Punzo et al., 2002; Blanco et al., 2010; Wang et al., 2010; Brockmann et al., 2011). hh is expressed within the ocellar field in the prospective interocellar region (Royet and Finkelstein, 1996; Royet and Finkelstein, 1997). The retinal determination (RD) genes eyes absent (eya) and sine oculis (so) are required for the formation of the ocelli (Bonini et al., 1993; Cheyette et al., 1994; Serikaku and O’Tousa, 1994; Bonini et al., 1998; Blanco et al., 2009; Blanco et al., 2010; Brockmann et al., 2011). Both eya and so are expressed with identical patterns in two domains flanking hh, and mark the prospective ocelli in late third larval (L3) stage discs (Blanco et al., 2009). The expression of eya and so depends on secreted Hh and on their mutual positive feedback (Pauli et al., 2005; Blanco et al., 2009). In addition to the activation of eya and so, expression of the TALE-homeodomain transcription factor homothorax (hth) is concomitantly repressed in the ocellar domains. Otherwise, maintenance of hth expression prevents ocellar development (Brockmann et al., 2011).

In this paper, we have investigated how the single domain of hh expression is capable of generating the ocellar pattern through the regulation of a downstream gene network using both genetic and modeling approaches.

oc2-GAL4 (Blanco et al., 2009) was used to drive UAS lines specifically in the developing dorsal anterior region of the eye-antennal imaginal disc (EAD), where the ocellar region derives from. In the case of UASdsRNAi strains, crosses were raised at 29°C, to maximize the penetrance of the knock-downs. Other crosses were set at 25°C. UAS lines used were as follows: UAS-GFP-ptcΔloop2 [UAS-ptcDN (Briscoe et al., 2001)], UAS-ci (Alexandre et al., 1996), UAS-mamDN (dominant negative) (Kumar and Moses, 2001) and UAS-en (Guillén et al., 1995; Tabata et al., 1995); the hedgehog transcriptional reporter line hh^P30 (referred to herein as hh-Z) and the engrailed transcriptional reporter line en^Xho25 (referred to herein as en-Z), which are from the Bloomington Stock Center (http://flystocks.bio.indiana.edu); and UAS-DlRNAi (28032) and UAS-Su(H)RNAi (103597), which are from the VDRC (http://stockcenter.vdrc.at/control/main). The GFP:Ptc strain [CB02030 from Flytrap (http://flytrap.med.yale.edu/) (Buszczak et al., 2007)] is a GFP-protein trap that tags the Ptc product. Two reference strains, Oregon-R (Or-R) and w¹¹¹⁸, and stocks carrying mutant alleles for patched [ptc: yw; FRT42 ptc^S2/CyO (BL # 6332)], smoothened (smo: w; smo³ FRT40A/CyO), cubitus interruptus (ci: ywhsflp; Ci+ FDT40A /CyO;; Ci⁹⁴) and Notch [N: w[ch2], N^264-39/FM4, B[+] (BL # 730)] are described in FlyBase.

The flip-out method (Basler and Struhl, 1994) was used to induce gain-of-function clones. Clones were induced 48-72 hours after egg laying (AEL) by a 10 heat-shock at 35.5°C in larvae from the cross of yw, hs-flp, act>hsCD2>Gal4;UAS-lacZ females with UAS-hthGFP males (hth+ clones) or yw, hs-flp, act >hsCD2 >Gal4;UAS-GFP females with UAS-en males (en+ clones). en loss-of-function clones were generated through mitotic recombination (Xu and Rubin, 1993) in yw, hs-flp; FRT42D Df(2R)en^E/FRT42D, ubiGFP larvae. Df(2R)enE deletes both the engrailed and invected paralogous genes (described in FlyBase). Clones were induced 48-72 hours AEL by a 45 heat-shock at 37°C. Clones were marked in larval tissues by the absence of GFP. Adult heads from this experiment were mounted and their dorsal head examined for ocellar field defects. Clones were not marked in the adult.

Dorsal head cuticle pieces were dissected from adult or late pharate heads in PBS, and mounted in Hoyers solution:acetic acid (1:1), as described previously (Casares and Mann, 2000). Images were obtained in a Leica DM500B microscope with a Leica DFC490 digital camera and processed with Adobe Photoshop. Ocellar (longest axis) and interocellar lengths were measured with ImageJ (http://imagej.nih.gov/ij/) on digital images and expressed in pixels.

Immunofluorescence was carried out as described previously (Bessa and Casares, 2005). Antibodies used were: guinea pig anti-Hth (Casares and Mann, 2000), guinea pig anti-So (Mutsuddi et al., 2005), rabbit anti β-galactosidase (Cappel), mouse anti-Ptc (Nakano et al., 1989), mouse anti-Eya (10H6), mouse anti-En (4D9), mouse anti-Dl (C594.9B) and rat anti-CiA (2A1) [which detects the activator form of Ci (Aza-Blanc et al., 1997; Méthot and Basler, 1999); all from the Developmental Studies Hybridoma Bank, University of Iowa (http://dshb.biology.uiowa.edu)]. Appropriate Alexa-conjugated secondary antibodies were used. Nuclei were counterstained with DAPI. Image acquisition was carried out in an Apotome Zeiss Axio Imager M2 fluorescence microscope and a Leica SME confocal system. Images were processed with Adobe Photoshop.

Confocal sections of GFP:ptc eye discs at different stages of L3 development, and co-stained with anti-Eya, were selected. To ensure that only signal coming from the prospective ocellar regions were analyzed, each of these regions was outlined with the lasso tool, copied and pasted on a black background and saved as a TIFF file, using Adobe Photoshop. The expression profiles were obtained from these TIFF files using ImageJ.

A temperature fluctuation assay was carried out essentially as previously (Li et al., 2009). Embryos were collected for 24 hours (at 25°C) and grown for an additional 24 hours at 25°C and then transferred to 31°C until larvae reached third instar. Then, larvae were subjected to five cycles of temperature pulses (1.5 hours at 18°C + 1.5 hours at 31°C). After these pulses, cultures were maintained at 25°C until eclosion. As controls, the same strains were grown at constant 25°C throughout development.

A simplified one-cell model (13 equations) was implemented using Vensim software, a visual tool for solving ODEs that allows parameter values modification in run-time (Vensim PLE version 5.11, Ventana Systems, http://www.vensim.com/software.html). This program contains a fourth order Runge Kutta method (RK4) to solve ODE systems. The 31-cell full model was implemented using MATLAB and solved with the integrator ode45.

Each of the 13 differential equations of the reaction-diffusion type describes the behavior of one system variable (gene transcription and protein production) in a row of 31 cells with a symmetrical distribution of cells centered on the morphogen source (five middle cells). Equation 1 describes the classical evolution of a morphogen (Hh) gradient with production and diffusion terms. In this model, the level of complexity was increased by adding a negative regulation, as formation of Ptc/Hh complexes reduces dynamically the concentration of free Hh. The production term is limited to the hh-expressing cells, as expressed in Eqn 14.

Following von Dassow et al. (von Dassow et al., 2000), all other equations distinguish between mRNA transcription and translation. Translation is described using linear terms of production and degradation. Transcriptional regulation is described using non-linear terms, either positive or negative, in the form of compound Hill equations. The specific form of these type of terms is ϕ(Xψ(Y,k2,n2),k1,n1), where ϕ(X,k,n)=Xⁿ/(kⁿ+Xⁿ) and ψ(Y,k,n)=1–Yⁿ/(kⁿ+Yⁿ). The ocellar model contains autoregulations. In these cases, the equation term is described as a simple sigmoid in the form ϕ(X,k,n) (see supplementary material Appendix S1 for further details). For each species, the equation takes specific forms, depending on its specific regulatory relationships (for example, with or without autoregulation term).

The model contains different parameter types: α_x for the basal transcription rates, β_x for the degradation rates, k_x for the Hill equation transcriptional regulators, n_x for the Hill coefficients, θ_x for the translation rates and γ_x for protein complex formation. The non-dimensional parameters k₀, k_Ci, k_En and k_Ciptc are used for changing the scale of different terms and D is the diffusion coefficient. Subscript X-Y, with X and Y system variables, indicates regulation from X to Y. For example, k_En-ptc is the Hill transcriptional regulation parameter of the interaction from En to ptc. All the reaction-diffusion equations contain a degradation term.

In order to understand how the ocellar pattern (Fig. 1C) was generated, we analyzed the expression of ptc (a Hh signaling readout) and eyes absent (eya) (a Hh target) in the prospective ocellar region throughout L3. To monitor Ptc expression, we used a GFP:Ptc protein trap line. In early L3 discs, we detected a single domain of GFP:Ptc expression and uniformly low levels of Eya (Fig. 1D). By mid-L3, levels of Eya rise within the GFP:Ptc-expressing region, also as a single domain (Fig. 1E). From mid to late-L3, the final pattern arises through the repression of Ptc/Eya expression in the prospective interocellar cuticle (Fig. 1F). This pattern suggests the existence of a repressor capable of attenuating the hh pathway in the middle of the ocellar field [see also Brockmann et al. (Brockmann et al., 2011)] and whose expression and/or activity should build up during L3.

engrailed (en) is a candidate hh repressor. It encodes a homeodomain transcription factor with an additional transcriptional repressor domain (Jaynes and O’Farrell, 1991). En is known to repress transcription of two major hh signaling components, ptc and ci, in embryos and wing imaginal discs (Eaton and Kornberg, 1990; Hidalgo and Ingham, 1990; Sanicola et al., 1995; Schwartz et al., 1995; Domínguez et al., 1996; Biehs et al., 2010). In the wing, hh expression in its posterior compartment depends on en. However, hh signaling from the posterior compartment induces en expression in anterior cells at a short range. Therefore, en is a low sensitivity hh target in the anterior wing (Guillén et al., 1995; Ohlmeyer and Kalderon, 1998; Méthot and Basler, 1999). In the ocellar region, hh expression precedes that of en, which is expressed in a hh-like pattern in late L3 (Royet and Finkelstein, 1996). These results suggest that en could be a hh target in the ocellar region. To test this point, we first checked the relative expression of hh, using the hh-Z strain as a hh transcriptional reporter, and of en. Their domains almost completely overlapped, with some En-only cells adjacent to the hh-Z domain (Fig. 2A). Second, when we knocked down the hh signaling pathway in oc2 >GFP-PtcΔloop2 discs, en expression was lost. The expression of so, which follows that of eya, was also lost from the ocellar region, confirming the effectiveness of the knock-down (Fig. 2B). This result indicated that en is a hh pathway target in the ocellar region. To test for en function, we carried out three experiments. First, we verified the status of the hh signaling pathway in the en-expressing cells by examining ci expression. In the en domain, ci is repressed (Fig. 3A), a fact that is consistent with the repressor role of en in other developmental contexts. Second, to test directly this repressing role, we induced marked clones of cells homozygous for the Df(2R)enE. This deficiency removes both en and its paralog invected (inv), thus avoiding potential functional redundancy between both genes. In Df(2R)enE clones spanning the ocellar field the expression of Ptc and Ci is now continuous, lacking the characteristic gap in the prospective interocellar region (Fig. 3B). In adult mosaics, the anterior and posterior ocelli are often fused (Fig. 3C). The fact that the area of the fused ocelli is larger than the sum of the wild-type anterior and posterior ones suggests that the increase in ocellar surface is at the expense of interocellar cuticle. And third, we checked the effects of en overexpression on the hh target eya. In GFP-marked en-expressing clones, eya is repressed in a cell-autonomous manner (Fig. 3D). This eya loss could be explained either by en directly repressing eya or, indirectly, by en blocking the hh pathway. To distinguish between these two possibilities, we overexpressed ci throughout the ocellar field in oc2>ci larvae, therefore making ci transcription insensitive to en regulation. In these larvae, the expression of both eya and en is detected in most of the ocellar field (Fig. 3E,F). Therefore, Ci can activate eya even in the presence of en. In oc2>ci adults, the resulting ocellar complex is composed of a large, single ocellus, without interocellar cuticle (Fig. 3G). This indicates that eya is functionally epistatic over en, and suggests that the primary role of en is as a hh pathway regulator. In all, these results indicate that high concentrations of Hh result in high en expression, which in turn attenuates hh signaling in the middle of the ocellar field. As a consequence, RD expression and ocellar specification can only occur in regions that flank the en domain, which becomes the interocellar domain.

en lays both downstream and upstream of the Hh signaling pathway, being activated by Hh and also repressing the pathway components ci and ptc. Therefore, these genetic relationships should lead to an unstable en expression (indeed, this conjecture was confirmed by our mathematical modeling, see below). Therefore, after its induction by hh signaling, an additional mechanism was required to stably maintain high levels of en expression in a hh-independent manner. It had been reported that in individuals mutant for a Notch temperature-sensitive (N^ts) allele, raised at the restrictive temperature during late larval life, the ocelli fuse (Amin, 2004), generating a ‘cyclopic’ ocellus similar to that observed in Df(2R)enE mosaics. To confirm the involvement of Notch signaling in ocellar development, we genetically manipulated several Notch pathway components. Ocellar-specific knock-downs of the nuclear transducer Su(H) [Su(H)KD] and the Dl ligand (Dl KD), or the overexpression of a dominant-negative form of mastermind (mamDN), a Notch co-activator, resulted in expanded or fused ocelli (Fig. 4A; and not shown). Interestingly, similar knock down of the other Notch ligand, Ser, does not affect ocellar complex development (not shown). The similarity between the Notch pathway mutant phenotypes and the loss of en pointed to Notch signaling being required for en expression. Indeed, in Su(H)KD discs, the en domain is reduced in size and expression intensity and, concomitantly, the two RD domains extend, contacting each other (Fig. 4B,C).

In principle, the input of Dl/Notch in the network could be upstream of hh (maintaining its expression or signaling) or parallel to hh. However, the incomplete activation of en should persist in the former scenario. When we checked the signaling status of the hh pathway in Su(H)KD discs by analyzing ptc expression, we detected an unsplit domain of strong Ptc signal, indicative of sustained hh production and signaling when the activity of the Notch pathway is reduced (Fig. 4D). In addition, the fact that knocking down the Notch pathway still allowed specification of ocelli, which is a hh-controlled fate, agrees with Dl/Notch acting parallel to, or downstream of, hh.

To determine the developmental window in which Notch signaling was required for the maintenance of en, we performed the following experiment. Dl expression was knocked down at different times during third instar, taking advantage of the temperature sensitivity of the GAL4/UAS system (supplementary material Fig. S1; see also Materials and methods), and the size of the interocellar cuticle in adults was analyzed, the fate of which depends on stable and high levels of en expression. Interocellar cuticle surface was estimated by the number of interocellar bristles formed (from 0 in its absence, to 6-8 in the wild type). Disconnecting Dl/Notch signaling using a UAS-Dl-RNAi (oc2>Dl-RNAi or ‘Dl KD’) prior to 80 hours post fertilization (hpf) results in almost total absence of interocellar cuticle. Dl KD during the 80-85 hpf interval results in intermediate phenotypes with incomplete and variable interocellar regions (supplementary material Fig. S1). Interestingly, this developmental window coincides with the establishment of a strong domain of en expression and the split of the eya/so domain in the disc (not shown). Knocking down Dl after 85 hpf no longer precludes the generation of the interocellar cuticle. This result shows that Notch signaling activity is required to establish the interocellar fate during a short developmental interval (coinciding with upregulation of en in the ocellar field), after which, it remains stable. As the interocellar fate depends on en, we interpret this result as en expression becoming fixed by Notch during the 80-85 hour interval.

In order to test whether our genetic reasoning was capable of generating the ocellar pattern, we developed a mathematical model incorporating all known genetic interactions (Fig. 5A). Several simplifications were made. First, the two-dimensional ocellar region is modeled as one dimensional (i.e. Fig. 5B, as a row of 31 cells). Second, the hh transcription domain (the central five cells) is set as a de facto in our model. Third, our model assumes that there is no proliferation during the developmental interval considered (see supplementary material Appendix S1). Hh production and diffusion have been modeled as in Eqn 1 (see Materials and methods), similar to the formalism used by Nahmad and Stathopoulos to model Hh gradient formation in the wing primordium, considering a diffusion coefficient of D=0.5μm²/s (Nahmad and Stathopoulos, 2009). Downstream of the Hh gradient, transcription and translation of all genes have been modeled using ordinary differential equations (ODEs), essentially following the modeling of the Drosophila embryonic segment polarity network by von Dassow and colleagues (von Dassow et al., 2000). Gene transcription may generally be affected by basal (b) and regulated transcription (T), and autoregulation (a), plus a decay term (Fig. 5C). Autoregulation is relevant only for en and eya. Transcriptional regulation terms have been modeled as sigmoids, allowing for potential cooperativity in transcriptional activation and repression. The general form of the transcription and translation equations, as well as the full set of equations are described in Materials and methods. In what follows, we explain how the new regulatory steps have been modeled. Further details on the specific biology underlying other equations (Eqns 2, 3, 4, 5, 6, 7) are described in supplementary material Appendix S1. We have shown that Dl/Notch signaling is required for maintaining high En levels in the interocellar region. We have modeled en maintenance as an autoregulation (Eqn 8), as en has been shown to autoregulate during embryo segmentation (Heemskerk et al., 1991). The contribution of Dl/Notch signaling would be to facilitate the autoregulation of en by lowering k in the autoregulatory term (which indicates the En concentration for which half autoregulatory activation is reached). Because of this, it has been named k_DlEn. This implementation is the simplest form of representing the role of Dl/Notch in allowing en autoregulation we could think of. It considers a constant and uniform Dl/Notch input and that the hh and Notch pathway act independently of one another. The en autoregulation adds on top of a positive Hh signaling input on en transcription (Eqn 8). The expression of eya has been shown to depend only on CiA (Blanco et al., 2009), so no CiR input on eya regulation has been included. In addition, the eya-so positive-feedback loop (Pauli et al., 2005; Brockmann et al., 2011) has been collapsed into a direct eya autoregulation for simplicity (Eqn 10). In addition, previous results had suggested a mutual repression between hth and eya (Brockmann et al., 2011), which is probably direct (supplementary material Fig. S2). Therefore, hth has been modeled as a repressor input on eya (Eqn 10). In addition, hth transcription is modeled as being positively regulated by a constant term (α_wg) (Eqn 12), which represents the likely action of Wnt1/wingless (wg) (Azpiazu and Morata, 2000; Casares and Mann, 2000; Pichaud and Casares, 2000).

The working model has 61 free parameters. For a few, prior biological knowledge is helpful in defining at least some ranges. For example, the basal transcription rates of ptc and ci are positive, as these genes are widely transcribed. In order to generate a working set of parameter values that result in the target ‘wild-type’ pattern (Fig. 5B, seven interocellar and two patches of nine ocellar cells in a 31 row), we first built a one-cell model in which to carry out the first parameter exploration. Then, this parameter set was used as a starting point to manually fine-tune the parameter values on the full model to reach a control pattern (see Materials and methods for further details). With this set of parameter values (supplementary material Table S1), the model accurately recapitulates the target pattern, including the dynamics of Eya, Ptc, CiA, En and Hth expression (Fig. 6; supplementary material Fig. S3). Modeling indicates that before reaching a steady state, the Hh gradient undergoes a transient expansion or ‘overshoot’ (Fig. 6A,B). This early dynamics depends on the non-linear Ptc-mediated feedback (Casali and Struhl, 2004; Nahmad and Stathopoulos, 2009). The model also predicts observed mutant behaviors, including the expansion of the ocellar tissue at the expense of interocellar cuticle in Dl and en loss-of-function mutants, or the effects of hth on ocellar size (supplementary material Fig. S4). Another computational experiment, the overexpression of Dl, predicted the expansion of the interocellar region at the expense of the ocelli. When this prediction was tested experimentally, by overexpressing Dl in the developing ocellar region (oc2>Dl), the interocellar region enlarged and the anterior ocellus disappeared (supplementary material Fig. S4).

An important test for any systems behavior is the stability of its solution and whether this solution is unique or not. To test this point, the initial condition of every system variable was randomized (up to a 10-fold change) in each individual cell (supplementary material Table S2). The solution obtained for the system is stable, as the resulting patterns are the same as the wild type (supplementary material Fig. S5A). Only when the initial condition for En exceeds the concentration determined by the parameter k_DlEn, which is responsible for En autoregulation (k_DlEn>0.2), is en expression fixed throughout, which precludes the establishment of the eya pattern, as expected (supplementary material Fig. 5B,C). In fact, En expression is not detected in the ocellar region until mid-L3, after Eya expression has increased uniformly in the ocellar region (not shown).

Next, and to test whether the network topology is robust to fluctuations, we perturbed all parameters related with production and degradation rates (α_X,θ_X,β_X) with a uniform random signal (white noise). The noise amplitude was 20% for each corresponding rate (Fig. 7A). This fluctuation alters the evolution of the network elements as shown in the Eya time series of Fig. 7B. Under these conditions, the system reproduces the ocellar pattern with a slight deviation (widening) of the interocellar region. This experiment shows that indeed the network model is robust. To further test the robustness of the biological system, we subjected several Drosophila strains to temperature fluctuations during early L3, as a means to increase the noise in the system (Li et al., 2009). We included two reference strains as controls (the wild-type strain Oregon-R and w¹¹¹⁸) and stocks in which the gene dose of smo, ptc, ci and Notch is halved (see Materials and methods), and measured the longer axes of the anterior and posterior ocelli and the interocellar distance. These different genotypes can be thought of as representing the same gene network in which the parameter values may have different, genotype-specific, values. First, we found that different strains showed differences in ocellar and interocellar sizes, indicating that the genotype has a significant influence in the precise size and proportions within the ocellar complex (supplementary material Fig. S6). Second, and directly related to the aim of the experiment, we found that for some genotypes, the temperature fluctuation regime results in size deviations from the control. However, these deviations are smaller than the differences between genotypes. For example, while the difference in ocellar size between w¹¹¹⁸ and Notch–/+ (at 25°C) is about 12%, the temperature fluctuations alter ocellar size in Notch–/+ by only 5%. Furthermore, the external noise introduced did not result in a significantly ‘noisier’ phenotype, measured as the coefficient of variation of ocellar and interocellar sizes (supplementary material Fig. S6). Therefore, the ocellar complex is robust against noise, as the model predicts. In addition, different genotypes, equivalent to the same gene network with a (presumably slightly) different set of parameter values, give rise to ocellar patterns that are also quantitatively different, even if these differences are small. The next section explores the properties of the model across the parameter space.

Next, we carried a parameter sensitivity assay. First, and starting with our working parameter set, we varied one parameter at a time in a two orders of magnitude range around its wild-type value, and measured the departure this variation caused from the wild-type pattern. As a metric for this deviation, we calculated a goodness score as a distance (λ) between the wild-type Eya pattern and that obtained from the varied parameter run, both normalized (see supplementary material Appendix S1). We established three thresholds for the complementary of this distance (1-λ): 1-λ≥0.8 (‘good’), 0.8>1-λ≥0.6 (‘medium’) and 0.6>1-λ≥0.4 (‘bad’), with 1-λ≥0.8 being the most accurate fits (Fig. 7D,E; supplementary material Table S1). Twenty-eight out of 61 parameters gave goodness values above 0.8 over several orders of magnitude and were classified as ‘insensitive’ (i.e. large variations in their value did not result in major deviations from the target pattern) (Fig. 7E; supplementary material Fig. S9.1-9.4; supplementary material Table S1). The remaining ones (33) were sensitive (Fig. 7D), although the range of values for each specific parameter that resulted in a score of at least 0.8 varied among them (see supplementary material Table S1). The most sensitive parameter is the basal transcription rate constant α_en, which is also the parameter that does not allow variation in the noise experiments. This suggests that en expression has to be kept strictly off in the absence of patterned Hh signal. Not surprisingly, other sensitive parameters are those related to ptc and ci expression, which affect the major feedbacks in the network (supplementary material Table S1).

Second, we carried out an analysis in which all 33 sensitive parameters were simultaneously randomized at each run inside one of the goodness intervals (the remaining parameters were left fixed at their wild-type values). A total of 10,000 runs were obtained, distributed among ‘good’ (6000), ‘medium’ (3000) and ‘bad’ (1000) randomized parameter values. For all 10,000 parameter sets, the distance for all the patterns of the system (one per variable) to the wild type was calculated. In this way, each parameter set defines a point in a 13-dimensional space, each dimension being one of the model variables (see supplementary material Appendix S1). For visualization, this information was reduced to three dimensions, two of them being the projections of the distance distributions two by two, and another representing the density of multidimensional points in such projections. The resulting 2D histogram (Fig. 7F; supplementary material Fig. S7) plots the density of patterns from the 10,000 randomized runs distributed relative to their distance from the wild type (point 0,0) (see supplementary material Appendix S1 for details). Therefore, it represents a map of the phenotypic space generated by the GRN (gene regulatory network) using randomized sets of parameters. Several conclusions can be derived from this analysis. First, the distribution of solutions (patterns) was not evenly dispersed. Instead, the solutions tended to concentrate in clusters or ‘islands’. Second, the wild-type pattern was placed inside a big and dense cluster, so this solution is highly probable, which indicates that the pattern is stable. Finally, when we analyzed the patterns of Eya and En (the two major readouts of the GRN) located in dense islands far apart from the wild type, we found that those patterns were still qualitatively similar to the wild type [Fig. 7F; see, for example, the high stability island around (–0.7,0.7) in supplementary material Fig. S7B]. This is interesting, because in this experiment we used parameter values coming not only from the ‘good’ interval, but also from ‘medium’ and ‘bad’ ones, as derived from our previous single parameter analysis. In summary, these analyses indicate that the GRN is robust, because wide and random parameter variation results in specific phenotypic clusters, all of them similar to the wild type.

Many gene networks are described as static regulatory (activating and inhibitory) relationships between network components (genes and their products), disregarding essential dynamic and spatial aspects. In such descriptions, all genetic interactions are given as though happening at once and without spatial context. Mathematical modeling allows us to test the logic consistency of such networks, and whether or not they are capable of explaining the spatial and temporal dynamics of the biological system.

Using an integration of experimentation and mathematical network modeling, our results help to explain how several alternative fates are controlled by the Hh morphogen (Fig. 8). Previous descriptions of the genetic interactions involved in the specification and patterning of the ocellar complex structures did not offer satisfactory explanations for this fate choice decision.

A first important point is the addition to the GRN of en as a hh target with self-maintaining capability. The transduction of the Hh gradient generates an initial asymmetry, with only the cells receiving the highest Hh concentrations being able to maintain en expression. This in turn sets in motion the dynamics of the GRN. The evolution of some key system components is shown in Fig. 8.

A second important point is the action of en as a Hh-pathway repressor. The fact that en expression is sustained just in cells receiving the highest Hh concentrations (the Hh-producing cells and their adjacent neighbors) makes these cells read the Hh signal only transiently, as the signaling pathway is blocked as en expression builds up. This means that it would be impossible for en to reach sufficient expression levels to shut off the pathway – and therefore, to inactivate eya – unless additional mechanisms were considered. In fact, the inactivation of eya is necessary for the specification of the interocellar region: thus, uniform and high ci expression results in the co-expression of eya and en throughout the ocellar field. In this situation, eya is functionally epistatic over en, and the only tissue type specified is ocellus. Therefore, a stable interocellar region can be established only if the initiation of en expression is followed by a hh signaling-independent phase. Such a transition from signal-induced expression to independent mode of maintenance has been reported for en during Drosophila embryonic segmentation (Heemskerk et al., 1991). In the ocellar field, we propose that this transition requires Notch signaling, specifically activated by its ligand Dl (but not by Ser). The molecular mechanisms of this en maintenance are not yet clear, but might involve PREs (polycomb response elements) in the en locus (Kwon et al., 2009).

Our model includes another repressor, hth, which enters the network as a direct repressor of RD. Its contribution seems limited to restricting the external extent of the eya/so expression domain. In hthKD animals, the ocelli are larger, but the interocellar cuticle is still present. In our model, en repressive action suffices to turn off hh signaling pathway, thereby precluding RD activation. We have tested, through modeling, the possibility of hth being required for en activity, as it has been shown to be the case during embryonic segmentation (Kobayashi et al., 2003). In this case, though, making the repression function of en dependent on hth does not allow the network to reach any steady state in which the interocellar domain is established – i.e. en does not reach the maintenance threshold. To verify this prediction, we checked en expression and the activity status of the hh pathway in hthKD discs. As predicted, en is expressed at normal levels in a domain where ci is off, as in wild type (supplementary material Fig. S8).

Recently, a similar situation to the one we detailed here has been described during the dorsoventral patterning of the vertebrate neural tube by Shh (Ribes et al., 2010). The floor plate (FP), the ventral-most region of the neural tube, expresses Shh and requires maximal Shh concentrations for its specification (Chiang et al., 1996; Ericson et al., 1996). However, this requirement is transient and followed by an attenuation of the pathway. This attenuation is necessary for FP specification (Ribes et al., 2010). This process of FP specification is reminiscent of the specification of the interocellar cuticle in our system. This similarity raises the possibility that a negative-feedback loop in the Shh pathway, of the type we have described here, could be part of the neural tube GRN. However, after the initial asymmetry within the ocellar field has been established, an external input, the Dl/Notch signal, is needed to maintain it. In our model, there is no need for a localized Notch signal: uniform signaling suffices, provided that en expression reaches a specific concentration threshold. In fact, using an anti-Dl monoclonal antibody, we detect uniform levels of Dl expression in the ocellar field in mid-L3 (not shown), the developmental period when we start to see the distinct domains emerging. Interestingly, a recent report finds an association between mutations in Dll1, a vertebrate Dl-like ligand, and holoprosencephaly (Dupé et al., 2011). Holoprosencephaly, a developmental defect caused by abnormal specification of the ventral midline structures of the anterior neural tube, is frequently associated with malfunction of the Shh pathway. In fact, work in vertebrates indicates that Notch signaling is indeed required for FP fate acquisition parallel to Shh (le Roux et al., 2003; Peyrot et al., 2011). It is therefore tempting to speculate that the Notch pathway is required to fix the FP fate in the vertebrate neural tube by stabilizing gene expression during the phase of Shh signaling attenuation, as we propose here for Notch in the Drosophila ocellar GRN.

The model we have described explains how mutually exclusive gene expression domains are produced under the control of a Hh gradient by connecting a repressive gene circuitry. These expression domains underlie the morphology of the structures that will later form and cannot be explained unless the information of Hh gradient is integrated with the logic of the GRN. They can therefore be considered as the emerging properties of the system we have described.

The structure of this GRN confers robustness to the patterning mechanism, buffering variations in the initial conditions, as well as absorbing noise. Although the model predicts an early overshoot of the Hh gradient, this might not be crucial, as variations in the initial conditions converge to the same pattern.

Interestingly, random variation of parameter values results in the system deviating from the ‘wild-type’ pattern in a non-random manner, but instead falling into specific ‘islands’ of the phenotypic space. That is, variations in the control parameters of the GRN generate phenotypes that maintain certain rules of proportionality. Still these phenotypic ‘variants’ are robust against noise, as is the wild-type pattern (supplementary material Fig. S7). These mathematical properties of the ocellar network might ensure the phenotypic stability of the ocellar structures in wild flies exposed to varying environmental conditions during development, as well as constraining the phenotypic variability of the ocelli during evolution.

We thank I. Guerrero, I. Rebay and R. S. Mann for antibodies; P. Rojas, J. Blanco, C. M. Luque, the Bloomington, VDRC and Flytrap Stock Centers, and the Consolider ‘From Genes to Shape’ Drosophila Collection for fly strains; and O. Abreu for help with experiments.