Core regulatory network motif underlies the ocellar complex patterning in Drosophila melanogaster

Core regulatory network motif underlies the ocellar complex patterning in Drosophila melanogaster

D. Aguilar-Hidalgo Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Strasse 38, 01187 Dresden, Germany.
Departamento de Física de la Materia Condensada, Universidad de Sevilla.
41012 Sevilla, Spain.
CABD (CSIC-UPO-Junta de Andalucía), 41013 Sevilla, Spain.
M.C. Lemos A. Córdoba

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 3-nodes 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.

Pattern formation, ocellar development, Reaction-diffusion model, Boolean networks

1 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 dorsal-anterior 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 one-dimensional description is then proposed where two ocelli regions are separated by an interocellar cuticle. A physical model based on reaction-diffusion equations has already been proposed in Aguilar-Hidalgo 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:

  1. Robustness with respect to changes in the model parameters.

  2. Robustness with respect to transient changes.

  3. Robustness with respect to changes in the structure of the equation itself.

Types 1 and 2 have been tested in Aguilar-Hidalgo 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 Aguilar-Hidalgo 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 Aguilar-Hidalgo 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 ligand-receptor 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 Aguilar-Hidalgo 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 non-linear 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 three-nodes network motif that satisfies the correct patterning by itself. This core network is biologically very relevant as it defines the main morphogen-signaling 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 activator-repressor 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 activator-repressor 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 one-dimensional 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 POC-OC-IOC-OC-POC (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 anterior-posterior 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 self-regulates 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 Aguilar-Hidalgo 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 Aguilar-Hidalgo 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 Aguilar-Hidalgo et al. (2013) for the GRN in (Fig. 1c) is implemented using a non-linear formulation and is given by the following equations:


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 non-dimensional parameters , , and are used for changing the scale of different terms and is the diffusion coefficient. Subscript X-Y, 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 Aguilar-Hidalgo 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 Aguilar-Hidalgo 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 non-linearly represented. These two terms were repressed by a non-linear 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:


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 re-inserting 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 non-linear 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. (6-7). In a coarse-grain 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 Hh-like profile repressed by En when inserting a contribution:


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 Hh-like profile in eq. (18), one can assume that CiA is self-stable and thus neglect the transformation from CiA to CiR in eqs. (6-7), and so, eq. (7) for CiR can be eliminated. CiA dynamics is shown in eq. (19).


At this point, the system stands for another change in a coarse-grain view. One can avoid the description of the receptor-ligand 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 Hh-signalling 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:


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. (20-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 switch-like 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 coarse-grain 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).

Figure 1: (a) Ocellar complex scheme showing the two posterior ocelli (PO) and the anterior ocellus (AO). The region among the ocelli is called the interocellar region. The black triangle represents the triangular cuticle that embraces the ocellar complex. (b) Scheme of the one-dimensional simplification of the ocellar complex. The two ocelli (AO and PO) are considered identical (OC) so the model is symmetric. The one dimension ocellar complex is formed by three different types of tissues represented by POC, or periocellar region (in blue), OC, or ocellar region (in green) and IOC, or interocellar region (in purple). (c-d) Network diagrams of the ocellar model representing the system’s variables (circle nodes) and constant inputs (square nodes) linked by positive (green arrows) and negative (red arrows) interactions. (c) GRN in Aguilar-Hidalgo et al. (2013). Dashed lines from Hh nodes implies that Hh contributes to the formation of the PtcHh complex by repressing the number of free receptor Ptc. (d) Simplified GRN.

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 Aguilar-Hidalgo 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 re-written 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).


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 (27-32) are the reaction-diffusion type differential equations of the model. The description of the system is the same as in the last section. Some non-linearities 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 non-linear 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 signaling-independent manner. As CiA is the only En activator, En should self-maintain. 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 self-maintenance should work only under certain circumstances, and these are Hh signaling-dependent. In order to treat this feature in a simplistic way En is self-maintained 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 self-regulate, 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); Aguilar-Hidalgo 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 self-regulation term will not cause any effect and RD will appear in only one domain over the whole tissue.

Figure 2: (a) Representation of RD and En normalized profiles in the case that En is not self-regulating. In this case, the RD is present in the whole ocellar complex appearing as one unsplit domain. The parameter values used to obtain this result are shown in Table 2 of Supplementary Material, modifying parameter value to . (b) Representation of RD and En normalized profiles in case En is self-regulated just in the cells where its concentration is over a determined threshold (). Now RD defines the ocellar domains while En characterizes the IOC.

The pruned GRNs both linearly and non-linearly 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 self-regulating 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 Hh-signaling 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 non-linearly done one. This refers to the use of a step function for the En self-regulation (linear model) instead of a Hill’s equation (non-linear model). As in Aguilar-Hidalgo 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).

Figure 3: Spatio-temporal representation of the signal intensity of the non-linear (left - eqs. (20-25)) and linear (right - eqs. (27) to (32))) simplified models. From up to down: Hh, PtcHh, CiA, En, Eya/RD, Hth. The brighter the color, the higher the concentration value.

2.3 A 3-nodes 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 3-node 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 gradient-like 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.


These equations (34-36) 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 3-nodes 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.

Figure 4: (a) Network diagram of the 3-nodes model. The network dynamics is the following: Hh is the input to the network in the form of a static Gaussian profile that activates PtcHh. PtcHh activates CiA while CiA activates both PtcHh and En. En represses both PtcHh and CiA. En is self-maintained in the case Dl function is activated. Dl activates En self-regulation if En surpasses a certain concentration. (b) Hh gradient steady state in the model from the pruned GRN (in black line) and its Gaussian fit (in dashed red line). (c-e) Solutions of the 3-nodes model of 3 equations for (c) PtcHh, (d) En and (e) CiA.

2.4 Reduction of the system to an activator-inhibitor dynamics

The 3-nodes 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 activator-inhibitor interaction. The patterning is correctly achieved if this specific set is maintained in the core of the network.

Following the 3-nodes model formulation it is possible to write this system in terms of ODEs as follows:


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:


A simple observation of eq. (40) shows that in the case of En self-regulation, 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 activator-repressor 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 mean-field 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 Aguilar-Hidalgo et al. (2013) (model ), the one in Fig. 1d (model ), the described in Fig. 4a (model ) and the activator-repressor 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 self-regulates, the Hh-signaling 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 self-regulate. 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 non-steady 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 switch-like 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 RD-Hth 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 sub-network 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.

Figure 5: Time series showing the behavior of the different network models from their Boolean perspective (black = 1, grey = 0). (a) Situation in which the Dl/Notch signaling pathway trigger the En self-regulation. The four models show a steady-state solution where the Hh-signaling pathway is blocked by En defining the IOC. (b) Scenario where En is not self-regulated and the OC region is defined by CiA and Eya in an oscillatory regime. Only model #2 shows a steady situation for RD as Hth persistently represses it. Red lines show the oscillation period.

3 Materials and Methods

The total one-dimension 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 one-node model (for the six equations (27-32) case) 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, This program contains a fourth order Runge Kutta method (RK4) to solve ODE systems.

The 117-nodes 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 non-linear 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); Ben-Zvi 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 non-linear formulation Aguilar-Hidalgo et al. (2013). We have systematically modified this model by considering simpler assumptions, which keep the patterning, while maintaining the basic biology from a coarse-grain 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 hh-signaling pathway is blocked by the expression of its inhibitor, En, which is also Hh-dependent. This inhibitor becomes Hh-independent when its concentration exceeds a certain threshold, self-regulating 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 receptor-bound 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 non-linear 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 non-linear model allows a precise control on the activation and repression processes by using Hill functions. This permits, for example, to completely block the Hh-signalling 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 non-linear than in the linear formulation due to the use of a step function in the self-regulation 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 steady-state 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 activator-repressor regulatory mechanism. This activator-repressor paradigm is a well known mechanism in pattern formation since Turing (1952). In our case of study, the activator-repressor mechanism is local for each cell, representing a triggering structure for a classical French flag model, where the activator-repressor network motif behavior responds to the morphogen concentration input in each cell. Though mathematically speaking, this activator-repressor 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 extra-intracellular communication. This way, we consider the minimal representation of the GRN, with biological significance, the (core) 3-nodes sub-network 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 three-component 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 3-nodes activator-repressor 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 activator-repressor mechanism is responsible for the emergence of these oscillations. The inclusion of En as a repressor of the hh signaling-pathway establishes a negative feedback within the network that might be capable of changing the system dynamics non-linearly 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 switch-like mechanism between two transient cell-states building spatial borders among oscillating and non-oscillating cells.


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 (FQM-122) to A. Córdoba and by the Spanish Ministry for Science and Innovation (MICINN/MINECO) and Feder Funds through grants BFU2012-34324 to F. Casares (CABD, Seville, Spain) and Consolider ÔFrom Genes to Shape,Õ of which F. Casares was a participant researcher.

Node Boolean updating function
for all
for all
for all
for all
for all

Table 1: Boolean functions

Appendix A Analysis of the cooperation reactions in CiA and En

Equations (41) and (42) show the temporal evolution of CiA and En concentrations respectively.


In both equations, the production rate term is written as a Hill function, typically represented as:


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 non-cooperative 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 self-regulation 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 .

Figure 6: (a) States plot showing the changes in En behavior for different values of the Hill coefficients and in equations (41) and (42). The vertical axis is the number of stable states that En can reach as a function of the Hill coefficient ( is either as ). Black line represents the number of stable states that En reaches when modifying . When , En shows just two states while for higher values, En just shows one stable state without transition to its self-regulation. Red line represents the number of stable nodes that En reaches when modifying . In this case, En shows the transition to a second state from on. The combination of the two Hill coefficients valid range gives a system correct behavior for and (b) RD profiles for different values of . For values higher than 4 () En is not able to self-regulate and the RD expression patterns covers the whole ocellar complex at high level. Lower values show a correct RD profile lowering its intensity while decreasing . (c) RD profiles for different values of . Only for the RD profile fulfill the correct pattern. As the value increases the RD pattern sharpens

Appendix B Stability analysis for non-linear 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 self-regulation, 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.


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 non-linear repression addition. Due to the high non-linearity 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


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:


where is the production rate and is the degradation rate. Fig. 7a shows a bi-stable 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 Hh-signaling pathway can either be closed (low CiA concentration unable to regulate RD) or open (high CiA concentration that up-regulates RD).

In this stability analysis, it was used the same equation but adding En in two possible ways: as a non-linear repression (equation 41) and as a linear repression (equation 47):


Repeating the same procedure to these two equations:


With expression 48 for the non-linear repression form, and


being the expression 49 for the linear repression form.

Two different values were chosen for En concentration depending on its state: without self-regulation (1.51) or with self-regulation (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:


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 self-regulated (black curve in Fig.7b), and CiA concentration is very low or null for high concentration of En (when En is self-regulated; blue curve in Fig. 7b). So the closure of the Hh-signaling pathway can be achieved when En represses CiA in a non-linear way. The second derivative is negative at the calculated and value, ensuring its convexity.


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 Hh-signaling pathway can remain open with a high concentration of CiA in the cells where En is self-regulated. 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 Hh-signaling 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 Aguilar-Hidalgo et al. (2013). It is known that CiA reaches a high concentration positive stable node when En is not self-maintained and that it is overridden when En self-regulation is activated. In the case of non-linear repression (eq. 48) shows a bi-stable state when En is not self-regulated. One stable node at and a second one at . When En self-regulation starts loses its bi-stability 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 self-maintains is the closure of Hh-signaling pathway, this is, CiA concentration falls down.

Figure 7: Stability analysis for non-linear contribution equations. (a) This graph represents the production rate (black line) and the degradation rate (red line) of (vertical axis: , in eq. 46) as a function of the concentration of (horizontal axis), when En is not considered as a CiA repressor. The points at which the production rate is equal to the degradation rate are considered steady states. This configuration shows two different steady states (two stable nodes as blue dots). (b) Representation of the system steady states when En is repressing CiA in a non-linear way. Black line shows the production rate when En is in its lower chosen value (En = 1.51, before beginning its positive self-regulation). Blue line (very close to horizontal axis) represents the production rate when En has reached its higher chosen value (En = 3502.49, after En self-regulation response). Red line represents the degradation rate. In the first case the system maintains its bi-stability and in the second case the system loses the bi-stability maintaining just one stable node at . (c) Representation of the system steady states when En is repressing CiA in a linear way. Black line represents the production rate when En remains in its lower chosen value. Orange line represents the production rate when considering En on its higher chosen value. Red line represents the degradation. The system shows two steady states on each configuration, one stable at high value and one unstable at low value (red dot). If the system stays at the unstable steady state, small perturbations will be amplified letting the system reach the stable node. The linear repression makes the system lose its bi-stability and does not allow any stable state at low value to consider closing the Hh signaling pathway. (d) Graphical determination of the first stable steady state for En evolution (equation 52) showing the intersection between the production rate (black line) and the degradation rate (red line). The blue line represents the production rate when En self-regulates (equation 55). In this case, the system loses its bi-stability remaining just one stable node at high En concentration. (e) Representation of the En production rates (black line) considering the chosen value for , and the degradation rate (red line). The system shows a bi-stable state with null and high levels, also considering one unstable state at low level

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 self-regulated (), this equation can be expressed as:


The z(x) first derivative of equation 52 is:


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.