Core regulatory network motif underlies the ocellar complex patterning in Drosophila melanogaster
CCore regulatory network motif underlies the ocellar complexpatterning in
Drosophila melanogaster
D. Aguilar-Hidalgo a,b,c, ∗ , M.C. Lemos b , A. C´ordoba b a Max Planck Institute for the Physics of Complex Systems,N¨othnitzer Strasse 38, 01187 Dresden, Germany. b Departamento de F´ısica de la Materia Condensada, Universidad de Sevilla.41012 Sevilla, Spain. c CABD (CSIC-UPO-Junta de Andaluc´ıa), 41013 Sevilla, Spain.
Abstract
During organogenesis, developmental programs governed by Gene Regulatory Networks(GRN) define the functionality, size and shape of the different constituents of living or-ganisms. Robustness, thus, is an essential characteristic that GRNs need to fulfill in orderto maintain viability and reproducibility in a species. In the present work we analyze the ro-bustness 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 systemto obtain the minimum pathway able to satisfy this pattern. We found that the mecha-nism underlying the patterning obeys to the dynamics of a 3-nodes network motif with adouble negative feedback loop fed by a morphogenetic gradient that triggers the inhibitionin a
French flag problem fashion. A Boolean modeling of the GRN confirms robustness inthe patterning mechanism showing the same result for different network complexity levels.Interestingly, the network provides a steady state solution in the interocellar part of thepatterning and an oscillatory regime in the ocelli. This theoretical result predicts that theocellar pattern may underlie oscillatory dynamics in its genetic regulation.
Keywords:
Pattern formation, ocellar development, Reaction-diffusion model, Booleannetworks
1. Introduction
The study of biological systems from a physical perspective has substantially increasedin the last decade as a result of the same increment in the synergy between theoretical andexperimental researchers in this field. Physical models help on a better understanding of theprinciples of biological systems. One of the biological areas that counts with a higher number ∗ Corresponding author
Email addresses: [email protected] (D. Aguilar-Hidalgo), [email protected] (M.C. Lemos), [email protected] (A. C´ordoba)
Preprint submitted to Elsevier October 4, 2018 a r X i v : . [ q - b i o . T O ] D ec f theoretical collaborations is developmental biology. This biology field, in the case ofmulticellular organisms, studies how these organisms develop from a single cell to a complexstructure with many differentiated cell types capable of performing distinct functions. Here,the time parameter gives an extra dimension able to show complex emergent properties, andthis is indeed very attractive for physicists interested in the study of complex systems.Evolutionarily speaking, every developmental process is subjected to mutations that, insome 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 toperturbations in its developmental program. Organisms containing a high level of complexityare able to maintain certain modifications between individuals while maintaining robustnessin the essential processes. Low complex organisms usually present a clear robustness in allthe developmental processes. A noteworthy example is that of the Caenorhabditis elegans worm. This worm has exactly 959 somatic cells in adult hermaphrodites and 1031 somaticcells in adult males [1], 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 specifictissue patterning, the ocellar complex. The ocellar complex is part of the visual system ofmany insects. In this fly, it is formed by three ocelli (or simple eyes), situated at the verticesof a triangular patch of cuticle, and their interocellar region, placed in the dorsal head. Thissingular patterning is actually formed by the fusion of the dorsal-anterior domains of thetwo eye discs [2], so it can indeed be studied as the two ocelli pattern of one of the two eyeimaginal discs. A one-dimensional description is then proposed where two ocelli regions areseparated by an interocellar cuticle. A physical model based on reaction-diffusion equationshas already been proposed in [3] to fulfill, not only the tissue patterning for the Drosophilaocellar complex, but also the correct experimentally tested behavior of the gene regulatorynetwork (GRN) that drives its formation.This specific pattern follows the
French flag problem paradigm, where different devel-opmental programs are run depending on the concentration of an extracellular diffusivemolecule (morphogen) [4].Following Kang et al. [5], in order to test robustness in a physical model we can identifythree 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 [3]. The equations there described provide the samepattern when the system is subject to noise perturbation and changes in the initial condi-tions. Moreover, an extensive parameter sensitivity analysis was performed, where pertur-bations in the parameter values retrieved the same patterning with changes in size of itscomponents [5], what suggests that the same GRN determines a certain phenotypic varia-tion. 2he third type of perturbations considers structural changes in the equations describinga physical model. In a networked system, an exercise to test this type of robustness isto modify the wiring of the network by adding/removing links like in [6]. This, in a GRNmodifies/removes regulatory interactions. As mentioned above, many regulatory interactionsrelated to the ocellar complex developmental GRN were tested in [3], and agreed with themodel. However, there is an important issue related to structural robustness of the systemremaining untested. This is the model description level [7]. As the patterning retrievedby this model is robust facing parameter values modification, noise, changes in the initialconditions and mutations, we are interested to know the reason for such robustness in thepatterning retrieval.In the present work we focus in perturbations on the structure of the equations withrespect to the description level of the physical model. The model presented in [3] describesthe 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 anda nonlinear formulation. Here, we have reduced the description level of the model whilemaintaining the biological sense according to the results in [3], and thus, we have tried tofind out the minimum network expression capable of showing the correct pattern. To dothis, we have systematically removed nodes and links and, then, changed the formulationfrom a non-linear to a linear paradigm.The patterning seems at this point to be independent on the description level and formu-lation complexity. Both the patterning and the individual behaviors of each GRN componentis correctly fulfilled. An observation on the resulting GRN leads to the extraction from thisGRN of a core formed by a three-nodes network motif that satisfies the correct patterningby 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 repressorin a double negative feedback loop. A further analysis shows that even the reduction of thiscore network to a simple activator-repressor feedback loop retrieves the desired patterning.All the different models described are then simulated using Boolean networks. This changesthe description level of the dynamical system to its simplest representation [7]. This exer-cise confirms that this patterning is independent on the description level but dependent onspecific regulatory topology of the activator-repressor motif.
2. Results
During fly development, the ocellar complex is formed describing a triangular patch witha 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 twoocelli, one PO and one AO. In a late larval stage the two anterior ocelli are fused intoone and only anterior ocellus. In order to describe the ocellar complex conformation in asimplified way, it can be reduced to a one-dimensional description by analyzing one row ofcells situated in a cross line over the two ocelli of one of the two eye discs. This conformationshows a 1D field compound of five blocks of three different types of tissues. The external3issue will be referred hereinafter as periocellar region (POC). The following inner type oftissue corresponds to the region where the ocelli will be placed; this tissue will be referredas 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 mor-phogen profile. One of the most studied morphogens is the molecule Hedgehog (Hh), whichhas already been considered to successfully describe intercellular communication in devel-opmental processes [8]. This morphogen is known to be essential for the formation of theocellar complex in
D. mel. [9]. Here, Hh is expressed in a central stripe of cells in the tar-get tissue diffusing in the anterior-posterior axis. The readout of this system is the retinaldetermination (RD) genes group. These genes, eyes absent ( eya ) and sine oculis ( so ), arerequired for the formation of the ocelli [10–16]. Experimental work in Drosophila and invertebrates indicates that the Hh signaling pathway is subject to extensive feedback amongelements of the pathway, most notably that of the Hh receptor patched ( ptc ) [17, 18]. Thegene engrailed ( en ) is activated by the Hh signaling pathway and it self-regulates only inthe cells where its protein (En) concentration is high. This En high concentration regionoverlaps with the Hh expression cells though the En domain can vary depending on externalconditions and variations in the GRN [3]. En, then, attenuates the Hh signaling pathway toestablish a region with no RD expression (IOC) and thus, the ocellar pattern.An extensive experimental and theoretical work in the ocellar complex shows robustnessin the patterning against noise, changes in the initial conditions and parameter values [3].There, some structural modifications were also tested by removing links in the regulatorynetwork. As example, the ocellar pattern is maintained with enlarged ocelli if the repressionbetween Homothorax (Hth) and eya is removed; on the contrary, it is noteworthy to mentionan example by which the ocellar pattern is not fulfilled, this takes place when removing thelink between Delta (Dl) and En; here the anterior and posterior ocelli fuse and just one bigocellus is developed. These structural modifications were experimentally validated and arealso considered in the present work. The proposed model in [3] for the GRN in (Fig. 1c) isimplemented using a non-linear formulation and is given by the following equations: ∂Hh∂τ = D ∂ Hh∂x + δ ( x ) α hh − γ P tc − Hh P tc · Hh − β Hh Hh (1) ∂ptc∂τ = κ β ptc (( α ptc + κ Ciptc φ ( CiA ψ ( CiR, k
CiR , n
CiR ) , k CiA , n
CiA )) ψ ( En, k En − ptc , n En − ptc ) − ptc )(2) ∂P tc∂τ = θ ptc ptc − γ P tc − Hh P tc · Hh − β P tc
P tc (3) ∂P tcHh∂τ = γ P tc − Hh P tc · Hh − β P tcHh
P tcHh, (4)4 ci∂τ = κ β ci ( α ci φ ( ψ ( En, k En − ci , n En − ci ) , k ci , n ci ) − ci ) (5) ∂CiA∂τ = κ β CiA ( θ ci ci − CiA ) − κ Ci φ ( CiA ψ ( P tcHhP tc , k
P H , n
P H ) , k CiA , n
CiA ) (6) ∂CiR∂τ = κ Ci φ ( CiA ψ ( P tcHhP tc , k
P H , n
P H ) , k CiA , n
CiA ) − κ β CiR
CiR (7) ∂en∂τ = κ β en ( α en + φ ( CiAψ ( CiR, k
CiR , n
CiR ) , k CiA , n
CiA ) + κ En φ ( En, k
DlEn , n En ) − en )(8) ∂En∂τ = θ en en − β En En (9) ∂eya∂τ = κ β eya ( α eya + α T oy φ ( CiAψ ( Hth, k
Hth − eya , n Hth ) , k CiA − eya , n CiA ) + φ ( Eya, k
Eya , n
Eya ) − eya )(10) ∂Eya∂τ = θ eya eya − β Eya
Eya (11) ∂hth∂τ = κ β hth ( α hth + α W g φ ( ψ ( Eya, k
Eya − hth , n Eya ) , k W g , n
W g ) − hth ) (12) ∂Hth∂τ = θ hth hth − β Hth
Hth (13) δ ( x ) = (cid:26) x ∈ hh -expressing cells0 if x / ∈ hh -expressing cells (14)where φ ( X, k, n ) = X n k n + X n and ψ ( X, k, n ) = 1 − φ ( X, k, n ), with X any system variable.The model contains different parameter types: α X for the basal transcription rates, β X for the degradation rates, k X for the Hill equation transcriptional regulators, n X for theHill coefficients, θ X for the translation rates; γ X for protein complex formation. The non-dimensional parameters κ , κ Ci , κ En and α ci are used for changing the scale of differentterms and D is the diffusion coefficient. Subscript X-Y, with X and Y system variables,indicates a regulation from X to Y. 5 .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[3] in order to test robustness in a different way as done there and to try to infer the basicmechanism responsible of the patterning. These changes include pruning of the GRN byeliminating certain nodes and links. These modifications assume changes in the behavioralparadigm of some regulatory elements that, while maintaining the overall concept, theysimplify the corresponding mechanisms. Additionally, changes in the formulation paradigmsare also being considered.
In order to find the basic mechanism to satisfy the ocellar pattern it is possible tosystematically reduce the GRN in [3] (Fig. 1c) by removing nodes and links as the model issimplified and the system of equations is reduced.To simplify this system we firstly focus on the receptor dynamics. This previous modelimplemented 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 CiArestricted 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 ofreceptors [19] in each cell with a certain binding ( γ P tc − Hh ) and unbinding ( γ − P tc − Hh ) rates.This way eq. (2) for ptc can be eliminated, the translational rate θ ptc = 0 in eq. (3) andeqs. (1, 3 and 4) can be written as follows: ∂Hh∂τ = D ∂ Hh∂x + δ ( x ) α hh − γ P tc − Hh P tc · Hh − β Hh Hh + γ − P tc − Hh P tcHh (15) ∂P tc∂τ = − γ P tc − Hh P tc · Hh − β P tc
P tc + γ − P tc − Hh P tcHh (16) ∂P tcHh∂τ = γ P tc − Hh P tc · Hh − β P tcHh
P tcHh − γ − P tc − Hh P tcHh (17)The result of this simplification is a replication of the ocellar pattern, as given in theprevious model, for every species involved except for Ptc and PtcHh due to the absence ofEn interaction. This modification reveals the need to preserve the interaction between Enand the receptor dynamics as, in order to maintain the full biological sense, PtcHh needs toshow the correct patterning. Note that this will be regulated in a subsequent modificationto the model by re-inserting En repression and CiA feedback.The next downstream step in the network is the activation of ci dependent proteins. Inthe former model, the gene ci is modeled to be expressed dependent on a non-linear Enrepression, eq. (5). Then, CiA is produced with a linear dependence on ci expression, eq.(6), which is turned into CiR form inversely dependent on the ratio of bound and unboundPtc to Hh, eqs. (6-7). In a coarse-grain view of this process, the activation of CiA correlateswith 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 all6aintain the desired pattern. Firstly, eq. (5) can turn from a constant production repressedby En to a Hh-like profile repressed by En when inserting a P tcHh/P tc contribution: ∂ci∂τ = κ β ci (cid:18) α ci φ ( P tcHhP tc ψ ( En, k En − ci , n En − ci ) , k ci , n ci ) − ci (cid:19) (18)This means that the ci expression is dependent on both the concentration of receptorsbound 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 thelatter case, an unknown species downstream the receptor dynamics could sense both, theconcentration of free and bound receptors and, accordingly use this information to activate ci expression.Now, as ci ’s output is already in a Hh-like profile in eq. (18), one can assume that CiAis 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). ∂CiA∂τ = κ β CiA ( θ ci ci − CiA ) (19)At this point, the system stands for another change in a coarse-grain view. One can avoidthe description of the receptor-ligand binding process if the downstream targets maintaintheir profile, as done in [20]. As it has already been commented, ci spatial form correlateswith PtcHh and Hh profiles. Therefore, the next proposed simplification assumes that thebound Hh can behave as an internalized molecule dependent on Hh concentration insteadof following the binding process. This, in the real system, can stand for the description ofa molecule situated downstream of PtcHh and upstream of ci in the Hh-signalling pathway.This other molecule senses the effect of Hh and can be altered by other elements as Enand CiA. To simplify notation, PtcHh is used as this internalized molecule, now sensitiveto En repression. This analysis resumes the need of including En interaction to the PtcHhspecies commented above. With this assumption, the first part of the pathway is modeledas follows: ∂Hh∂τ = D ∂ Hh∂x + δ ( x ) α hh − α P tcHh
P tcHh − β Hh Hh (20) ∂P tcHh∂τ = α Hh − P H φ ( Hh, k Hh , n Hh ) φ ( CiA ψ ( En, k En − P H , n En − P H ) , k CiA − P H , n
CiA − P H ) − β P tcHh
P tcHh (21)Note that an equation for unbound Ptc can be eliminated from the system as the bindingprocess is no longer described. Most of the constant basal productions were set to zero inthe original model. It has also been tested that the system keeps the same behavior whenneglecting them so they are removed from the equations to simplify the equations structureat this point of the network pruning.Now a straight forward simplification is performed to the system. The original modelshowed a detailed dynamics involving both genes expression and protein dynamics. One can7otice that in the steady state, the genes expression and the protein concentration differ in aproportionality constant. According to this, one can represent the system using just one ofthe element types adjusting this proportionality constant in each case. The system is thenshown in eqs. (20-25). ∂CiA∂τ = κ β CiA ( θ ci α ci φ ( P tcHh ψ ( En, k
Enci , n
Enci ) , k ci , n ci ) − CiA ) (22) ∂En∂τ = θ en ( φ ( CiAψ ( CiR, k
CiR , n
CiR ) , k CiA , n
CiA ) + κ En φ ( En, k
DlEn , n En )) − β En En (23) ∂Eya∂τ = θ eya ( α T oy φ ( CiAψ ( Hth, k
Hth − eya , n Hth ) , k CiA − eya , n CiA ) + φ ( Eya, k
Eya , n
Eya )) − β Eya
Eya (24) ∂Hth∂τ = θ hth ( α W g φ ( ψ ( Eya, k
Eya − hth , n Eya ) , k W g , n
W g )) − β Hth
Hth (25)The proposed GRN to fulfill the ocellar complex pattern is shown in (Fig. 1d). A sumup description of this network is the following. The morphogen Hh diffuses extracellularlyand forms the input profile of the system. Hh enhances an intracellular PtcHh concentrationwhich causes a decrease in Hh concentration. Then, PtcHh results in signaling for CubitusInterruptus (CiA) in its active form. CiA activates the inhibitor Engrailed (En) that appliesnegative feedback loops at two different points on the Hh signaling pathway, PtcHh andCiA. Finally CiA activates the formation of Eya which interacts with Homothorax (Hth) ina switch-like form, so Hth is present where Eya is not, and vice versa. Hth is activated byWingless (Wg). Though Wg is a well known morphogen, in these models it is considered as aconstant input. Wg is, then, not behaving as a graded morphogen but a constant input. Thistakes into account that its effect in the system would be qualitatively similar if its expressionis flat in the whole imaginal disc [21]. After taking the simplifications where many of theinteraction mechanisms are referred in a coarse-grain description, we will refer hereinafterto Eya as RD proteins, without particularly considering any of them but showing a generaloverview of the RD complex dynamics. Then, the spatial domains in the model where RD ispresent correspond to the ocellar cuticle. (See complete results in Supplementary Material).
In order to perturb the structure of the model’s equations, the system has been substan-tially pruned from that in [3]. Now one can question whether the patterning can be affectedby the formulation used. In order to overcome this issue, the model has been re-writtenusing a mostly linear formulation (eqs. (26) to (33)). This way, Hill’s function formulationor any other kind of relaxation kinetics used to simulate more realistic biological behaviorsare avoided except in necessary cases (see Supplementary Material). δ ( x ) = (cid:26) x ∈ hh -expressing cells0 if x / ∈ hh -expressing cells (26)8 igure 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 triangularcuticle that embraces the ocellar complex. (b) Scheme of the one-dimensional simplification of the ocellarcomplex. The two ocelli (AO and PO) are considered identical (OC) so the model is symmetric. The onedimension ocellar complex is formed by three different types of tissues represented by POC, or periocellarregion (in blue), OC, or ocellar region (in green) and IOC, or interocellar region (in purple). (c-d) Networkdiagrams of the ocellar model representing the system’s variables (circle nodes) and constant inputs (squarenodes) linked by positive (green arrows) and negative (red arrows) interactions. (c) GRN in [3]. Dashedlines from Hh nodes implies that Hh contributes to the formation of the PtcHh complex by repressing thenumber of free receptor Ptc. (d) Simplified GRN. Hh∂τ = D ∂ Hh∂x + δ ( x ) α hh − α P tcHh − Hh P tcHh − β Hh Hh (27) ∂P tcHh∂τ = Θ ( α Hh − P tcHh Hh + α CiA − P tcHh
CiA − α En − P tcHh En ) − β P tcHh
P tcHh (28) ∂CiA∂τ = γ PtcHh − CiA α En − CiA En α P tcHh − CiA
P tcHh n CiA k n CiA m CiA + α P tcHh − CiA
P tcHh n CiA − β CiA
CiA (29) ∂En∂τ = γ CiA − En CiA n En k n En m En + CiA n En + α En − En EnDl ( ζ ) − β En En (30) ∂RD∂τ = Θ ( α CiA − RD CiA − α Hth − RD Hth + α RD − RD RD ) − β RD RD (31) ∂Hth∂τ = Θ ( α W g − Hth − α RD − Hth RD ) − β Hth
Hth (32) Dl ( ζ En , ζ τ , τ ζ En ) = (cid:26) Dl if En ≥ ζ En , ∀ τ ∈ [ τ ζ En , τ ζ En + ζ τ ]0 otherwise (33)The coefficients α X − Y drive an action (activation or repression, depending on sign, + or − respectively) from source element X to target element Y; β X is the degradation coefficientof the element X; γ X − Y is the maximum Hill term value when X regulates Y.The proposed linear formulation is not conservative so some parameter sets could leadto negative concentrations. To avoid this issue Heaviside functions Θ are introduced. Inthis new system, equations (27-32) are the reaction-diffusion type differential equations ofthe 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 CiAis shown as a reduction of the maximum level that CiA can reach ( γ P tcHh − CiA ). A stabilitystudy 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 SupplementaryMaterial). Secondly, in order to achieve the desired pattern, En needs to be maintainedin 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 completeEn domain with no response to Hh influence (Fig. 2a); so the En self-maintenance shouldwork only under certain circumstances, and these are Hh signaling-dependent. In order totreat this feature in a simplistic way En is self-maintained if its concentration exceeds acertain threshold, ζ En , for a time interval, ζ τ , starting from the time point when En reachesthe concentration threshold, τ ζ En , (eqs. A.2 and 33). This requirement refers to the need ofintegrating a protein quantity over time to trigger some targets [2].These premises form an En domain in the IOC allowing the production of the RD proteinsin the OC region (Fig. 2b). If En transcription is eliminated ( γ CiA − En = 0) or substantially10 igure 2: (a) Representation of RD and En normalized profiles in the case that En is not self-regulating. Inthis case, the RD is present in the whole ocellar complex appearing as one unsplit domain. The parametervalues used to obtain this result are shown in Table 2 of Supplementary Material, modifying Dl parametervalue to Dl = 0. (b) Representation of RD and En normalized profiles in case En is self-regulated just in thecells where its concentration is over a determined threshold ( Dl = 1). Now RD defines the ocellar domainswhile En characterizes the IOC. reduced so that it cannot reach enough concentration to self-regulate, then RD appears ina single unsplit domain that follows the Hh profile. There is experimental evidence showingthat 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, RDexpression and Hh signaling are detected in two domains. These domains are separated by aregion in which RD and Ptc expressions are very low or absent [2, 3]. This pattern suggestedthe existence of a repressor capable of shutting down the hh pathway. The same pattern isobserved if Dl signaling were shut off ( Dl = 0), the mathematical effect would be identical,En self-regulation term will not cause any effect and RD will appear in only one domainover the whole tissue.The pruned GRNs both linearly and non-linearly formulated are able to reach the correctpattern formed by two OC domains separated by the IOC (Fig. 3). In this interocellarregion Eya/RD values are null or very low while En concentration is high; likewise, in theocellar domains En concentration is null or very low and RD concentration is high. Theregulation of the different domains depends on the Hh gradient concentration value on eachspatial point. En concentration is dependent on CiA concentration which also depends onHh concentration as output of its signaling pathway. Therefore, the spatial point whereEn is starting self-regulating depends on Hh concentration. This behavior responds to theclassical 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 showshigh concentration in the interocellar and periocellar regions. PtcHh and CiA behaviors11re also shown as expected. They reflect the Hh gradient as being part of the Hh-signalingpathway and feel the negative effect of En repression in its domain (IOC). One must notethat the IOC domain is set in a more abrupt way in the linearly formulated model than inthe 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 [3], the simplified modelis robust against noise perturbations. This model also shows a size control of the OC whenperturbing the activation of Hth by Wg (see the details of these results in SupplementaryMaterial).
The study made to this point shows that the correct patterning seems to be extremelycorrelated to the use of a negative feedback loop with a dependence on the morphogenconcentration to trigger this inhibition. To finally prove this, a single 3-node network motifhas been extracted maintaining the negative feedback loop. The description of this motifis the following: PtcHh activates CiA, CiA activates both En and PtcHh, and En repressesthe previous two network nodes in a double negative feedback loop (Fig. 4a). The input tothis network motif is, in an adiabatic approximation, a stationary Hh gradient-like profile tosimplify the calculations, avoiding this way the PDE equation for Hh. Both the analyticalsolution for equation (27) in its stationary state or its Gaussian fit (see SupplementaryMaterial.), provide the same input to the system as in the dynamical case in a steady statesituation. Fig. 4b shows the comparison between the numerical stationary solution of theHh morphogen gradient and the Gaussian distribution fit (see details in SupplementaryMaterial). Then, the same system’s behavior is maintained. By doing this, the systembecomes a system of ODEs of the reaction type.The resulting system is formed then by three equations. ∂P tcHh∂τ = Θ ( α Hh − P tcHh Hh ( x ) + α CiA − P tcHh
CiA − α En − P tcHh En ) − β P tcHh
P tcHh (34) ∂CiA∂τ = γ PtcHh − CiA α En − CiA En α P tcHh − CiA
P tcHh n CiA k n CiA m CiA + α P tcHh − CiA
P tcHh n CiA − β CiA
CiA (35) ∂En∂τ = γ CiA − En CiA n En k n En m En + CiA n En + α En − En En Dl ( ζ En , ζ τ , τ ζ En ) − β En En (36)These equations (34-36) have the same formal expression as in eqs. (28) (PtcHh), (A.1)(CiA) and (A.2) (En) but considering a constant Hh profile in PtcHh equation (Fig. 4b).This 3-nodes ocellar model is still capable of showing qualitatively identical expressionsfor 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 noteworthyto emphasize that both a dynamic morphogen profile and its stationary solution retrievesthe correct pattern. Thus, one can provide two important results. The first one implies theuse of the dynamical description. Here, the behavior of the system is positively described at12 igure 3: Spatio-temporal representation of the signal intensity of the non-linear (left - eqs. (20-25)) andlinear (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. The 3-nodes network model raises the question of whether an ultimate simplification onthe system will lead to a correct patterning. In this case, though biologically unrealistic, itis considered that Hh is directly feeding CiA production with no intermediate; CiA is thenactivating 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 maintainedin the core of the network.Following the 3-nodes model formulation it is possible to write this system in terms ofODEs as follows: ∂CiA∂τ = γ Hh − CiA Hh ( x )1 + α En − CiA En − β CiA
CiA (37) ∂En∂τ = γ CiA − En CiA n En k n En m En + CiA n En + α En − En En Dl ( ζ En , ζ τ , τ ζ En ) − β En En (38)This simple system allows an analytical treatment to find the steady state solutions forEn and CiA. In order to ease the analysis, the steady state solutions for En and CiA in thesteady state are written in a dependent manner: CiA = 1 β CiA γ Hh − CiA Hh ( x )1 + α En − CiA En (39) En = 1 β En − α En − En Dl ( ζ En , ζ τ , τ ζ En ) γ CiA − En CiA n En k n En m En + CiA n En (40)A simple observation of eq. (40) shows that in the case of En self-regulation, this is Dl ( ζ En , ζ τ , τ ζ En ) = 1, if the parameters α En − En and β En are equal or very similar, Enconcentration goes to infinity or to a very high value; and this leads to a null or very lowCiA value. So the Hh signaling pathway is closed and no ocellar cuticle is developed. Inthe case of Dl ( ζ En , ζ τ , τ ζ En ) = 0, depending on the parameter values, positive or imaginaryvalues are obtain for both En and CiA. Though imaginary solutions can be found if n En > igure 4: (a) Network diagram of the 3-nodes model. The network dynamics is the following: Hh is theinput to the network in the form of a static Gaussian profile that activates PtcHh. PtcHh activates CiAwhile CiA activates both PtcHh and En. En represses both PtcHh and CiA. En is self-maintained in thecase Dl function is activated. Dl activates En self-regulation if En surpasses a certain concentration. (b) Hhgradient steady state in the model from the pruned GRN (in black line) and its Gaussian fit (in dashed redline). (c-e) Solutions of the 3-nodes model of 3 equations for (c) PtcHh, (d) En and (e) CiA. .5. Boolean networks analysis confirms robustness in the patterning circuit It has been demonstrated that several formulations of the same GRN with different levelsof complexity are able to determine the same patterning. This pattern seems indeed to berobust 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 parametervalues and the dynamics of a mean-field formulation, the GRN has been studied from theperspective of Boolean networks. For this, four models have been tested, the one in Fig. 1cfor the model in [3] (model i is determined by a Boolean (logical) function Φ i of itsstate and the states of those nodes that have incoming edges on it. In general, a Booleanor logical function is written as a statement acting on the inputs using the logical operators and , or and not and its output is 1(0) if the statement is true (false). For the constructionof the Φ i functions it is assumed that the inhibitors are dominant in the networks, so theywill be treated as AND logical operators [22]. This type of network modeling has alreadybeen widely used to study genetic networks [22–28]. Table 1 gives the Boolean functions Φ i for every genetic element of the four considered models.Models en , En, hth , Wg, Hth and Dl. Model eya , Eya, hth and Hthin model ci has also been set active as initial condition in16odel ci is set inactiveas initial condition, the Eya species are set to 0 as Wg is always 1 and it activates Hth beforeEya can be turned on. Therefore, the switch is kept this way in time. In model
3. Materials and Methods
The total one-dimension ocellar complex is built on a linear lattice of 117 nodes. It hasbeen measured for this work that the IOC and each ocelli is approximate 30 µm long. Theocellar complex is approximately 30 cells long and the cell diameter is ∼ µm . A goodapproximation in the model is to set each lattice node as a 1 µm of tissue length. So, itcan be assumed that each element can reach 33 nodes long, thus leaving 9 nodes free oneach side for fitting purposes. This way, a linear space of 117 nodes is enough to considerthe development of a satisfactory ocellar complex pattern. The spatial dimension on thegraphical results are, then, shown in µm units.A simplified one-node model (for the six equations (27-32) case) was implemented usingVensim software, a visual tool for solving ODEs that allows parameter values modification inrun-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 solvedwith the integrator ode45.For all the numerical experiments, the differential equations system has been discretizedusing a finite differences method.The values for all the parameters used in the simulations and the initial conditions foreach variable can be found respectively in tables 1, 2 for the non-linear models and 3 and 4for the linear models in the Supplementary Material. Some parameter values are cited andothers have been mathematically derived. The rest of the parametric set has been manuallyadjusted to obtain a similar ocellar pattern to the experimental one.The Gaussian fitting for Hh steady state profile was performed by Microcal Origin soft-ware. The Boolean networks analysis was performed using application BooleSim [29].
4. Discussion
How to establish a pattern in biological tissues is one of the most hot topics in recentdevelopmental biology, though its interest has been shown from several decades from both17
Model
Steps Hh ptc Ptc
PtcHh ci Ci A Ci R en En eya Eya hth
Hth Wg Dl Model Hh PtcHh Ci A En RD Hth Dl Model Hh PtcHh Ci A En Dl Model ! Hh Ci A En Dl (b) (a)
10 15 16
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 definingthe IOC. (b) Scenario where En is not self-regulated and the OC region is defined by CiA and Eya in anoscillatory regime. Only model
Drosophila melanogaster fly, in terms of perturbations in the sys-tem on the description level of the GRN structure representation. We started from a highdescription level PDE model using a non-linear formulation [3]. We have systematicallymodified this model by considering simpler assumptions, which keep the patterning, whilemaintaining the basic biology from a coarse-grain point of view. This study reveals thatthis pattern follows one of the most known patterning mechanisms, the
French flag model [44]. This mechanism is driven by a morphogen gradient (Hh) where the cellular fate isdetermined by the concentration of this diffusive molecule, triggering different parts of theGRN. The hh-signaling pathway is blocked by the expression of its inhibitor, En, which isalso Hh-dependent. This inhibitor becomes Hh-independent when its concentration exceedsa 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 isnot able to inhibit the expression of the RD genes and the OC region is fulfilled. This basicmechanism seems to be robust in this GRN even after being substantially pruned. We havefound that the patterning itself is independent, for example, on the receptor dynamics andon how the downstream signaling molecule CiA forms a graded profile. In this latter, wefound that ci expression could be directly or indirectly dependent on the concentration ofreceptor-bound morphogen.In order to test if the patterning is significantly affected by the type of formulationused to simulate the model, we have implemented this model using both a non-linear anda linear paradigm. This result showed that the patterning seems to be also robust againstthe choice of formulation paradigm. Comparing the results of both types of formulation onecan see that the non-linear model allows a precise control on the activation and repressionprocesses 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. Thiscontrol is less clear in the linear formulation where the parameters contribute to the rateof change of concentration by adding/subtracting concentration dependant terms. On thecontrary, 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 therepression transition is smoother in the non-linear than in the linear formulation due tothe use of a step function in the self-regulation of En in the linear case, eq. (33). Wealso found that the correct patterning is fulfilled when considering both a dynamical and asteady morphogen profile. In the dynamical description, the system can positively describedifferent developmental stages during Hh profile transition to its steady state from a nullconcentration scenario. On the other hand, as an adiabatic approximation also retrieves thecorrect pattern, we can infer that the morphogen dynamics is indeed not necessary to fulfillthe pattern. This implies that the morphogen kinetics could be faster than the intracellularresponse 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, wehave found that the patterning is dependent on the specific topology of a core regulatory net-19ork motif containing an activator-repressor regulatory mechanism. This activator-repressorparadigm is a well known mechanism in pattern formation since Turing [45]. In our caseof study, the activator-repressor mechanism is local for each cell, representing a triggeringstructure for a classical French flag model , where the activator-repressor network motif be-havior responds to the morphogen concentration input in each cell. Though mathematicallyspeaking, 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 oneintermediate to describe an extra-intracellular communication. This way, we consider theminimal 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 andCiA. This is, then, a network motif with a double negative feedback loop. It is also knownthat this type of three-component systems are appropriate to account for a wide class ofbiological phenomena [46]. Some examples are found in the enrichment of microRNA tar-gets [47, 48]. A closer example to the system studied here is the developmental patterningof neural subtype specification in mouse [49]. There, the morphogen responsible for thispatterning is the mouse Hh homologue Sonic Hedgehog (Shh). It is likely to think thatthis specific type of 3-nodes activator-repressor network motif can be responsible for thedevelopment of similar patterns in other developmental systems.We have also confirmed the patterning independence on the dynamical model implemen-tation 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 asthe core network topology is maintained. We found that the IOC region is established in asteady state fashion. Interestingly, the OC region is defined as an oscillatory regime in thefour models. Specifically, the activator-repressor mechanism is responsible for the emergenceof these oscillations. The inclusion of En as a repressor of the hh signaling-pathway estab-lishes a negative feedback within the network that might be capable of changing the systemdynamics non-linearly causing oscillations in gene expression. An experimental research forvalidating these oscillations during ocellar development has not been explored yet, thoughit is a tantalizing one. Oscillatory regimes have been proposed to increase the efficiencyof reactions [50], and to be involved in cell patterning by a morphogen graded profile [51].Even this oscillatory behavior could act as a switch-like mechanism between two transientcell-states building spatial borders among oscillating and non-oscillating cells.
Acknowledgements
The authors would like to thank Fernando Casares (CABD) for introducing them to thisfascinating system and for many fruitful discussions. This work is partially financed byJunta de Andaluc´ıa (FQM-122) to A. C´ordoba and by the Spanish Ministry for Scienceand Innovation (MICINN/MINECO) and Feder Funds through grants BFU2012-34324 toF. Casares (CABD, Seville, Spain) and Consolider From Genes to Shape, of which F. Casareswas a participant researcher. 20 ppendix A. Analysis of the cooperation reactions in CiA and En
Equations (A.1) and (A.2) show the temporal evolution of CiA and En concentrationsrespectively. ∂CiA∂τ = γ PtcHh − CiA α En − CiA En α P tcHh − CiA
P tcHh n CiA k n CiA m CiA + α P tcHh − CiA
P tcHh n CiA − β CiA
CiA (A.1) ∂En∂τ = γ CiA − En CiA n En k n En m En + CiA n En + α En − En EnDl ( ζ ) − β En En (A.2)In both equations, the production rate term is written as a Hill function, typically rep-resented as: F ( x ) = C x n k nm + x n , (A.3)where k m is the x concentration producing half F ; C is the maximum value for F and n is the Hill coefficient which describes the cooperativity or affinity to continue thereaction defined with the Hill equation once started. It is known that n < n = 1 means non-cooperative reaction and n > n CiA and n En ∈ [2 ,
4] (Fig. A.6a). Thisis, just in the defined Hill coefficients intervals En is able to achieve its self-regulation statetransiting from one to two stable states. When increasing n En over 4, the RD expressionprofile increases its intensity extending across the entire ocellar region (Fig. A.6b). Givingvalues for n CiA >
2, the RD ocellus profile sharpens acquiring a more similar to squareshape (Fig. A.6c). Hereinafter, the Hill coefficient parameters values will be n En = 4 and n CiA = 2.
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 validrange of values for k m parameters. And secondly, this analysis will show the system behaviorbefore and after the beginning of the En self-regulation, from a stability point of view. Thereare two differential equations in which a formulation of the Hill’s equation form has beenused: CiA’s (equation A.1) and En’s (equation A.2). Appendix B.1. Stability analysis for CiA’s equation
In order to achieve a stability analysis on CiA behavior we performed a simplification onequation A.1. ∂CiA∂τ = γ P tcHh − CiA α P tcHh − CiA
P tcHh n CiA k n CiA m CiA + α P tcHh − CiA
P tcHh n CiA − β CiA
CiA (B.1)21 igure A.6: (a) States plot showing the changes in En behavior for different values of the Hill coefficients n En and n CiA in equations (A.1) and (A.2). The vertical axis is the number of stable states that En can reach asa function of the Hill coefficient ( n is either n En as n CiA ). Black line represents the number of stable statesthat En reaches when modifying n En . When n En ≤
4, 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 ofstable nodes that En reaches when modifying n CiA . In this case, En shows the transition to a second statefrom n CiA = 2 on. The combination of the two Hill coefficients valid range gives a system correct behaviorfor n En ∈ [0 ,
4] and n CiA ≥ n En . For n En values higher than 4( n En >
4) En is not able to self-regulate and the RD expression patterns covers the whole ocellar complexat high level. Lower n En values show a correct RD profile lowering its intensity while decreasing n En . (c)RD profiles for different values of n CiA . Only for n CiA ≥ n CiA value increases the RD pattern sharpens x , that will represent the full output of the Hh signaling pathway. This way,equation (B.1) can be reduced to ∂x∂τ = γ P tcHh − CiA x n CiA k n CiA m CiA + x n CiA − β CiA x where x ∈ [0 , .
4] (B.2)The x value range was taken from the numerical simulations and it is built from theunion of the PtcHh and CiA value intervals. This equation can also be expressed as: ∂x∂τ = f ( x ) − β CiA x (B.3)where f ( x ) is the production rate and β CiA x is the degradation rate. Fig. B.7a showsa 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 (inblue). So, the Hh-signaling pathway can either be closed (low CiA concentration unable toregulate 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 possibleways: as a non-linear repression (equation A.1) and as a linear repression (equation B.4): ∂CiA∂τ = γ P tcHh − CiA α P tcHh − CiA
P tcHh n CiA k n CiA m CiA + α P tcHh − CiA
P tcHh n CiA − β CiA
CiA − α En − CiA En (B.4)Repeating the same procedure to these two equations: ∂CiA∂τ = γ PtcHh − CiA α En − CiA En x n CiA k n CiA m CiA + x n CiA − β CiA x = g ( x ) − β CiA x = z ( x ) (B.5)With expression B.5 for the non-linear repression form, and ∂CiA∂τ = γ P tcHh − CiA x n CiA k n CiA m CiA + x n CiA − β CiA
CiA − α En − CiA En = h ( x ) − β CiA x = z ( x ) (B.6)being the expression B.6 for the linear repression form.Two different values were chosen for En concentration depending on its state: withoutself-regulation (1.51) or with self-regulation (3502.49). These values were taken from nu-merical simulations being closed to a steady state. In order to obtain a value for k m CiA thatensures at least one stable steady state, the second derivative of z ( x ) must be negative ina { x, k m CiA } valid value pair. This value pair is the tangent point between g ( x ) and β CiA x .The first z ( x ) derivative from equation B.5 is:23 zdx = − γ PtcHh − CiA α En − CiA En x (cid:0) k m CiA + x (cid:1) + γ PtcHh − CiA α En − CiA En xk m CiA + x − β CiA (B.7)From equations (B.5) and (B.7) is obtained, as the only real and positive solution, thatthe x production and degradation rate are tangent at x = k m CiA = 2000. For these values thesystem ensures a positive stable node and a valid CiA behavior. This is, CiA concentrationis high when En is not self-regulated (black curve in Fig.B.7b), and CiA concentration isvery low or null for high concentration of En (when En is self-regulated; blue curve in Fig.B.7b). So the closure of the Hh-signaling pathway can be achieved when En represses CiAin a non-linear way. The second z ( x ) derivative is negative at the calculated x and k m CiA value, ensuring its convexity. d zdx = γ PtcHh − CiA α En − CiA En x (cid:0) k m CiA + x (cid:1) − γ PtcHh − CiA α En − CiA En x (cid:0) k m CiA + x (cid:1) + γ PtcHh − CiA α En − CiA En k m CiA + x (B.8)For k m CiA = 2000, d zdx ( x = 2000) = − · − .When evaluating the En repression action to CiA in a linear way (eq. B.6), we find thesame behavior independently of En activation. Either En being at low or high concentrationit is shown one unstable steady state at low x values and a stable steady state at high x levels (Fig. B.7c). With this linear repression the Hh-signaling pathway can remain openwith a high concentration of CiA in the cells where En is self-regulated. Moreover, as thesteady state found at low x values is unstable, perturbations applied to the system canlead its behavior to the stable, and high CiA concentration node. In other words, if thesystem is able to close the Hh-signaling pathway due to a linear En repression, it can bespontaneously open again if the system is perturbed. This behavior does not correspond tothe CiA experimentally observed in [3]. It is known that CiA reaches a high concentrationpositive 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. B.5) x shows a bi-stablestate when En is not self-regulated. One stable node at x = 0 and a second one at x = 2000.When En self-regulation starts x loses its bi-stability maintaining just the stable steady stateat null x levels (Fig. B.7b). This behavior fits with the experimental one. The only possiblestable node of CiA when En self-maintains is the closure of Hh-signaling pathway, this is,CiA concentration falls down. Appendix 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 A.2). If En is not self-regulated ( Dl = 0), this equation can be expressed as: ∂x∂τ = γ CiA − En x n En k n En m En + x n En − β En x = f ( x ) − β En x = z ( x ) (B.9)The z(x) first derivative of equation B.9 is:24 igure B.7: Stability analysis for non-linear contribution equations. (a) This graph represents the productionrate (black line) and the degradation rate (red line) of x (vertical axis: z = f ( x ), z = β CiA x in eq. B.3)as a function of the concentration of x (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. Thisconfiguration shows two different steady states (two stable nodes as blue dots). (b) Representation of thesystem steady states when En is repressing CiA in a non-linear way. Black line shows the x productionrate when En is in its lower chosen value (En = 1.51, before beginning its positive self-regulation). Blueline (very close to horizontal axis) represents the x production rate when En has reached its higher chosenvalue (En = 3502.49, after En self-regulation response). Red line represents the x degradation rate. Inthe first case the system maintains its bi-stability and in the second case the system loses the bi-stabilitymaintaining just one stable node at x = 0. (c) Representation of the system steady states when En isrepressing CiA in a linear way. Black line represents the x production rate when En remains in its lowerchosen value. Orange line represents the production rate when considering En on its higher chosen value.Red line represents the x degradation. The system shows two steady states on each configuration, one stableat high x value and one unstable at low x 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 makesthe system lose its bi-stability and does not allow any stable state at low x value to consider closing the Hhsignaling pathway. (d) Graphical determination of the first stable steady state for En evolution (equationB.9) 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 B.12). In this case, the systemloses its bi-stability remaining just one stable node at high En concentration. (e) Representation of the Enproduction rates (black line) considering the chosen value for k m En = 89 .
56, and the degradation rate (redline). The system shows a bi-stable state with null and high x levels, also considering one unstable state atlow x level zdx = − γ CiA − En x (cid:0) k m Em + x (cid:1) + 4 γ CiA − En x k m En + x − β En (B.10)The tangent point in this case is x = 7500 and k m En = 5698 . d zdx = − · − as corresponds to a maximum. d zdx = 32 γ CiA − En x (cid:0) k m En + x (cid:1) − γ CiA − En x (cid:0) k m En + x (cid:1) + 12 γ CiA − En x k m En + x (B.11)This way the system shows again a bi-stable state with a steady state at x = 0 andanother at x = 7500 (Fig. B.7d - black curve). Every value for the parameter k m En ≤ . k m En = 89 .
56, thisway the RD pattern is closer to the experimental one. By using this value, the systemmaintains its bi-stability with an intermediate unstable steady state at a very low x value(Fig. B.7e). Once En self-regulates ( Dl = 1, eq. B.12) the system loses its bi-stabilitymaintaining just one stable steady state at high x value (Fig. B.7d - blue curve). In otherwords, once En self-regulates, its concentration can only be high, and, as commented, highenough to make CiA fall down triggering the closure of the Hh-signaling pathway. ∂x∂τ = γ CiA − En x n En k n En m En + x n En + α En − En En − β En x = f ( x ) − β En x = z ( x ) (B.12) Appendix C. The ocellar regulatory network is robust facing noise
Complex networks usually include network motifs that may infer robustness to the sys-tem capable of withstanding the presence of noise and maintain its original behavior [8].This well known behavior raises the question of whether the ocelli GRN is capable to with-stand noisy perturbations and still maintain its qualitative behavior. In order to in silico demonstrate that the ocelli GRN is robust against noise we have implemented some compu-tational experiments in which a random component was inserted in the system to simulatewhite noise on the Hh signaling-pathway. These experiments applied a white noise ( ν ) (uni-form random distribution centered at one) to CiA production term (eq. C.1). The standarddeviation ( σ ) of the distribution was varied from 5% to 20% by 5%. The highest noiseproportion should be large enough to be considered superior to biological noise. ∂CiA∂τ = γ PtcHh − CiA α En − CiA En α P tcHh − CiA
P tcHh n CiA k n CiA m CiA + α P tcHh − CiA
P tcHh n CiA ν ( σ ) − β CiA
CiA (C.1)By adding this noisy signal to CiA, it is intended to simulate conditions in which vari-ations in the concentration of CiA, as the readout of the Hh-signaling pathway, may inferdifferences in the real ocellar patterning due to changes in the environmental conditions dur-ing the developmental process. The obtained results show that, even in the most dramaticsituation, the behavior of the system is qualitatively the same as in the case without noise.26 igure C.8: (a) Time series of the systems elements showing the difference of each signal when applying awhite noise (flat random distribution) to the end of the Hh-signaling pathway (CiA node) with an ampli-tude of the 20% of the original production term value. This plot shows that the system is robust againstnoise as the signals remain qualitatively, and nearly quantitatively, the same after applying the white noisedistribution. (b) RD profile showing one single ocellus pattern applying white noise and without applyingwhite noise. This plot shows that the pattern remains qualitatively the same and that the ocellus is slightlysmaller when applying the white noise distribution
Only at the maximum noise proportion, the system shows small quantitative differencesin PtcHh and Hh concentration, but never qualitative (Fig. C.8a only shown for 20% ofnoise). For this maximum noise proportion it is also observed a quantitative difference in RDconcentration, showing a narrower profile when considering CiA noisy signal (Fig. C.8b). Itwas also tested if the presence of Hth alters the noise effect in the RD profile. The resultwas that no changes appeared in the noise effect due to Hth presence, so Hth does not inferany buffer effect to RD signal.
Appendix D. Analysis of the Hh gradient steady state
It is known that hh expression is limited to the interocellar region named as hh -expressingzone. This region has a width of ω = 33 µm and it is centered in x = 0 in this model. It hasbeen described that in this region δ ( x ) = 1, so ω is the size of the morphogen source. Oncethe morphogen gradient reaches its steady state ( τ = τ std ) we have: ∂Hh∂τ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) τ = τ std = 0 (D.1)so, D ∂ Hh∂x + δ ( x ) α hh − α P tcHh − Hh P tcHh − β Hh Hh = 0 (D.2)27n the steady state, we can consider that PtcHh concentration is also constant in time, sothis can be treated as a constant value for this calculation. To solve this diffusion equationwith degradation and a source of width ω centered in x = 0 it is needed to formalizethe contour conditions. In the first place, it will be considered vanishing concentrations for x = ± L , being 2 L the whole one-dimensional field of cells that we are working on. We defined2 L = 117 µm . This assumption make sense if the characteristic length of the gradient λ ismuch smaller than L . Secondly, it will be considered that both the morphogen concentration, Hh ( x ), and its flux, j = − D ∂Hh∂x + β Hh Hh must be continuous at the boundaries of themorphogen source ( x = ± ω/ Hh ( x ) = Hh b e ( x + ω ) λ − − S b if x ≤ − ω α hh β Hh (cid:16) − e ( x − ω ) λ − − e ( − x − ω ) λ − (cid:17) − S b if x − ω ≤ x ≤ ω Hh b e ( − x + ω ) λ − − S b if x ≥ ω (D.3)With Hh b = α hh β Hh (cid:0) − e − ωλ (cid:1) and S b = α PtcHh − Hh P tcHhβ Hh , being Hh b − S b the Hh concen-tration at the boundaries of the source. λ = (cid:112) D/β Hh is also known as decay length. Thislast corresponds to the distance at which the morphogen concentration decays by a factorof 1 /e . It can be calculated from the parameters used in the simulation (see Table D.2) thatthe value for λ = 0 . µm , so λ << L . Appendix D.1. Hh morphogen gradient fits a Gaussian distribution
It has been shown that the analytical solution for the Hh morphogen equation on itsstationary state can be written as a combination of exponential functions (equation D.3).From an experimental perspective, it is known that morphogen gradients profiles, whichspatial distribution is similar to a 1D profile with a non-punctual source, can follow aGaussian distribution or a sum of Gaussian distributions. An example of this is the Dppmorphogen gradient in the Drosophila wing disc [32] (see Appendix D). The same way, we cantreat the numerical output as an experimental profile and fit it to a Gaussian distribution.This distribution expression is the following: Hh ( x ) = A + Aσ √ π e − ( x − xcσ ) (D.4)In Fig. 4a is shown the comparison between the numerical stationary solution of theHh morphogen gradient and the Gaussian distribution fitting with A = − ± A =6 . · ± · , σ = 23 . ± .
15 and x c = 0. The errors in these parameters come fromthe fitting calculation with a R-squared value of R = 0 . References [1] J. E. Sulston, H. R. Horvitz, Post-embryonic cell lineages of the nematode, Caenorhabditis elegans.,Developmental Biology 56 (1) (1977) 110–156.
2] E. Dessaud, V. Ribes, N. Balaskas, L. L. Yang, A. Pierani, A. Kicheva, B. G. Novitch, J. Briscoe,N. Sasai, Dynamic assignment and maintenance of positional identity in the ventral neural tube by themorphogen sonic hedgehog., PLoS Biology 8 (6) (2010) e1000382.[3] D. Aguilar-Hidalgo, M. A. Dom´ınguez-Cejudo, G. Amore, A. Brockmann, M. C. Lemos, A. C´ordoba,F. Casares, A Hh-driven gene network controls specification, pattern and size of the Drosophila simpleeyes., Development 140 (1) (2013) 82–92.[4] S. Kondo, T. Miura, Reaction-diffusion model as a framework for understanding biological patternformation., Science 329 (5999) (2010) 1616–1620.[5] H.-W. Kang, L. Zheng, H. G. Othmer, The effect of the signalling scheme on the robustness of patternformation in development., Interface focus 2 (4) (2012) 465–486.[6] F. Li, T. Long, Y. Lu, Q. Ouyang, C. Tang, The yeast cell-cycle network is robustly designed., Proceed-ings of the National Academy of Sciences of the United States of America 101 (14) (2004) 4781–4786.[7] T. Schlitt, A. Brazma, Modelling gene networks at different organisational levels., FEBS letters 579 (8)(2005) 1859–1866.[8] M. Nahmad, A. Stathopoulos, Dynamic interpretation of hedgehog signaling in the Drosophila wingdisc., PLoS Biol. 7 (9) (2009) e1000202.[9] J. Royet, R. Finkelstein, hedgehog, wingless and orthodenticle specify adult head development inDrosophila., Development 122 (6) (1996) 1849–1858.[10] N. M. Bonini, W. M. Leiserson, S. Benzer, The eyes absent gene: genetic control of cell survival anddifferentiation in the developing Drosophila eye., Cell 72 (3) (1993) 379–395.[11] B. N. Cheyette, P. J. Green, K. Martin, H. Garren, V. Hartenstein, S. L. Zipursky, The Drosophilasine oculis locus encodes a homeodomain-containing protein required for the development of the entirevisual system., Neuron 12 (5) (1994) 977–996.[12] M. A. Serikaku, J. E. O’Tousa, sine oculis is a homeobox gene required for Drosophila visual systemdevelopment., Genetics 138 (4) (1994) 1137–1150.[13] N. M. Bonini, W. M. Leiserson, S. Benzer, Multiple roles of the eyes absent gene in Drosophila., Dev.Biol. 196 (1) (1998) 42–57.[14] J. Blanco, M. Seimiya, T. Pauli, H. Reichert, W. J. Gehring, Wingless and Hedgehog signaling pathwaysregulate orthodenticle and eyes absent during ocelli development in Drosophila, Dev. Biol. 329 (1) (2009)104–115.[15] J. Blanco, T. Pauli, M. Seimiya, G. Udolph, W. J. Gehring, Genetic interactions of eyes absent, twinof eyeless and orthodenticle regulate sine oculis expression during ocellar development in Drosophila.,Dev. Biol. 344 (2) (2010) 1088–1099.[16] A. Brockmann, M. A. Dom´ınguez-Cejudo, G. Amore, F. Casares, Regulation of ocellar specificationand size by twin of eyeless and homothorax, Dev. Dyn. 240 (1) (2010) 75–85.[17] Y. Chen, G. Struhl, Dual roles for patched in sequestering and transducing Hedgehog., Cell 87 (3)(1996) 553–563.[18] V. Marigo, C. J. Tabin, Regulation of patched by sonic hedgehog in the developing neural tube., Proc.Nat. Acad. Sci. USA 93 (18) (1996) 9346–9351.[19] D. M. Umulis, Analysis of dynamic morphogen scale invariance., Journal of the Royal Society, Interface/ the Royal Society 6 (41) (2009) 1179–1191.[20] S. Zhou, W.-C. Lo, J. L. Suhalim, M. A. Digman, E. Gratton, Q. Nie, A. D. Lander, Free extracellulardiffusion creates the Dpp morphogen gradient of the Drosophila wing disc., Current biology : CB 22 (8)(2012) 668–675.[21] C. Alexandre, A. Baena-Lopez, J.-P. Vincent, Patterning and growth control by membrane-tetheredWingless., Nature 505 (7482) (2014) 180–185.[22] R. Albert, H. G. Othmer, The topology of the regulatory interactions predicts the expression pattern ofthe segment polarity genes in Drosophila melanogaster, Journal of Theoretical Biology 223 (1) (2003)1–18.[23] S. A. Kauffman, Metabolic stability and epigenesis in randomly constructed genetic nets., Journal ofTheoretical Biology 22 (3) (1969) 437–467.
24] R. Thomas, Boolean formalization of genetic control circuits., Journal of Theoretical Biology 42 (3)(1973) 563–585.[25] R. Thomas, R. D’Ari, Biological Feedback, CRC Press, 1990.[26] S. A. Kauffman, The Origins of Order, Self Organization and Selection in Evolution, Oxford UniversityPress, 1993.[27] A. Gonz´alez, C. Chaouiya, D. Thieffry, Dynamical analysis of the regulatory network defining thedorsal-ventral boundary of the Drosophila wing imaginal disc., Genetics 174 (3) (2006) 1625–1634.[28] J. Buceta, H. Herranz, O. Canela-Xandri, R. Reigada, F. Sagu´es, M. Mil´an, Robustness and Stabilityof the Gene Regulatory Network Involved in DV Boundary Formation in the Drosophila Wing, PLoSONE 2 (7) (2007) e602.[29] M. Bock, T. Scharp, C. Talnikar, E. Klipp, BooleSim: an interactive Boolean network simulator.,Bioinformatics (Oxford, England) 30 (1) (2014) 131–132.[30] L. Wolpert, Positional information and pattern formation in development, Developmental Genetics15 (6) (1994) 485–490.[31] L. Lum, The Hedgehog Response Network: Sensors, Switches, and Routers, Science 304 (5678) (2004)1755–1759.[32] A. Kicheva, P. Pantazis, T. Bollenbach, Y. Kalaidzidis, T. Bittig, F. J¨ulicher, M. Gonz´alez-Gait´an,Kinetics of morphogen gradient formation., Science 315 (5811) (2007) 521–525.[33] O. Wartlick, P. Mumcu, A. Kicheva, T. Bittig, C. Seum, F. Julicher, M. Gonzalez-Gaitan, Dynamicsof Dpp Signaling and Proliferation Control, Science 331 (6021) (2011) 1154–1159.[34] S. Restrepo, J. J. Zartman, K. Basler, Coordination of Patterning and Growth by the Morphogen DPP.,Current biology : CB 24 (6) (2014) R245–R255.[35] J. R. Collier, N. A. Monk, P. K. Maini, J. H. Lewis, Pattern formation by lateral inhibition withfeedback: a mathematical model of delta-notch intercellular signalling., Journal of Theoretical Biology183 (4) (1996) 429–446.[36] H. Meinhardt, Models of biological pattern formation: common mechanism in plant and animal devel-opment, Int. J. Dev. Biol. (1996) 123–134.[37] N. Barkai, S. Leibler, Robustness in simple biochemical networks., Nature 387 (6636) (1997) 913–917.[38] D. Ben-Zvi, N. Barkai, Scaling of morphogen gradients by an expansion-repression integral feedbackcontrol., Proceedings of the National Academy of Sciences of the United States of America 107 (15)(2010) 6924–6929.[39] M. Nahmad, A. Stathopoulos, Establishing positional information through gradient dynamics: a lessonfrom the Hedgehog signaling pathway, Fly 4 (4) (2010) 273–277.[40] Z. Burda, A. Krzywicki, O. C. Martin, M. Zagorski, Motifs emerge from function in model gene reg-ulatory networks., Proceedings of the National Academy of Sciences of the United States of America108 (42) (2011) 17263–17268.[41] M. Nahmad, Steady-state invariant genetics: probing the role of morphogen gradient dynamics indevelopmental patterning., Journal of the Royal Society, Interface / the Royal Society 8 (63) (2011)1429–1439.[42] A. Kicheva, M. Cohen, J. Briscoe, Developmental pattern formation: insights from physics and biology.,Science 338 (6104) (2012) 210–212.[43] W. Li, H. Hu, Y. Huang, H. Li, M. R. Mehan, J. Nunez-Iglesias, M. Xu, X. Yan, J. Zhou, Xianghong,Frequent pattern discovery in multiple biological networks: Patterns and algorithms, Stat. Biosci.(2012) 157–176.[44] L. Wolpert, The French Flag problem: a contribution to the discussion on pattern development and reg-ulation. In: Waddington CH, ed. Towards a Theoretical Biology: Prolegomena, Edinburgh UniversityPress, pp. 125-133, 1968.[45] A. M. Turing, The Chemical Basis of Morphogenesis, Philosophical Transactions of the Royal Societyof London 237 (1952) 36.[46] H. Meinhardt, Out-of-phase oscillations and traveling waves with unusual properties: the use of three-component systems in biology, Physica D (2004) 264–277.
47] Q. Cui, Z. Yu, E. O. Purisima, E. Wang, Principles of microRNA regulation of a human cellularsignaling network., Molecular Systems Biology 2 (2006) 46.[48] C. Marr, F. J. Theis, L. S. Liebovitch, M.-T. H¨utt, Patterns of Subnet Usage Reveal Distinct Scales ofRegulation in the Transcriptional Regulatory Network of Escherichia coli, PLoS Computational Biology6 (7) (2010) e1000836.[49] N. Balaskas, A. Ribeiro, J. Panovska, E. Dessaud, N. Sasai, K. M. Page, J. Briscoe, V. Ribes, GeneRegulatory Logic for Reading the Sonic Hedgehog Signaling Gradient in the Vertebrate Neural Tube,Cell 148 (1-2) (2012) 273–284.[50] P. H. Richter, J. Ross, Oscillations and efficiency in glycolysis., Biophysical chemistry 12 (3-4) (1980)285–297.[51] A. Koch, H. Meinhardt, Biological pattern formation: from basic mechanisms to complex structures,Rev. Mod. Phys. 66 (1994) 1481. able 1: Boolean functions Node Boolean updating function Φ i Model Hh i Hh i = 1 for all tptc i ptc t +1 i = CiA i AND NOT
CiR i AND NOT En i P tc i P tc t +1 i = ptc i P tcHh i P tcHh t +1 i = Hh i AND
P tc i ci i ci t +1 i =NOT En i CiA i CiA t +1 i = ( P tcHh i AND NOT
P tc i ) OR ci i CiR i CiR t +1 i = CiA i AND (NOT
P tcHh i AND
P tc i ) en i en t +1 i = ( CiA i AND NOT
CiR i ) OR ( En i AND Dl i ) En i En t +1 i = en i eya i eya t +1 i = Eya i OR (
CiA i AND NOT
Hth i ) Eya i Eya t +1 i = eya i hth i hth t +1 i = W g i AND NOT
Eya i Hth i Hth t +1 i = hth i W g i W g i = 1 for all tDl i Dl i = (cid:26) i ∈ IOC domain0 otherwise for all t Model Hh i Hh i = 1 for all tP tcHh i P tcHh t +1 i = ( Hh i OR CiA i ) AND NOT En i CiA i CiA t +1 i = P tcHh i AND NOT En i En i En t +1 i = CiA i OR ( En i AND Dl i ) RD i RD t +1 i = ( CiA i OR RD i ) AND NOT Hth i Hth i Hth t +1 i = NOT RD i Dl i Dl i = (cid:26) i ∈ IOC domain0 otherwise for all t Model Hh i Hh i = 1 for all tP tcHh i P tcHh t +1 i = ( Hh i OR CiA i ) AND NOT En i CiA i CiA t +1 i = P tcHh i AND NOT En i En i En t +1 i = CiA i OR ( En i AND Dl i ) Dl i Dl i = (cid:26) i ∈ IOC domain0 otherwise for all t Model Hh i Hh i = 1 for all tCiA i CiA t +1 i = Hh AND NOT En i En i En t +1 i = CiA i OR ( En i AND Dl i ) Dl i Dl i = (cid:26) i ∈ IOC domain0 otherwise for all t able D.2: Parameter values Parameter Description Value D Hh diffusion constant 0 . µm s − , [8] α hh Hh transcription rate 10 µM s − α Hh − P tcHh
Association rate for PtcHh complex formation 1 s − α P tcHh − Hh Ptc efficiency to capture Hh 5 · − s − α P tcHh − CiA
PtcHh activating CiA rate 1 α CiA − P tcHh
CiA activating PtcHh rate 0 . s − α CiA − RD CiA activating RD rate 1 s − α En − P tcHh
En repressing PtcHh rate 1 s − α En − CiA
En repressing CiA rate 0 . α En − En En translation rate 1 s − α RD − Hth
RD repressing Hth rate 0 . s − α Hth − RD Hth repressing RD rate 0 . s − α RD − RD RD auto-regulation rate 0 . s − α W g − Hth
Wg activating Hth rate 4 . µM s − Dl En self-activation switch 1-ON, 0-OFF ζ En En threshold for Dl activation 256 µMζ t Time interval with En over ζ En for Dl activation 1400s γ P tcHh − CiA
Maximum CiA production rate 10000 µM s − γ CiA − En Maximum En activation by CiA 10000 µM s − β Hh Hh degradation rate 1 s − β P tcHh
PtcHh degradation rate 1 s − β CiA
CiA degradation rate 1 s − β En En degradation rate 1 s − β RD RD degradation rate 1 s − β Hth
Hth degradation rate 0 . s − k m CiA
CiA half-maximal activation concentration 2000 µM (sec. Appendix B.1) k m En En half-maximal transcriptional activation 89 . µM (sec. Appendix B.2) n CiA
Hill coefficient for CiA transcription 2 (sec.Appendix A) n En Hill coefficient for En transcription 4 (sec. Appendix A)33 able D.3: Variables initial conditions
Variable Description Initial conditionHh Hh concentration 0PtcHh PtcHh complex concentration 100 µM CiA CiA concentration 0En En concentration 1 µM RD RD concentration 1 µM Hth Hth concentration 156 . µM aa Initial Hth concentration corresponds to its stationary concentration in the absence of RDrepression ([