Core regulatory network motif underlies the ocellar complex patterning in Drosophila melanogaster
Abstract
During organogenesis, developmental programs governed by Gene Regulatory Networks (GRN) define the functionality, size and shape of the different constituents of living organisms. Robustness, thus, is an essential characteristic that GRNs need to fulfill in order to maintain viability and reproducibility in a species. In the present work we analyze the robustness of the patterning for the ocellar complex formation in Drosophila melanogaster fly. We have systematically pruned the GRN that drives the development of this visual system to obtain the minimum pathway able to satisfy this pattern. We found that the mechanism underlying the patterning obeys to the dynamics of a 3nodes network motif with a double negative feedback loop fed by a morphogenetic gradient that triggers the inhibition in a French flag problem fashion. A Boolean modeling of the GRN confirms robustness in the patterning mechanism showing the same result for different network complexity levels. Interestingly, the network provides a steady state solution in the interocellar part of the patterning and an oscillatory regime in the ocelli. This theoretical result predicts that the ocellar pattern may underlie oscillatory dynamics in its genetic regulation.
keywords:
Pattern formation, ocellar development, Reactiondiffusion model, Boolean networks1 Introduction
The study of biological systems from a physical perspective has substantially increased in the last decade as a result of the same increment in the synergy between theoretical and experimental researchers in this field. Physical models help on a better understanding of the principles of biological systems. One of the biological areas that counts with a higher number of theoretical collaborations is developmental biology. This biology field, in the case of multicellular organisms, studies how these organisms develop from a single cell to a complex structure with many differentiated cell types capable of performing distinct functions. Here, the time parameter gives an extra dimension able to show complex emergent properties, and this is indeed very attractive for physicists interested in the study of complex systems.
Evolutionarily speaking, every developmental process is subjected to mutations that, in some cases, may lead to a modification in the organism that can imply a certain improvement. In many other cases, mutations can lead to a failure and stop of the developmental process. In order to avoid this type of failure, the biological system needs to be as robust as possible to perturbations in its developmental program. Organisms containing a high level of complexity are able to maintain certain modifications between individuals while maintaining robustness in the essential processes. Low complex organisms usually present a clear robustness in all the developmental processes. A noteworthy example is that of the Caenorhabditis elegans worm. This worm has exactly 959 somatic cells in adult hermaphrodites and 1031 somatic cells in adult males Sulston and Horvitz (1977), which is an impressive sign of robustness in its growth process.
In the present work we are focusing in another organism, Drosophila melanogaster (D. mel.), also known as the fruit fly. In D. mel. our interest relies on the robustness of a specific tissue patterning, the ocellar complex. The ocellar complex is part of the visual system of many insects. In this fly, it is formed by three ocelli (or simple eyes), situated at the vertices of a triangular patch of cuticle, and their interocellar region, placed in the dorsal head. This singular patterning is actually formed by the fusion of the dorsalanterior domains of the two eye discs Dessaud et al. (2010), so it can indeed be studied as the two ocelli pattern of one of the two eye imaginal discs. A onedimensional description is then proposed where two ocelli regions are separated by an interocellar cuticle. A physical model based on reactiondiffusion equations has already been proposed in AguilarHidalgo et al. (2013) to fulfill, not only the tissue patterning for the Drosophila ocellar complex, but also the correct experimentally tested behavior of the gene regulatory network (GRN) that drives its formation.
This specific pattern follows the French flag problem paradigm, where different developmental programs are run depending on the concentration of an extracellular diffusive molecule (morphogen) Kondo and Miura (2010).
Following Kang et al. (2012), in order to test robustness in a physical model we can identify three classes of perturbations that lead to three types of robustness:

Robustness with respect to changes in the model parameters.

Robustness with respect to transient changes.

Robustness with respect to changes in the structure of the equation itself.
Types 1 and 2 have been tested in AguilarHidalgo et al. (2013). The equations there described provide the same pattern when the system is subject to noise perturbation and changes in the initial conditions. Moreover, an extensive parameter sensitivity analysis was performed, where perturbations in the parameter values retrieved the same patterning with changes in size of its components Kang et al. (2012), what suggests that the same GRN determines a certain phenotypic variation.
The third type of perturbations considers structural changes in the equations describing a physical model. In a networked system, an exercise to test this type of robustness is to modify the wiring of the network by adding/removing links like in Li et al. (2004). This, in a GRN modifies/removes regulatory interactions. As mentioned above, many regulatory interactions related to the ocellar complex developmental GRN were tested in AguilarHidalgo et al. (2013), and agreed with the model. However, there is an important issue related to structural robustness of the system remaining untested. This is the model description level Schlitt and Brazma (2005). As the patterning retrieved by this model is robust facing parameter values modification, noise, changes in the initial conditions and mutations, we are interested to know the reason for such robustness in the patterning retrieval.
In the present work we focus in perturbations on the structure of the equations with respect to the description level of the physical model. The model presented in AguilarHidalgo et al. (2013) describes the behavior of the GRN for the development of the ocellar complex in a high level of detail, including both genes and products representations, complex ligandreceptor dynamics and a nonlinear formulation. Here, we have reduced the description level of the model while maintaining the biological sense according to the results in AguilarHidalgo et al. (2013), and thus, we have tried to find out the minimum network expression capable of showing the correct pattern. To do this, we have systematically removed nodes and links and, then, changed the formulation from a nonlinear to a linear paradigm.
The patterning seems at this point to be independent on the description level and formulation complexity. Both the patterning and the individual behaviors of each GRN component is correctly fulfilled. An observation on the resulting GRN leads to the extraction from this GRN of a core formed by a threenodes network motif that satisfies the correct patterning by itself. This core network is biologically very relevant as it defines the main morphogensignaling pathway in the ocellar development with the addition of its main network repressor in a double negative feedback loop. A further analysis shows that even the reduction of this core network to a simple activatorrepressor feedback loop retrieves the desired patterning. All the different models described are then simulated using Boolean networks. This changes the description level of the dynamical system to its simplest representation Schlitt and Brazma (2005). This exercise confirms that this patterning is independent on the description level but dependent on specific regulatory topology of the activatorrepressor motif.
2 Results
2.1 Description of the system: the ocellar complex pattern
During fly development, the ocellar complex is formed describing a triangular patch with a simple eye on each vertex: one anterior ocellus (AO) and two posterior ocelli (PO) (Fig. 1a). It must be clarified that the fly develops two eye discs, and each eye disc forms two ocelli, one PO and one AO. In a late larval stage the two anterior ocelli are fused into one and only anterior ocellus. In order to describe the ocellar complex conformation in a simplified way, it can be reduced to a onedimensional description by analyzing one row of cells situated in a cross line over the two ocelli of one of the two eye discs. This conformation shows a 1D field compound of five blocks of three different types of tissues. The external tissue will be referred hereinafter as periocellar region (POC). The following inner type of tissue corresponds to the region where the ocelli will be placed; this tissue will be referred as ocellus (OC). The tissue between the two ocelli will be named interocellar cuticle (IOC). The distribution of these types of tissue in the 1D field will be POCOCIOCOCPOC (Fig. 1b).
As mentioned above, the ocellar complex pattern in D. mel. is dependent on a morphogen profile. One of the most studied morphogens is the molecule Hedgehog (Hh), which has already been considered to successfully describe intercellular communication in developmental processes Nahmad and Stathopoulos (2009). This morphogen is known to be essential for the formation of the ocellar complex in D. mel. Royet and Finkelstein (1996). Here, Hh is expressed in a central stripe of cells in the target tissue diffusing in the anteriorposterior axis. The readout of this system is the retinal determination (RD) genes group. These 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, 2010); Brockmann et al. (2010). 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). The gene engrailed (en) is activated by the Hh signaling pathway and it selfregulates only in the cells where its protein (En) concentration is high. This En high concentration region overlaps with the Hh expression cells though the En domain can vary depending on external conditions and variations in the GRN AguilarHidalgo et al. (2013). En, then, attenuates the Hh signaling pathway to establish a region with no RD expression (IOC) and thus, the ocellar pattern.
An extensive experimental and theoretical work in the ocellar complex shows robustness in the patterning against noise, changes in the initial conditions and parameter values AguilarHidalgo et al. (2013). There, some structural modifications were also tested by removing links in the regulatory network. As example, the ocellar pattern is maintained with enlarged ocelli if the repression between Homothorax (Hth) and is removed; on the contrary, it is noteworthy to mention an example by which the ocellar pattern is not fulfilled, this takes place when removing the link between Delta (Dl) and En; here the anterior and posterior ocelli fuse and just one big ocellus is developed. These structural modifications were experimentally validated and are also considered in the present work. The proposed model in AguilarHidalgo et al. (2013) for the GRN in (Fig. 1c) is implemented using a nonlinear formulation and is given by the following equations:
(1) 
(2) 
(3) 
(4) 
(5) 
(6) 
(7) 
(8) 
(9) 
(10) 
(11) 
(12) 
(13) 
(14) 
where and , with X any system variable. The model contains different parameter types: for the basal transcription rates, for the degradation rates, for the Hill equation transcriptional regulators, for the Hill coefficients, for the translation rates; for protein complex formation. The nondimensional parameters , , and are used for changing the scale of different terms and is the diffusion coefficient. Subscript XY, with X and Y system variables, indicates a regulation from X to Y.
2.2 Robustness in the patterning circuit against changes in the structure of the system
Here we propose a series of modifications in the structure of the system developed in AguilarHidalgo et al. (2013) in order to test robustness in a different way as done there and to try to infer the basic mechanism responsible of the patterning. These changes include pruning of the GRN by eliminating certain nodes and links. These modifications assume changes in the behavioral paradigm of some regulatory elements that, while maintaining the overall concept, they simplify the corresponding mechanisms. Additionally, changes in the formulation paradigms are also being considered.
2.2.1 Reduction of the GRN leads to proper patterns
In order to find the basic mechanism to satisfy the ocellar pattern it is possible to systematically reduce the GRN in AguilarHidalgo et al. (2013) (Fig. 1c) by removing nodes and links as the model is simplified and the system of equations is reduced.
To simplify this system we firstly focus on the receptor dynamics. This previous model implemented a detailed and complex receptor dynamics involving some downstream targets. The receptor production was modelled as the gene ptc expression with two different kinetics. The first one is a basal constant production and the second one is an input from CiA restricted by CiR, both nonlinearly represented. These two terms were repressed by a nonlinear En contribution. To attempt this simplification it is proposed a constant number of receptors Umulis (2009) in each cell with a certain binding () and unbinding () rates. This way eq. (2) for ptc can be eliminated, the translational rate in eq. (3) and eqs. (1, 3 and 4) can be written as follows:
(15) 
(16) 
(17) 
The result of this simplification is a replication of the ocellar pattern, as given in the previous model, for every species involved except for Ptc and PtcHh due to the absence of En interaction. This modification reveals the need to preserve the interaction between En and the receptor dynamics as, in order to maintain the full biological sense, PtcHh needs to show the correct patterning. Note that this will be regulated in a subsequent modification to the model by reinserting En repression and CiA feedback.
The next downstream step in the network is the activation of ci dependent proteins. In the former model, the gene is modeled to be expressed dependent on a nonlinear En repression, eq. (5). Then, CiA is produced with a linear dependence on expression, eq. (6), which is turned into CiR form inversely dependent on the ratio of bound and unbound Ptc to Hh, eqs. (67). In a coarsegrain view of this process, the activation of CiA correlates with Hh signalling inside the cell, and so, with the amount of Hh bound to its receptor. To perform this simplification it is possible to go through different small steps which all maintain the desired pattern. Firstly, eq. (5) can turn from a constant production repressed by En to a Hhlike profile repressed by En when inserting a contribution:
(18) 
This means that the ci expression is dependent on both the concentration of receptors bound to ligands and free receptors, which is a plausible result to be experimentally tested. Note that this interaction in the real system could be done either direct or indirect. In the latter case, an unknown species downstream the receptor dynamics could sense both, the concentration of free and bound receptors and, accordingly use this information to activate ci expression.
Now, as ’s output is already in a Hhlike profile in eq. (18), one can assume that CiA is selfstable and thus neglect the transformation from CiA to CiR in eqs. (67), and so, eq. (7) for CiR can be eliminated. CiA dynamics is shown in eq. (19).
(19) 
At this point, the system stands for another change in a coarsegrain view. One can avoid the description of the receptorligand binding process if the downstream targets maintain their profile, as done in Zhou et al. (2012). As it has already been commented, spatial form correlates with PtcHh and Hh profiles. Therefore, the next proposed simplification assumes that the bound Hh can behave as an internalized molecule dependent on Hh concentration instead of following the binding process. This, in the real system, can stand for the description of a molecule situated downstream of PtcHh and upstream of in the Hhsignalling pathway. This other molecule senses the effect of Hh and can be altered by other elements as En and CiA. To simplify notation, PtcHh is used as this internalized molecule, now sensitive to En repression. This analysis resumes the need of including En interaction to the PtcHh species commented above. With this assumption, the first part of the pathway is modeled as follows:
(20) 
(21) 
Note that an equation for unbound Ptc can be eliminated from the system as the binding process is no longer described. Most of the constant basal productions were set to zero in the original model. It has also been tested that the system keeps the same behavior when neglecting them so they are removed from the equations to simplify the equations structure at this point of the network pruning.
Now a straight forward simplification is performed to the system. The original model showed a detailed dynamics involving both genes expression and protein dynamics. One can notice that in the steady state, the genes expression and the protein concentration differ in a proportionality constant. According to this, one can represent the system using just one of the element types adjusting this proportionality constant in each case. The system is then shown in eqs. (2025).
(22) 
(23) 
(24) 
(25) 
The proposed GRN to fulfill the ocellar complex pattern is shown in (Fig. 1d). A sum up description of this network is the following. The morphogen Hh diffuses extracellularly and forms the input profile of the system. Hh enhances an intracellular PtcHh concentration which causes a decrease in Hh concentration. Then, PtcHh results in signaling for Cubitus Interruptus (CiA) in its active form. CiA activates the inhibitor Engrailed (En) that applies negative feedback loops at two different points on the Hh signaling pathway, PtcHh and CiA. Finally CiA activates the formation of Eya which interacts with Homothorax (Hth) in a switchlike form, so Hth is present where Eya is not, and vice versa. Hth is activated by Wingless (Wg). Though Wg is a well known morphogen, in these models it is considered as a constant input. Wg is, then, not behaving as a graded morphogen but a constant input. This takes into account that its effect in the system would be qualitatively similar if its expression is flat in the whole imaginal disc Alexandre et al. (2014). After taking the simplifications where many of the interaction mechanisms are referred in a coarsegrain description, we will refer hereinafter to Eya as RD proteins, without particularly considering any of them but showing a general overview of the RD complex dynamics. Then, the spatial domains in the model where RD is present correspond to the ocellar cuticle. (See complete results in Supplementary Material).
2.2.2 Changes in formulation paradigm maintain the correct patterning
In order to perturb the structure of the model’s equations, the system has been substantially pruned from that in AguilarHidalgo et al. (2013). Now one can question whether the patterning can be affected by the formulation used. In order to overcome this issue, the model has been rewritten using a mostly linear formulation (eqs. (26) to (33)). This way, Hill’s function formulation or any other kind of relaxation kinetics used to simulate more realistic biological behaviors are avoided except in necessary cases (see Supplementary Material).
(26) 
(27) 
(28) 
(29) 
(30) 
(31) 
(32) 
(33) 
The coefficients drive an action (activation or repression, depending on sign, or respectively) from source element X to target element Y; is the degradation coefficient of the element X; is the maximum Hill term value when X regulates Y.
The proposed linear formulation is not conservative so some parameter sets could lead to negative concentrations. To avoid this issue Heaviside functions are introduced. In this new system, equations (2732) are the reactiondiffusion type differential equations of the model. The description of the system is the same as in the last section. Some nonlinearities are applied in order to maintain the correct pattern: Firstly, En repression of CiA is shown as a reduction of the maximum level that CiA can reach (). A stability study shows that this kind of nonlinear repression makes this system have a stable point, while linear repression leaves the system without any stable solution (see Supplementary Material). Secondly, in order to achieve the desired pattern, En needs to be maintained in a Hh signalingindependent manner. As CiA is the only En activator, En should selfmaintain. This new behavior should not work in every cell, as the result would be a complete En domain with no response to Hh influence (Fig. 2a); so the En selfmaintenance should work only under certain circumstances, and these are Hh signalingdependent. In order to treat this feature in a simplistic way En is selfmaintained if its concentration exceeds a certain threshold, , for a time interval, , starting from the time point when En reaches the concentration threshold, , (eqs. 42 and 33). This requirement refers to the need of integrating a protein quantity over time to trigger some targets Dessaud et al. (2010).
These premises form an En domain in the IOC allowing the production of the RD proteins in the OC region (Fig. 2b). If En transcription is eliminated () or substantially reduced so that it cannot reach enough concentration to selfregulate, then RD appears in a single unsplit domain that follows the Hh profile. There is experimental evidence showing that at an early developmental stage (mid L3) RD is appears as a single unsplit domain, similar to the model’s output. However, this pattern evolves so that by the end of L3, RD expression and Hh signaling are detected in two domains. These domains are separated by a region in which RD and Ptc expressions are very low or absent Dessaud et al. (2010); AguilarHidalgo et al. (2013). This pattern suggested the existence of a repressor capable of shutting down the hh pathway. The same pattern is observed if Dl signaling were shut off (), the mathematical effect would be identical, En selfregulation term will not cause any effect and RD will appear in only one domain over the whole tissue.
The pruned GRNs both linearly and nonlinearly formulated are able to reach the correct pattern formed by two OC domains separated by the IOC (Fig. 3). In this interocellar region Eya/RD values are null or very low while En concentration is high; likewise, in the ocellar domains En concentration is null or very low and RD concentration is high. The regulation of the different domains depends on the Hh gradient concentration value on each spatial point. En concentration is dependent on CiA concentration which also depends on Hh concentration as output of its signaling pathway. Therefore, the spatial point where En is starting selfregulating depends on Hh concentration. This behavior responds to the classical French flag model. It is also noteworthy that Hth is only present where RD is not, due to the double negative feedback control between these two elements. So Hth shows high concentration in the interocellar and periocellar regions. PtcHh and CiA behaviors are also shown as expected. They reflect the Hh gradient as being part of the Hhsignaling pathway and feel the negative effect of En repression in its domain (IOC). One must note that the IOC domain is set in a more abrupt way in the linearly formulated model than in the nonlinearly done one. This refers to the use of a step function for the En selfregulation (linear model) instead of a Hill’s equation (nonlinear model). As in AguilarHidalgo et al. (2013), the simplified model is robust against noise perturbations. This model also shows a size control of the OC when perturbing the activation of Hth by Wg (see the details of these results in Supplementary Material).
2.3 A 3nodes network motif retrieves the patterning: The core
The study made to this point shows that the correct patterning seems to be extremely correlated to the use of a negative feedback loop with a dependence on the morphogen concentration to trigger this inhibition. To finally prove this, a single 3node network motif has been extracted maintaining the negative feedback loop. The description of this motif is the following: PtcHh activates CiA, CiA activates both En and PtcHh, and En represses the previous two network nodes in a double negative feedback loop (Fig. 4a). The input to this network motif is, in an adiabatic approximation, a stationary Hh gradientlike profile to simplify the calculations, avoiding this way the PDE equation for Hh. Both the analytical solution for equation (27) in its stationary state or its Gaussian fit (see Supplementary Material.), provide the same input to the system as in the dynamical case in a steady state situation. Fig. 4b shows the comparison between the numerical stationary solution of the Hh morphogen gradient and the Gaussian distribution fit (see details in Supplementary Material). Then, the same system’s behavior is maintained. By doing this, the system becomes a system of ODEs of the reaction type.
The resulting system is formed then by three equations.
(34) 
(35) 
(36) 
These equations (3436) have the same formal expression as in eqs. (28) (PtcHh), (41) (CiA) and (42) (En) but considering a constant Hh profile in PtcHh equation (Fig. 4b).
This 3nodes ocellar model is still capable of showing qualitatively identical expressions for PtcHh and En (Figs. 4c and 4d respectively) and the correct ocelli pattern in CiA element (Fig. 4e) without changing the parametric values from the former model. It is noteworthy to emphasize that both a dynamic morphogen profile and its stationary solution retrieves the correct pattern. Thus, one can provide two important results. The first one implies the use of the dynamical description. Here, the behavior of the system is positively described at different developmental stages when generating the Hh profile from null concentration values (as seen in previous sections). The second result comes from the solution of the patterning when using an adiabatic approximation to the morphogen dynamics. Here, one can obtain the correct pattern when using a steady graded profile as the system’s input. In other words, the ocellar patterning can be fulfilled also in case the morphogen kinetics is much faster than the intracellular response.
2.4 Reduction of the system to an activatorinhibitor dynamics
The 3nodes network model raises the question of whether an ultimate simplification on the system will lead to a correct patterning. In this case, though biologically unrealistic, it is considered that Hh is directly feeding CiA production with no intermediate; CiA is then activating En and this last represses CiA formation. The correct pattern is still maintained. It is clear that this French flag type model forms the ocellar pattern due to this activatorinhibitor interaction. The patterning is correctly achieved if this specific set is maintained in the core of the network.
Following the 3nodes model formulation it is possible to write this system in terms of ODEs as follows:
(37) 
(38) 
This simple system allows an analytical treatment to find the steady state solutions for En and CiA. In order to ease the analysis, the steady state solutions for En and CiA in the steady state are written in a dependent manner:
(39) 
(40) 
A simple observation of eq. (40) shows that in the case of En selfregulation, this is , if the parameters and are equal or very similar, En concentration goes to infinity or to a very high value; and this leads to a null or very low CiA value. So the Hh signaling pathway is closed and no ocellar cuticle is developed. In the case of , depending on the parameter values, positive or imaginary values are obtain for both En and CiA. Though imaginary solutions can be found if , only real solutions are considered, and these are positive. So in this case, CiA expression is maintained and OC is formed.
2.5 Boolean networks analysis confirms robustness in the patterning circuit
It has been demonstrated that several formulations of the same GRN with different levels of complexity are able to determine the same patterning. This pattern seems indeed to be robust to the choice of model if the activatorrepressor module is maintained in the network. In this section, in order to test the robustness of this pattern independently of the parameter values and the dynamics of a meanfield formulation, the GRN has been studied from the perspective of Boolean networks. For this, four models have been tested, the one in Fig. 1c for the model in AguilarHidalgo et al. (2013) (model ), the one in Fig. 1d (model ), the described in Fig. 4a (model ) and the activatorrepressor motif network detailed in section 2.4 (model ).
To perform the Boolean models, the state of each node is 1 or 0, depending on the presence or not of the corresponding genetic element. The states of the nodes can change in time, and the next state of node is determined by a Boolean (logical) function of its state and the states of those nodes that have incoming edges on it. In general, a Boolean or logical function is written as a statement acting on the inputs using the logical operators and, or and not and its output is if the statement is true (false). For the construction of the functions it is assumed that the inhibitors are dominant in the networks, so they will be treated as AND logical operators Albert and Othmer (2003). This type of network modeling has already been widely used to study genetic networks Kauffman (1969); Thomas (1973); Thomas and D’Ari (1990); Kauffman (1993); Albert and Othmer (2003); González et al. (2006); Buceta et al. (2007). Table 1 gives the Boolean functions for every genetic element of the four considered models.
Models , , and show the same behavior. In case that the network is describing the behavior of a cell in the IOC region (Fig. 5a), Dl = 1 so En selfregulates, the Hhsignaling pathway is blocked and the RD genes are not expressed. This situation is given as a steady state regime in the four analyzed models. It is found the steady state takes longer to be reached as both the number of elements involved in the system and their relationships increase. In model , the active nodes in the steady state are Hh, , En, , Wg, Hth and Dl. Model shows the same behavior so the only active nodes are Hh, En, Hth and Dl. In models and the active nodes are then Hh, En and Dl. The IOC is well defined in the four models. In the case of the OC region (Fig. 5b), Dl is set to 0 and En does not selfregulate. Interestingly, this state in the four models also show the same behavior though reveals an oscillatory regime. One can observe that CiA and RD (depending on the studied model) are active in time, so the ocellar pattern is fulfilled though in a nonsteady situation. Models and show a transient, that is longer in the first one. Models , and show the same period length of 5 time steps while model present a shorter period of 4 time steps, this is due to the reduction of the number of nodes between the input, Hh, and the repressor feedback.
Models and show a dependence on the initial conditions to express a different behavior in nodes related to the outcome of the network. These are , Eya, and Hth in model , and RD and Hth in model . This dependence is due to the mutual double repression between the Eya/RD and Hth species. The switchlike mechanism will activate only one of the two family species or both in an oscillatory fashion depending on the initial conditions given. In Fig. 5 is shown the behavior where the only active node is Hh as initial condition for models and . For model , Wg is also considered active from the initial step. In the OC description of this model, has also been set active as initial condition in model #1 to allow the Eya and Hth species rule in an oscillatory manner. If is set inactive as initial condition, the Eya species are set to 0 as Wg is always 1 and it activates Hth before Eya can be turned on. Therefore, the switch is kept this way in time. In model #2, in order to maintain the switch always activating the RD species it is necessary to initialize RD as 1. If this is not considered, Hth will finally block RD in a persistent way as the oscillation period of CiA is larger than the number of links between RDHth in this model, so even initializing Hth will block RD if RD is 0 in its initial state. From this one can extract that, in the real network, either the activation of RD species is done earlier than Hth’s in the developmental program, or the subnetwork involving RD and Hth species is more complex than a simple mutual interaction. The full analysis of the oscillatory regime and the dynamical properties of the network species involved will be released by the authors in another work.
3 Materials and Methods
The total onedimension ocellar complex is built on a linear lattice of 117 nodes. It has been measured for this work that the IOC and each ocelli is approximate long. The ocellar complex is approximately 30 cells long and the cell diameter is . A good approximation in the model is to set each lattice node as a of tissue length. So, it can be assumed that each element can reach 33 nodes long, thus leaving 9 nodes free on each side for fitting purposes. This way, a linear space of 117 nodes is enough to consider the development of a satisfactory ocellar complex pattern. The spatial dimension on the graphical results are, then, shown in units.
A simplified onenode model (for the six equations (2732) case) was implemented using Vensim software, a visual tool for solving ODEs that allows parameter values modification in runtime (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 117nodes models were also implemented using MATLAB (MathWorks) and solved with the integrator ode45.
For all the numerical experiments, the differential equations system has been discretized using a finite differences method.
The values for all the parameters used in the simulations and the initial conditions for each variable can be found respectively in tables 1, 2 for the nonlinear models and 3 and 4 for the linear models in the Supplementary Material. Some parameter values are cited and others have been mathematically derived. The rest of the parametric set has been manually adjusted to obtain a similar ocellar pattern to the experimental one.
The Gaussian fitting for Hh steady state profile was performed by Microcal Origin software. The Boolean networks analysis was performed using application BooleSim Bock et al. (2014).
4 Discussion
How to establish a pattern in biological tissues is one of the most hot topics in recent developmental biology, though its interest has been shown from several decades from both experimental Wolpert (1994); Lum (2004); Dessaud et al. (2010); Kicheva et al. (2007); Nahmad and Stathopoulos (2009); Wartlick et al. (2011); Restrepo et al. (2014) and theoretical perspectives Collier et al. (1996); Meinhardt (1996); Barkai and Leibler (1997); BenZvi and Barkai (2010); Nahmad and Stathopoulos (2010); Burda et al. (2011); Nahmad (2011); Kicheva et al. (2012); Li et al. (2012). Several mechanisms for pattern formation have been proposed and this research is still ongoing with many unanswered questions. In the present work we are focusing on the robustness of the ocellar complex pattern in the Drosophila melanogaster fly, in terms of perturbations in the system on the description level of the GRN structure representation. We started from a high description level PDE model using a nonlinear formulation AguilarHidalgo et al. (2013). We have systematically modified this model by considering simpler assumptions, which keep the patterning, while maintaining the basic biology from a coarsegrain point of view. This study reveals that this pattern follows one of the most known patterning mechanisms, the French flag model Wolpert (1968). This mechanism is driven by a morphogen gradient (Hh) where the cellular fate is determined by the concentration of this diffusive molecule, triggering different parts of the GRN. The hhsignaling pathway is blocked by the expression of its inhibitor, En, which is also Hhdependent. This inhibitor becomes Hhindependent when its concentration exceeds a certain threshold, selfregulating and repressing the expression of the RD genes. This way, at high Hh concentration, the En domain is set and the IOC is formed. Otherwise, En is not able to inhibit the expression of the RD genes and the OC region is fulfilled. This basic mechanism seems to be robust in this GRN even after being substantially pruned. We have found that the patterning itself is independent, for example, on the receptor dynamics and on how the downstream signaling molecule CiA forms a graded profile. In this latter, we found that ci expression could be directly or indirectly dependent on the concentration of receptorbound morphogen.
In order to test if the patterning is significantly affected by the type of formulation used to simulate the model, we have implemented this model using both a nonlinear and a linear paradigm. This result showed that the patterning seems to be also robust against the choice of formulation paradigm. Comparing the results of both types of formulation one can see that the nonlinear model allows a precise control on the activation and repression processes by using Hill functions. This permits, for example, to completely block the Hhsignalling pathway by regulating the Hill function parameters on each regulatory term. This control is less clear in the linear formulation where the parameters contribute to the rate of change of concentration by adding/subtracting concentration dependant terms. On the contrary, the parameters in the linear model tell about effective rates that, in some cases, could be experimentally measured and give feedback to the model. One can note that the repression transition is smoother in the nonlinear than in the linear formulation due to the use of a step function in the selfregulation of En in the linear case, eq. (33). We also found that the correct patterning is fulfilled when considering both a dynamical and a steady morphogen profile. In the dynamical description, the system can positively describe different developmental stages during Hh profile transition to its steady state from a null concentration scenario. On the other hand, as an adiabatic approximation also retrieves the correct pattern, we can infer that the morphogen dynamics is indeed not necessary to fulfill the pattern. This implies that the morphogen kinetics could be faster than the intracellular response so we can consider the morphogen dynamics as a serial of steadystate profiles.
By analyzing the minimal expression of the GRN capable of maintaining this behavior, we have found that the patterning is dependent on the specific topology of a core regulatory network motif containing an activatorrepressor regulatory mechanism. This activatorrepressor paradigm is a well known mechanism in pattern formation since Turing (1952). In our case of study, the activatorrepressor mechanism is local for each cell, representing a triggering structure for a classical French flag model, where the activatorrepressor network motif behavior responds to the morphogen concentration input in each cell. Though mathematically speaking, this activatorrepressor mechanism is enough to fulfill the correct pattern by itself, biologically speaking, Hh cannot directly reach the activation of CiA but needs at least one intermediate to describe an extraintracellular communication. This way, we consider the minimal representation of the GRN, with biological significance, the (core) 3nodes subnetwork described in this work. Here PtcHh is defined as the intermediate between Hh and CiA. This is, then, a network motif with a double negative feedback loop. It is also known that this type of threecomponent systems are appropriate to account for a wide class of biological phenomena Meinhardt (2004). Some examples are found in the enrichment of microRNA targets Cui et al. (2006); Marr et al. (2010). A closer example to the system studied here is the developmental patterning of neural subtype specification in mouse Balaskas et al. (2012). There, the morphogen responsible for this patterning is the mouse Hh homologue Sonic Hedgehog (Shh). It is likely to think that this specific type of 3nodes activatorrepressor network motif can be responsible for the development of similar patterns in other developmental systems.
We have also confirmed the patterning independence on the dynamical model implementation complexity by analyzing the four studied network models from a Boolean perspective, and thus, the robustness of the GRN behavior and their patterning mechanism as long as the core network topology is maintained. We found that the IOC region is established in a steady state fashion. Interestingly, the OC region is defined as an oscillatory regime in the four models. Specifically, the activatorrepressor mechanism is responsible for the emergence of these oscillations. The inclusion of En as a repressor of the hh signalingpathway establishes a negative feedback within the network that might be capable of changing the system dynamics nonlinearly causing oscillations in gene expression. An experimental research for validating these oscillations during ocellar development has not been explored yet, though it is a tantalizing one. Oscillatory regimes have been proposed to increase the efficiency of reactions Richter and Ross (1980), and to be involved in cell patterning by a morphogen graded profile Koch and Meinhardt (1994). Even this oscillatory behavior could act as a switchlike mechanism between two transient cellstates building spatial borders among oscillating and nonoscillating cells.
Acknowledgements
The authors would like to thank Fernando Casares (CABD) for introducing them to this fascinating system and for many fruitful discussions. This work is partially financed by Junta de Andalucía (FQM122) to A. Córdoba and by the Spanish Ministry for Science and Innovation (MICINN/MINECO) and Feder Funds through grants BFU201234324 to F. Casares (CABD, Seville, Spain) and Consolider ÔFrom Genes to Shape,Õ of which F. Casares was a participant researcher.
Node  Boolean updating function 

Model  
for all  
AND NOT AND NOT  
AND  
NOT  
AND NOT OR  
AND (NOT AND )  
AND NOT ) OR AND  
OR ( AND NOT )  
AND NOT  
for all  
Model  
for all  
OR AND NOT  
AND NOT  
OR AND  
OR AND NOT  
NOT  
Model  
for all  
OR AND NOT  
AND NOT  
OR AND  
Model  
for all  
AND NOT  
OR AND  
Appendix A Analysis of the cooperation reactions in CiA and En
(41) 
(42) 
In both equations, the production rate term is written as a Hill function, typically represented as:
(43) 
where is the concentration producing half ; is the maximum value for and is the Hill coefficient which describes the cooperativity or affinity to continue the reaction defined with the Hill equation once started. It is known that means negative cooperation (the reaction affinity decreases), means noncooperative reaction and means positive cooperation (increasing the reaction affinity). The cooperation study was done by modifying the Hill coefficients value and observing the RD output profile, while maintaining unchanged the rest of the parameters. RD expression pattern is only able to correctly define the IOC and OC regions when both and (Fig. 6a). This is, just in the defined Hill coefficients intervals En is able to achieve its selfregulation state transiting from one to two stable states. When increasing over 4, the RD expression profile increases its intensity extending across the entire ocellar region (Fig. 6b). Giving values for , the RD ocellus profile sharpens acquiring a more similar to square shape (Fig. 6c). Hereinafter, the Hill coefficient parameters values will be and .
Appendix B Stability analysis for nonlinear contribution equations
The aim of this stability analysis is double. In the first place this study can give a valid range of values for parameters. And secondly, this analysis will show the system behavior before and after the beginning of the En selfregulation, from a stability point of view. There are two differential equations in which a formulation of the Hill’s equation form has been used: CiA’s (equation 41) and En’s (equation 42).
b.1 Stability analysis for CiA’s equation
In order to achieve a stability analysis on CiA behavior we performed a simplification on equation 41.
(44) 
This equation (44) represents the behavior of CiA before including En into the system. This way we can also analyze the behavior of the system when introducing this repressor. It is also possible to differentiate between a linear and a nonlinear repression addition. Due to the high nonlinearity of the introduced equations system, we are proposing to consider PtcHh as a direct contributor to the CiA production, including the effect of both variables into just one , that will represent the full output of the Hh signaling pathway. This way, equation (44) can be reduced to
(45) 
The value range was taken from the numerical simulations and it is built from the union of the PtcHh and CiA value intervals. This equation can also be expressed as:
(46) 
where is the production rate and is the degradation rate. Fig. 7a shows a bistable switch where the intersection of the two tendencies exhibit two steady states (production rate is equal to the degradation rate); these steady states are stable points (in blue). So, the Hhsignaling pathway can either be closed (low CiA concentration unable to regulate RD) or open (high CiA concentration that upregulates RD).
In this stability analysis, it was used the same equation but adding En in two possible ways: as a nonlinear repression (equation 41) and as a linear repression (equation 47):
(47) 
Repeating the same procedure to these two equations:
(48) 
With expression 48 for the nonlinear repression form, and
(49) 
being the expression 49 for the linear repression form.
Two different values were chosen for En concentration depending on its state: without selfregulation (1.51) or with selfregulation (3502.49). These values were taken from numerical simulations being closed to a steady state. In order to obtain a value for that ensures at least one stable steady state, the second derivative of must be negative in a valid value pair. This value pair is the tangent point between and . The first derivative from equation 48 is:
(50) 
From equations (48) and (50) is obtained, as the only real and positive solution, that the production and degradation rate are tangent at . For these values the system ensures a positive stable node and a valid CiA behavior. This is, CiA concentration is high when En is not selfregulated (black curve in Fig.7b), and CiA concentration is very low or null for high concentration of En (when En is selfregulated; blue curve in Fig. 7b). So the closure of the Hhsignaling pathway can be achieved when En represses CiA in a nonlinear way. The second derivative is negative at the calculated and value, ensuring its convexity.
(51) 
For , .
When evaluating the En repression action to CiA in a linear way (eq. 49), we find the same behavior independently of En activation. Either En being at low or high concentration it is shown one unstable steady state at low values and a stable steady state at high levels (Fig. 7c). With this linear repression the Hhsignaling pathway can remain open with a high concentration of CiA in the cells where En is selfregulated. Moreover, as the steady state found at low values is unstable, perturbations applied to the system can lead its behavior to the stable, and high CiA concentration node. In other words, if the system is able to close the Hhsignaling pathway due to a linear En repression, it can be spontaneously open again if the system is perturbed. This behavior does not correspond to the CiA experimentally observed in AguilarHidalgo et al. (2013). It is known that CiA reaches a high concentration positive stable node when En is not selfmaintained and that it is overridden when En selfregulation is activated. In the case of nonlinear repression (eq. 48) shows a bistable state when En is not selfregulated. One stable node at and a second one at . When En selfregulation starts loses its bistability maintaining just the stable steady state at null levels (Fig. 7b). This behavior fits with the experimental one. The only possible stable node of CiA when En selfmaintains is the closure of Hhsignaling pathway, this is, CiA concentration falls down.
b.2 Stability analysis for En’s equation
A similar analysis can be done to the differential equation that drives the behavior of En (equation 42). If En is not selfregulated (), this equation can be expressed as:
(52) 
The z(x) first derivative of equation 52 is:
(53) 
The tangent point in this case is and .
The second derivative (eq. 54) shows that the tangent point is a stable one as as corresponds to a maximum.