Analysis of the transcriptional logic governing differential spatial expression in Hh target genes

RESEARCH ARTICLE Analysis of the transcriptional logic governing differential spatial expression in Hh target genes Manuel Cambo ´n 1,2☯ , O ´scar Sa ´nchezID 1,2☯ * 1Applied Mathematics Department, University of Granada, Granada, Spain, 2Excellence Research Unit “Modeling Nature” (MNat), University of Granada, Granada, Spain ☯These authors contributed equally to this work. *[email protected] Abstract This work provides theoretical tools to analyse the transcriptional effects of certain biochemical mechanisms (i.e. affinity and cooperativity) that have been proposed in previous literature to explain the proper spatial expression of Hedgehog target genes involved in Drosophila development. Specifically we have focused on the expression of decapentaplegic, wingless,stripe and patched. The transcription of these genes is believed to be controlled by enhancer modules able to interpret opposing gradients of the activator and repressor forms of the transcription factor Cubitus interruptus (Ci). This study is based on a thermodynamic approach, which provides expression rates for these genes. These expression rates are controlled by transcription factors which are competing and cooperating for common binding sites. We have made mathematical representations of the different expression rates which depend on multiple factors and variables. The expressions obtained with the model have been refined to produce simpler equivalent formulae which allow for their mathematical analysis. Thanks to this, we can evaluate the correlation between the different interactions involved in transcription and the biological features observed at tissular level. These mathematical models can be applied to other morphogenes to help understand the complex transcriptional logic of opposing activator and repressor gradients. Introduction Hedgehog (Hh) is a morphogen, a signalling protein that induces several cellular responses. It is involved in the development of different biological systems, for example that of, the Drosophila melanogaster fly. In Drosophila’s wing imaginal disc the secretion of Hh from the Posterior compartment cells induces the expression of several target genes inside the cells in the Anterior compartment. Among them are decapentaplegic (dpp) and patched (ptc). Both give rise to the synthesis of their corresponding proteins, Dpp and Ptc, which are essential for the wing central domain development [1,2]. In the embryonic ectoderm Hh also regulates wingless (wg) and stripe (sr) genes. PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 1 / 25 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Cambo ´n M, Sa ´nchez O ´(2019) Analysis of the transcriptional logic governing differential spatial expression in Hh target genes. PLoS ONE 14(1): e0209349. https://doi.org/10.1371/journal. pone.0209349 Editor: Suzannah Rutherford, Fred Hutchinson Cancer Research Center, UNITED STATES Received: June 5, 2018 Accepted: December 4, 2018 Published: January 7, 2019 Copyright: ©2019 Cambo ´n, Sa ´nchez. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All relevant data are within the paper and its Supporting Information files. Funding: This work has been partially supported by: the Ministerio de Economı ´a y competitividad (Spain Government) research grants numbers MTM2014-53406-R and FPI2015/074837 http:// www.mineco.gob.es/portal/site/mineco/ and the Consejerı ´a de Economı ´a, Innovacio ´n, Ciencia y Empleo, Junta de Andalucı ´a (Andalucı ´a Government) Project FQM 954, https://www. juntadeandalucia.es/organismos/ However, it is known that the same signal of Hh produces different spatial expression of theses genes. That is to say, the expression of ptc is only limited to disc zones close to the Anterior/Posterior (A/P) border with high Hh concentrations, while dpp expresses in a broader disc range under low Hh concentrations. This poses a question: Why does the same signal give rise to different spatial expressions for different genes? The answer to this question is still under debate. The current understanding is that both genes respond, basically, to the same principles that we list below. Hh transcriptionally controls both Dpp and Ptc through the transcription factors (TFs) Cubitus interruptus (Ci). It dictates the activity of RNA polymerase enzymes (RNAP), which controls the genetic transcription via the synthesis of Ribonucleic Acid (mRNA). This process requires the binding of RNAP to some specific sites on the DNA chain called promoters. However, the transcription rate of the target genes not only depends on the total concentration of RNAP in the system, but also is controlled by the protein Ci. Ci is present in two opposite forms: activator and repressor. The activators, CiA, attempt to promote the transcription rate while the repressors, CiR, attempt to decrease it. Hh signal affects the balance between both forms of Ci, i.e., in the absence of Hh Ci appears in its repressed form but when Hh is absorbed by the cell, Cubitus changes it role presenting its activator form. So, the Hh gradient in the Anterior compartment creates opposing activator (CiA) and repressor (CiR) gradients. Furthermore both Ci forms need to bind specific DNA sites called enhancers or cis-regulatory sites which are different from the RNAP binding sites. A single promoter can be regulated by one or many relatively short enhancer modules, which are activated/repressed by binding of multiple TFs. Thus the expression pattern of a gene reflects the combined activity of all the enhancer modules that are capable of activating/repressing its transcription [3,4]. In this work we plan to analyse the control executed by a particular type of modules denominated Hh/Ci target enhancers that integrate competing inputs with opposing transcriptional functions [5]. Only a limited number of direct Hh/Ci target enhancers have been identified in Drosophila over the past quarter century in Hh target genes as dpp,ptc,wg or sr, see details and more examples in Table 1 in [5], [6] and references therein. Among them, we have that dpp, in Drosophila imaginal disc, is both activated and repressed by CiA and CiR that are in constant competition for the binding of a module of 3 enhancers (dppD) [7,8]. We have to remark that the control of these Hh/Ci modules does not probably justify the absolute genetic expression of the target gene although they have been described to be essential for their proper genetic spatial patterns [5]. Furthermore, the methodology and results provided in this work allow us to understand the contribution of a single Hh/Ci module even in the lack of knowledge about the control executed by other transcription factors or other enhancer modules. Some recent works [5,9–11] postulate that the reason for the proper spatial expression of these genes could be found in certain biochemical factors involved in the transcription process. Firstly, the binding of both RNAP and Ci in the promoter and enhancers is carried out by chemical reactions. These require some free energy that is commonly characterised by a magnitude called binding affinity. This affinity depends on several characteristics of the promoters and enhancers of each transcribed gene. In fact, in [5] it was observed that the enhancers with lower relative affinity seem to be necessary to obtain normal expression of dpp in regions of low signal. Secondly, it is possible that transcription factors that are already bound in some enhancers can modify the affinity of other binding elements. In this case, bound TFs may modify the free energy of a later binding reaction of either TF or RNAP. This process is generally termed cooperativity, however this can be positive or negative. If it facilitates the binding it is called (normal) cooperativity and if it impedes it, is called anti-cooperativity. In [12,13] it was proposed that the activator/repressor TFs modify the transcription rate by promoting or Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 2 / 25 empleoempresaycomercio.html. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist. blocking respectively the recruitment of RNA polymerase. This implies that cooperativity or anti-cooperativity with the RNAP changes the promoter binding affinity. The combination of all these biochemical factors (competition, cooperativity and binding affinities) gives rise to a very complex balance between the concentration of activators and repressors making it difficult to discern their interacting effects at tissular level. In [5,10], the spatial expression of some of the Hh target genes was related to the respective binding affinity between Cubitus proteins and Hh/Ci module enhancers. The relative in vitro affinities of Ci sites in the ptc and dpp enhancers have been measured by electrophoretic mobility shift assays, see Fig S2 in [10]. In Fig 4A in [11], the Ci binding affinity of four Ci sites in the wg embryonic ectoderm enhancer was measured by using the same methodology. In Table 1 in [5], a Ci binding sites rank of 9-mer in order of predicted binding affinity for Ci for genes in the Drosophila genome for which Hh/Ci-regulated enhancers have been functionally characterised. It has been observed that ptc is activated by Hh/Ci in larval imaginal discs via a module with high-affinity Ci sites, by contrast, with the relative low-affinity of dpp,wg and sr enhancers located in their corresponding Hh/Ci modules. The experiments developed in [10] confirm that, under moderate Hh signal, the wild type low-affinity sites in dppD produce activation, whereas if they are substituted by high-affinity sites produce repression. Similar results were obtained in embryonic enhancers of wg and sr in [5,11]. Let us remark that this activation/repression is relative to the transcription levels observed in the absence of any enhancer in the Hh/Ci module (Fig 2 in [10], Supplementary Fig 5 in [11], Fig 2E in [5] and [14]). It is clear the existence of other mechanisms for the Hh gradient interpretation in the rest of the locus. In the special case of dpp the expression levels in the absence of these three binding sites in dppD show a deep repression close to the A/P boundary. To the author’s knowledge, this repression is not fully understood (see Results and discussion (h) in [5]). Nevertheless, the contribution (activation/repression) of the single enhancer module dppD can be perceived by the analysis of relative expression levels. In order to discriminate between the mechanisms that could give rise to a differential spatial expression, in [10,11] the experiments were contrasted with fittings to a thermodynamical model based on the ideas of Shea, Ackers and coworkers [15,16]. Furthermore, by fitting a repressor cooperativity model in [10] they observed that CiR plays a substantial role in the response to a moderate Hh signal. In [11] the authors also proposed that the cooperativity between repressors may play an important role in the change of the genetic expression along the imaginal disc, by using a mathematical model of occupancy competition between repressors and activators. The large amount of biochemical variables that are present in the system calls for mathematical models [17,18] that can shed some light on the origins of the differential spatial expression in the target genes of Hh, among others. The thermodynamic model proposed by Shea, Ackers and coworkers [15,16], also known as BEWARE [19] (Binding Equilibrium Weighted Average Rate Expression), is a method frequently used in the mathematical modelling of genetic transcription processes. See [20] or [21] for a general discussion/comparison with other modelling approaches as for instance Boolean models. However, this model gives rise to long and complex mathematical expressions even when there are only a few transcription factors involved. For the analysis of independent and specific binding sites and the analysis of two non competitive transcription factors, only simple mathematical expressions have previously been proposed [22,23]. It is difficult to decipher the biological effects in the model even if they are supported by numerical tools [24] because the expressions inherently involve a great number of constants and variables. In this work we try to have a better understanding of the transcriptional logic of target genes controlled by a Hh/Ci module of enhancers, from a theoretical point of view by using a Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 3 / 25 thermodynamic model. Our analysis proposes that the transcriptional logic in the presence of opposing activator/repressor gradients can exhibit different versions depending on the cooperativity between the transcription factors. In fact, the theoretical methodology we developed is able to demonstrate how the combination of similar biochemical factors applied under different frameworks gives rise to completely different transcriptional effects. We have tested this general framework with experimental results for several Hh target genes, such as dpp,ptc,wg and sr. Among all the different frameworks deduced, we conclude that proper spatial expression of these genes is due to differences in affinities between their corresponding binding sites in combination with cooperativity between repressors (please refer to Section Results for more details). Our analysis confirms the results obtained in [4,5,10,11]. However, this does not exclude that the control executed by the Hh/Ci module can be modulated by other enhancer modules or signals depending on the biological context. Recently it has been described in [4] that the Hh/Ci module activity can be stage-tissue specific which could be related with the synergistic work of different modules. In the case of a single activator gradient our analysis represents the well accepted transcriptional logic which comes from the activator threshold model. This model explains the role of certain biochemical factors involved in the signalling interpretation. For example, differential affinities of activators for DNA elements [25] and cooperativity between activators [10,26]. High-affinity binding sites and cooperativity between activators benefit the binding of the activators to the enhancers, allowing the expression of genes at low activator concentrations, so here we observe a broader response within the activator gradient. In contrast, low-affinity sites and the absence of cooperativity between activators restrict the gene expression to high activator concentration regions. Although this rationale is well accepted in the single gradient scenario, it has not been succesfully applied in combinatorial interactions as, for instance, in opposing activator-repressor gradients such as those we have described in Hh signalling [10, 25]. In this case the balance between both gradients causes the existence of ranges of net activated/repressed cells, i.e., cells along the tissue that express higher/lower levels than basal level [10], and hence there is no global activation/repression. The cellular expression ranges (CERs) will not be determined only by signal intensity but also by net activated/repressed cellular ranges. Our analysis suggests that biochemical differences between genes can affect both signal modulation and changes in the net activated cellular ranges. That is, the variation of CERs can be explained by analysing the combination of both aspects. A very remarkable point in this work is that using the same biochemical characteristics with different cooperativity between activators and repressors will give drastically different expression rates. Different types of cooperative interactions between transcription factors will produce variations in the transcription logic. Another very interesting aspect in our analysis is how the number of enhancers could modify the gene expression. In [10] a transgenic fly line carrying a GFP reporter with a single high-affinity Ci site was constructed in order to detect whether cooperativity between TFs played a important role in the signalling process. Even in absence of cooperativity, this variation of the number of enhancers could modify the gene expression, at least theoretically. This aspect, will play a key role in our argument in order to understand Hh target genes. Although they seem to play a central role, affinity and cooperativity between TFs are not the only biochemical factors involved in the interpretation of general morphogen signalling. In the development of the chick/mouse embryo neural tube, another paradigmatic morphogenetic patterning example, cells are differentiated in response to Sonic Hedgehog (Shh) morphogenetic signals [27–30]. In this case, the Shh signal balances the concentration of different versions of activators and repressors of the Gli family. These Gli TFs recognise target sequences which are very similar, however it has been proposed that some TFs have more Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 4 / 25 potency as activators or repressors [26] than others. Thus, the different potency as activators and repressors could be another factor to be taken into account by using differential TFsRNAPs cooperativity. Results In this work we will study which of the previous factors can vary cellular expression ranges by analysing their combined effect on both, the modulation of the signal and the variation of the ranges of net activated/repressed cells. We begin with a modelling exercise that will provide the mathematical representation (operators) of expression levels to be later analysed. Since we are interested in the transcriptional effect of an Hh/Ci module we apply the BEWARE method in order to obtain expression rates of a gene controlled by two opposing general transcription factors: the activator CiA and the repressor CiR. The “BEWARE operators” obtained include, for each gen: • the transcription factors competing in order to bind one of ncommon cis-regulatory sites, • separate binding affinities of the TFs depending on their activator or repressor form, • cooperativity interactions between TFs and cooperativity/anti-cooperativity between TFs and RNAPs. We consider different process of TF-TF cooperativity that give rise to different expressions of the BEWARE operator: • non-cooperativity, • total cooperativity, when a bound TF modifies the binding affinity of any other transcription factor, • partial cooperativity, that only takes place between TFs that are of the same form, activators only interact with activators and repressors only interact with repressors. Panel A in Fig 1 shows all these interactions with n= 3 binding sites. We have simplified the mathematical expressions of the BEWARE operators because in their original form it would be impossible to make the subsequent mathematical analysis. From the analysis of these mathematical expressions, we can estimate which biochemical characteristics can modify the spatial expression of two genes controlled by the same TFs and in what way. The experimental evidences that motivate our analysis are mainly related with the affinity and the number nof the binding sites. By electrophoretic mobility shift assays Parker and coauthors found in [10] that Ci binding sites in the ptc enhancer have considerably higher affinity than dpp sites. The same authors constructed transgenic fly lines that allow them to compare the transcriptional activity of reporter genes containing different variants of these sites modifying their affinity. Similar results were described in embryo for wg in [11]. Since there are two opposing signals we have to determine when a cell is net activated or repressed. These notions have been adopted from [5,10] where using reporter genes, the activity of different versions of the dpp enhancer containing three low-affinity sites (dppD-Ci WT ), three high-affinity sites (dppD-3xCi ptc ) or three null-affinity sites (dppD-3xCi KO ) were compared. The reporter gene with null-affinity sites provided the basal expression, since it reflects the effects of all other factors which are different than Ci on the module dppD. We refer to Results and discussion (h) in [5] for a more detailed discussion. Specifically, in [5,10] the effects of Ci signalling with lowor high-affinity enhancers was measured comparing the gene activity versus the basal in any cell. Cells expressing a gene with higher expression rates than the basal level are called net activated cells. The set of all the net activated cells constitutes the Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 5 / 25 net activated cellular (NAC) range. Fig 1 shows how this range is determined by using a thermodynamic model in the same way as was done from measurements in [10]. We want to find out which concentrations of activators and repressors, [CiA] and [CiR], will provide more or less gene expression than the basal. So, we define a threshold separating concentrations of both TFs that would produce net activated cells or net repressed cells using the BEWARE operators. Then, we can use the threshold between net activation/repression concentrations to determine the limit between the ranges of net activated or repressed cells. Once we have defined the ranges of activated cells we can predict how biochemical differences will affect them as well as the signal intensity. To do this analysis we need to assume that the opposing activator and repressor gradients are monotone and do not change over time (see Panel B in Fig 1 for a graphical example). Let us mention that the approach we follow is independent of the specific values adopted by the TFs concentrations. Mathematical analysis is able Fig 1. Net activated cellular (NAC) range described by a thermodynamic model. A) Schematic of the experiment for NAC range determination. The arrows represent all the possible interactions captured by a thermodynamic model determining the transcription rates: double-headed straight arrows show protein-DNA binding site affinities while single-headed black and red arrows are TFs-RNAP and TFs-TFs cooperativities respectively. The net activated cellular range of dppD3xCi ptc , a reporter gene with a version of the dpp enhancer with three high-affinity binding sites, is obtained by comparing its theoretical transcriptional activity with the activity of dppD3xCi KO , a gen containing different version of the dpp enhancer containing three null-affinity sites. Both cases are represented in the upper and lower schemes respectively. TFs binding sites are represented by rounded rectangles filled in green (high-affinity) or black (null affinity). B) Theoretical transcription rates predicted for both genes in cells of the Anterior compartment. This compartment occupies the 60% of the Drosophila imaginal disc and the Posterior compartment the rest (60% to 100%). The expression levels given by the BEWARE operators are between 0nM/min and 1nM/min being the basal level equal to 0.5nM/min. These reference expression levels have been chosen for a proper appreciation of signal modulation. Since dppD3xCi KO has been modelled independent of external factors it is expressed at basal level anywhere. Cells expressing dppD3xCi ptc more than the basal level are in the NAC range. The expression of both genes in the wing imaginal disc is also indicated by using coloured bars. The blue circle inside the bar, indicate the position of a cell expressing dppD3xCi ptc at the basal level. The color scale used in these bars is shown in C) black meaning no expression (0nM/min), and full color meaning high expression (1nM/min). The inset in B) depicts the activator/repressors (CiA/CiR) gradients generated by Hh signalling: activator concentrations are higher close to the Anterior/Posterior border. A more detailed description can be found in Eq (16). https://doi.org/10.1371/journal.pone.0209349.g001 Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 6 / 25 to detect how biochemical differences provoke variations in transcription rates by only assuming that they respond to the same (unknown) opposing activator/repressor concentrations. This is in contrast with the great variability in TFs gradients determination exhibited by the theoretical models fitted in [10]. This methodology is explained in detail in Section Methods. Note that another theoretical approach has been proposed in the previous work [11]. With the help of the analysis performed we have determined which biochemical factors could be involved in the differential expression observed in Hh target genes. To do this we test and compare our theoretical results with the already existent experimental evidence, in particular, for dpp,wg,sr and ptc. In all cases we firstly describe the experiments and results observed. On a second stage we contrast these measurements with the different qualitative behaviours predicted by our modelling and finally we deduce the biochemical framework that explains the experimental evidences. We will first apply this methodology to dpp and the dppD enhancers module. In this case there exist several independent experimental approaches predicting the same biochemical framework which could be interpreted as a validation of our theoretical analysis. Experiment 1: Transcriptional effects of the reduction binding sites In [10], a transgenic fly line carrying a GFP reporter with a single high-affinity Ci site (dppD1xCi ptc , where the superscript ‘ptc’ stands for the high-affinity version of the dpp enhancers following the original nomenclature) was also constructed. In this case, the range of net activated cells was wider than the range of net activated cells for (dppD-3xCi ptc ) and in consequence a broader, but attenuated, expression was observed for the single enhancer case than in the 3 high-affinity sites case. This comparison can be reproduced by using the BEWARE operators. It can be theoretically proved that the cooperativity between the TFs determines the effect of the binding sites reduction: • In the presence of total cooperativity or non-cooperativity between the TFs the NAC range would essentially remain unaltered although a reduction of signal intensity could be observed, that is, less repression/activation in the repressed/activated cells. • If the activators CiA only cooperate between them a reduction in the NAC range would be observed, so some net activated cells for dppD-3xCi ptc would change to be net repressed for the gene dppD-1xCi ptc . • Finally, in the case of cooperativity only between repressors the NAC range would be incremented, in concordance with the measurements for dppD-3xCi ptc and dppD-1xCi ptc obtained in [10]. These results (sumarised in Table 1 row 3)) are balancing a twofold consequence of the reduction in the number of enhancers from 3 to 1. At one hand, the reduction implies the vanishing of any possible kind of cooperativity between TFs. In the case of total cooperativity between TFs these relations are symmetric for activators and repressors, so their disappearance reduces the signalling, that is less transcription in activated cells and more transcription in repressed cells, but it does not modify the balance between net activated or repressed cellular ranges. In the case of asymmetric cooperativity, that is, partial cooperativity either only between activators or repressors, the cooperative specie is loosing that advantage. This would provoke a global reduction of activation, in the case of activators cooperativity, and repression in the case of repressor cooperativity. This forces the NAC range variation: when the cooperativity between activators is removed the NAC range is reduced and it increases when the repressor cooperativity is abolished. The interpretation of these assertions on the contrary Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 7 / 25 allow us to affirm that one of the roles of asymmetric cooperativity between activators/repressors is to increase/decrease the NAC range with respect to the NAC range in the non cooperative case. A summary of these results can be found in Table 1 row 2). On the other hand, regardless of the cooperativity, the reduction in the number of enhancers also implies that signalling has to be weakened. Both considerations explains the theoretical transcriptional effects of the binding sites reduction. Fig 2 provides particular examples of the different behaviours in the NAC range under the same experiment in presence of the cooperativities previously mentioned. Under our modelling, the only transcriptional logic compatible with the experiments is the one which occurs in presence of partial cooperativity between repressors. Nevertheless we can find more concordances in this direction by using other experimental evidences. Experiment 2: Differential affinity effects The second experiment we focus on is the comparison of net transcriptional rates of reporter genes containing either three high-affinity sites version of the dpp enhancer (dppD-3xCi ptc ) or three low-affinity dpp sites (dppD-Ci WT ). It was observed that higher Ci affinity provides a reduction in the net activated cellular region, that is a relevant intermediate region where net activated cells for (dppD-Ci WT ) are net repressed for (dppD-3xCi ptc ). It was also observed that increased affinity provides stronger activation in the region close to the A/P border as well as a stronger repression in regions far from the same border (see Fig 2D in [10]). If we accept that the cooperativity between TFs is working in the same way in the binding to both, dpp and ptc, versions of the binding sites, our analysis suggests that the effect of the reduction in affinity could again depend on the type of cooperativity occurring between TFs: Table 1. Transcriptional logics in the presence of opposing A/R gradients. Biochemical characteristics a) Non/Total cooperativity b) Act. partial cooperativity c) Rep. partial cooperativity 1) #TFs affinity #Sig, NAC )#CER #NAC )#CER "NAC )"CER 2) #TFs-TFs coop. #Sig, NAC )#CER #Act, #NAC )#CER #Rep, "NAC )"CER 3) #no. enhancers no. = 3 to no. = 1 #Sig, NAC )#CER #NAC )#CER "NAC )"CER 4) #A-RNAP coop. #Act, #NAC )#CER #Act, #NAC )#CER #Act, #NAC )#CER 5) #R-RNAP coop. #Rep, "NAC )"CER #Rep, "NAC )"CER #Rep, "NAC )"CER Key: "increase, #decrease, no change, )produces, NAC net activated cellular range, CER cellular expression range. This is a simplified comparison of the transcriptional response to different biochemical characteristics between two genes controlled by opposing activator/repressor gradients. The column headings are the kind of TF cooperativity analysed: non/total cooperativity (TFs can cooperate with any other TF), partial cooperativity only between activators or partial cooperativity only between repressors. The biochemical characteristics are: row 1): affinity of TFs for their binding sites, rows 2): cooperativity between TFS, row 3): number of enhancers, row 4) and 5): cooperativity between TFs and RNAP. The variation of affinity considered in row 1) is proportionally equivalent for both activators and repressors. The table shows whether the cellular expression ranges increase or decrease ("CER, #CER) and how it works. Decreases in cooperativity between TFs and RNAP, rows 4) and 5), again produce globally higher/lower expression rates which cause the increase/decrease in the net activated cellular range ("NAC, #NAC resp.) and CER. The response to differences in the other analysed biochemical characteristics varies depending on the kind of cooperation between TFs. If activators and repressor do not cooperate or cooperate globally the net activated cellular region remains unaltered ( NAC) and the signal is weakened (#Sig) provoking the decrease of activation in net activated cells but also repression in net repressed cells. On the other hand, if partial cooperation occurs between TFs the same biochemical characteristics can produce either increase or decrease of the net activated cellular region ("#NAC) and in consequence broader or narrower expression ranges ("#CER) depending on the type of cooperation. https://doi.org/10.1371/journal.pone.0209349.t001 Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 8 / 25 Fig 2. Transcriptional responses to the experiment 1. First column: schematic of the experiment 1: comparison of the expression ranges of reporter genes with 3 high-affinity sites (dppD-3xCi ptc ) and a single high-affinity Ci site (dppD-1xCi ptc ). A) corresponds to the non/total cooperativity case where, if cooperativity holds, all the TFs cooperate between them, C) to the activators cooperativity case, only activators cooperate, and finally E) to the repressors cooperativity case where only repressors cooperate. Second column shows the different transcriptional responses that can be theoretically described depending on the case of cooperativity considered. The schemes and plots employ the same keys explained in Fig 1. https://doi.org/10.1371/journal.pone.0209349.g002 Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 9 / 25 in the presence of total cooperativity, while ½BAjARjR� ¼ ½B�cðjA1Þþ AcðjR1Þþ R½A� KA � �jA½R� KR � �jR ð4Þ if partial cooperativity for TFs occurs. Here, (�) + denotes the positive part function ((x) + =xif x>0 and (x) + = 0 if x�0). This is needed because the cooperativity will not take place unless two or more cooperative TFs are present in the configuration. In the rest of this paper, we will designate the cases when the TFs cooperate between them totally and partially as {{A,R} c } and ffAgcA;fRgcRgrespectively. Note that this notation covers the case of non cooperativity since it would correspond to the case {{A,R} 1 } or equivalently {{A} 1 , {R} 1 }. The binding sites are ordered spatially and, in general, there is not an unique spatial distribution for a configuration with j A activators, j R repressors and n−j A −j R free sites. For instance, if we consider j A =j R = 1 there are six possible spatial distributions with the same elements (ARO,RAO,AOR,ROA,OAR,ORA where Odenotes the empty space). In our description, spatial localisation of bound particles is not considered. In fact, for a specific configuration with j A activators, j R repressors and j 0 free sites n! j0!jA!jR!different spatial configurations are plausible, where k! denotes the factorial of k. Regarding the promoter’s RNA polymerase binding process, the TFs work together trying to promote or repress the binding process [8] by a mechanism known as recruitment [12,13]. Thus, we consider that the activators interact with RNAP with ‘adhesive’ interaction [22] that gives rise to a modification of the RNA polymerase binding affinity: KRP=ajAwhere ais a cooperativity constant greater than 1. In contrast, the effect of j R repressors is modelled in terms of a ‘repulsive’ interaction that modifies the binding affinity KRP=rjRwith an anti-cooperativity factor r<1 (repressor interaction). We will refer to these parameters as TF transcriptional activation/repression intensity. By using the previous guidelines we will now describe the concentrations of all possible configurations as was done in [15,16]: Step 1: Construction of the sample space. All the possible ways of obtaining an equilibrium concentration with j A ,j R and j P activators, repressors and RNA polymerases is given by the states ZðnÞðjA;jR;jP¼1; CÞ ¼ CðCÞn! j0!jA!jR!½B�½RNAP� KRP a½A� KA � �jAr½R� KR � �jR ;ð5Þ ZðnÞðjA;jR;jP¼0; CÞ ¼ CðCÞn! j0!jA!jR!½B�½A� KA � �jA½R� KR � �jR where j P = 1 means there is a bound RNA polymerase and j P = 0 there is none, j 0 =n−j A −j R �0, and the variable Cdescribes the relation of cooperativity between the TFs. Specifically, by using (3) and (4), the cooperativity function Ctakes the values CðC¼ fA;RgcÞ ¼ cðjAþjR1Þþð6Þ and CðC¼ ffAgcA;fRgcRgÞ ¼ cðjA1Þþ AcðjR1Þþ R:ð7Þ This allows us to describe the entire sample space, i.e. the space of all the possible Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 16 / 25 configurations, by O¼ fðjA;jR;jPÞ;jA;jR�0;jAþjR�n;jP¼0;1g: Step 2: Definition of the probability. Once we have described all the possible configurations in terms of the concentrations of activator, repressor and RNA polymerase, we easily obtain the probability of finding the promoter in a particular configuration of j P RNA polymerase and j A ,j R TFs related by a cooperativity relation Cas PðnÞðjA;jR;jP;CÞ ¼ ZðnÞðjA;jR;jP;CÞ X fj0 A;j0 R;j0 Pg2O ZðnÞðj0 A;j0 R;j0 P;CÞ;ð8Þ for all (j A ,j R ,j P )2O. Step 3: Definition of the BEWARE operator. In this last step, the BEWARE operator is obtained in terms of the probabilities P (n) . Following the work of Shea et al [16] the synthesis of a certain protein depends on the total probability of finding RNA polymerase in the promoter, specifically, the synthesis is proportional to the marginal distribution of the case j P = 1 [10, 22, 23]. This justifies the definition of the BEWARE operator as BEWAREð½A�;½R�;½RNAP�;CÞ ¼ CBX jAþjR�n jA;jR�0 PðnÞðjA;jR;jP¼1; CÞ where in definition (8) expression (5) is assumed and C B is a proportionality constant that could depend on other factors not considered in this work. Splitting the denominator in two sums, when RNA polymerase is bound or not bound to the configuration, this expression can be rewritten in terms of the regulation factor function, F reg : BEWAREð½A�;½R�;½RNAP�;CÞ ¼ CB 1þPj0 Aþj0 R�n j0 A;j0 R�0ZðnÞðj0 A;j0 R;j0 P¼0; CÞ Pj0 Aþj0 R�n j0 A;j0 R�0ZðnÞðj0 A;j0 R;j0 P¼1; CÞ ¼CB 1þKRP ½RNAP�Fregð½A�;½R�;CÞ : ð9Þ Doing some basic algebra, this regulation factor can be reduced to facilitate the understanding of the general process (see Section A in S1 File). This has been done by using a classic strategy employed for obtaining the General Binding Equation more than a century ago [32]. This, have not yet been applied, in this context, to the authors knowledge. In fact, we can prove that the regulation factor can be equivalently written as Fregð½A�;½R�;CÞ ¼ SðnÞða½A�K1 A;r½R�K1 R;CÞ SðnÞð½A�K1 A;½R�K1 R;CÞ;ð10Þ where the explicit expression of SðnÞðx;y;CÞdepends on the kind of cooperativity presumed, Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 17 / 25 that is SðnÞðx;y;ffA;Rg1gÞ ¼ ð1þxþyÞn;ð11Þ SðnÞðx;y;ffA;RgcgÞ ¼ 11 cþ1 cð1þcx þcyÞn;ð12Þ SðnÞðx;y;ffAgcA;fRgcRgÞ ¼ ð1þcAxþcRyÞn cAcRþ11 cR � �ð1þcAxÞn cA þ11 cA � �ð1þcRyÞn cRþ11 cA � � 11 cR � �;ð13Þ for the non cooperative, total and partial cooperative cases respectively. Note that the mathematical complexity in these expressions is mainly related to the assumed cooperativity. Versions of transcription logic in the presence of opposing gradients In this section we are going to describe what the transcriptional reaction of genes, controlled by the same opposing TFs, would be when there are biochemical differences between them. As we explained in the Section Results the consequences of such differences will depend on the type of cooperativity occurring between the TFs. The analysis of the case of single gradients can be found in Section K in S1 File. Transcriptional logic in the case of opposing gradients and non/total cooperativity between TFs. We observe that expression (12) coincides with (11) when c= 1 which allows us to use the same mathematical expression for both cases, non cooperativity and total cooperativity. So, in both cases the transcription rates are given by (9) with Fregð½A�;½R�;fA;RgcÞ ¼ 11 cþ1 c1þac ½A� KAþrc ½R� KR � �n 11 cþ1 c1þc½A� KAþc½R� KR � �nð14Þ being the regulation factor, where c= 1 if there is no cooperativity between TFs and c>1 total cooperativity occurs. In fact, we prove that the transcription logic will be basically the same in both cases. Determination of net activation/repression concentrations and cell activated ranges. Thanks to the BEWARE operator we can theoretically describe which concentrations of activators and repressors will cause higher or lower gene expressions than the basal level, that is, we can describe in great detail the effect of the balance of both signals. Note that the basal state, determined by the absence of TFs, that is [A] = [R] = 0, corresponds to F reg = 1 in expression (9). Thus, the regulation factor describes an effective increase (for F reg >1) or decrease (for F reg <1) of the number of RNAP molecules bound to the promoter, with respect to the basal level, as was stablished in [22]. We can see that in the case of (14) the threshold between activation/repression concentrations (that is F reg = 1) is determined by the linear relation ½R� ¼ a1 1r KR KA½A�ð15Þ dividing the plane [A]−[R] into two parts that we can denominate activation region if ½R�< a1 1r KR KA½A�and repression region if on the contrary ½R�>a1 1r KR KA½A�. See Panel A of S1 Fig where Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 18 / 25 examples of these thresholds are depicted for different values of the parameters. Obviously, the threshold (15) is a linear relation between concentrations of activators and repressors. The steepness of this straight line is: independent of the values cand n, depends on the TFs affinities through the ratio K R /K A , it increases with respect to aand decreases with respect to r (since r<1). This justifies the behaviour of the net activated cellular ranges in the case of non/ total cooperativity stated in Table 1 as we will now explain. We can take this information and by considering appropriate gradients of activators and repressors we can define the tissular regions of net activated and net repressed cells, which are, regions of cells expressing more or less than the basal expression level (see Fig 1 for a detailed explanation). For the sake of clarity, and taking into account that our main goal is to understand how these mechanisms could modify the expression of the Hh target genes, we are going to assume that transcription factors act in the same way as Cubitus works in the Drosophila system. Hh secreted from the posterior into the anterior compartment of the wing imaginal disc results in opposing gradients of activator and repressor Ci. The A/P boundary is located at around 60% of the dorso-ventral (D/V) axis. The influence of Hh gradient can be appreciated in the middle of the anterior compartment, specifically the cells located in the region between the 30% and 60% on the D/V axis, approximately. In order to model these TFs concentration distributions we are going to assume that they do not change over time and both are monotonic along the tissue, ½A� ¼ ½A�ðxÞstrictly decreasing and ½R� ¼ ½R�ðxÞstrictly increasing ð16Þ in terms of x, the distance from the A/P boundary. See example in inset in Panel B of Fig 1 where the [CiA]/[CiR] decreases/increases from the A/P border. Let us justify the monotone character assumed previously on the TFs gradients in Drosophila development. It is well known that CiA concentrations depend on the Hh gradient secreted from the posterior compartment. Hence, we will consider that this concentration decrease with the distance from the A/P border, that is, [A](x) is decreasing. On the other hand, the levels of Ci transcripts are high at the positions of A compartment distal to the A/P boundary whereas are low at the positions proximal to the A/P boundary. Assuming the proportionality between transcription levels and protein concentrations this would imply that the total amount of Ci hðxÞ ¼ ½R�þ½A�;ð17Þ is a non decreasing function in terms of the distance from the A/P boundary, x. Hence, the concentration of CiR must increase with the distance from the A/P border. A particular case of this situation is the conservation of the total amount of TFs proposed in [10] corresponding to ½A� ¼ hex=ffiffiffiD p;½R� ¼ h ½A�;ð18Þ being hthe TFs total concentration and Dis the steepness of the gradient. For the sake of simplicity, we are going to employ hypothesis (18) in all the Figures of this work, although the analytical results will consider the more general case (17). Under (18) the concentrations will be restricted to a straight line in the [A]−[R] plane (see Panel A in S1 Fig and Panels A, C and E in S2 and S3 Figs). Insets in Figs 1–3and S1 and S2 Figs, show the distributions [A], [R]. The intersection points between the straight line (18) and the thresholds, ([A] th , [R] th ) (represented by black circles in Panel A in S1 Fig and Panels A, C and E in S2 and S3 Figs) will determine a boundary between genetically activated and repressed cells. That is, repressed cells will be those containing concentrations ([A], [R]), verifying (18) and [A]<[A] th . For activated cells this would be [A]>[A] th . In consequence, they would express transcription rates lower/higher than the basal. Due to the monotonic nature of the TFs distributions, (16), activated cells are Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 19 / 25 closer to the A/P boundary and the limit of the percentage of the wing imaginal disc occupied by activated cells will be determined by the distance x th given by ½A�th ¼ ½A�ðxthÞ: This limit is represented by blue circles in Panels B, C and D in S1 Fig and Panels B, D and F in S2 and S3 Figs. In any case, it is verified that the lower [A] th the larger is the NAC. That is, ½A�th1<½A�th2!xth1>xth2ð19Þ by the monotone character of the activator gradient as can be seen in S3 Fig. Transcriptional consequences of differential biochemical characteristics. Now, by using previous considerations, we want to justify the transcription logic in presence of total cooperativity or in absence of any cooperativity, results collected in Table 1, column a) and represented in S1 Fig. In this case it is quite easy to see the behaviour of the net activated cellular range. Eq (15), which determines the threshold between activation/repression concentrations, does not depend on the number of enhancers nor the cooperativity cand depends on the TFs affinities in terms of the ratio K R /K A . Thus, these regions do not change ( NAC) for genes gene1 and gene2 such that: 1. The affinity for the binding sites of gene1 is smaller than the affinity for the binding sites of gene2 in a proportional manner. In terms of the dissociation constants, this would be expressed as Kg2 R¼dKg1 Rand Kg2 A¼dKg1 Abeing 0 <δ<1 and this occurs because we are considering proportional change of affinity for activator and repressors, Kg1 R=Kg1 A¼Kg2 R=Kg2 A. 2. The TFs cooperate in less intense manner for gene1 than for gene2, that is c g1 <c g2 . 3. gene1 has less binding sites than gene2, that is n g1 <n g2 . Here, the superscripts g1 and g2 stand for the parameters of the genes gene1 and gene2, respectively. However, in the case of differential affinities, where the proportionality is not verified, the net activated cellular range would change. For instance, if Kg2 R=Kg2 A>Kg1 R=Kg1 Athen the net activated cellular range for gene1 would be narrower than for gene2. The rest of the assertions in Table 1, column a) requiere some simple monotonicity properties that have been checked with Lemmas C and D in S1 File. The biochemical differences, numbered 1), 2) and 3), have been proven to verify that Fg2 reg >Fg1 regð>1Þ;if ½A�and ½R�belong to the activation region and Fg2 reg <Fg1 regð<1Þ;if ½A�and ½R�belong to the repression region: Note that the BEWARE operator is monotonic increasing, with respect to the regulation factor F reg which allow us to extrapolate these estimates to expression rates. These three results have been interpreted as a signal weakening. Biochemical differences 1), 2) and 3) can not modify the character of net activation/repression but are able to make signalling less efficient. That is, they do not change the NAC range but they cause less activation in the activated region and less repression in the repressed region (#Sig). In consequence, we can say that in situations 1), 2) and 3) the expression rates will decrease in the net activated cellular range and increase in the repressed cells, which will attenuate the cellular expression range. See S1 Fig where these variations have been depicted. In contrast, assertion 3) in Lemma D in S1 File implies that Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 20 / 25 • If Ais a weaker transcriptional activator in gene1 than in gene2, exhibiting lower cooperativity between the activators and the RNA polymerase (a g1 <a g2 ), then expression of gene1 will be smaller than in gene2 because Fg1 regð½A�;½R�;fA;RgcÞ<Fg2 regð½A�;½R�;fA;RgcÞ: Obviously, this involves globally lower transcription rates (#Act) and more restricted net activated cellular ranges (#NAC) for gene1 than for gene2. • If Ris a weaker transcriptional repressor in gene1 than in gene2, exhibiting lower anti-cooperativity between the repressors and the RNA polymerase (r g2 <r g1 <1), then expression rate of gene1 will be higher than gene2 expression because Fg1 regð½A�;½R�;fA;RgcÞ>Fg2 regð½A�;½R�;fA;RgcÞ: Obviously, this involves globally higher transcription rates (#Rep) and wider net activated cellular ranges ("NAC) for gene1 than for gene2. Transcription logic in the case of partial cooperativity between TFs. In the case of partial cooperativity, the expression rates are given by (9) where the regulation factor is defined by: Fregð½A�;½R�;ffAgcA;fRgcRgÞ ¼ 1þacA½A� KAþrcR½R� KR � �n cAcRþ11 cR � � 1þcAa½A� KA � �n cAþ11 cA � � 1þcRr½R� KR � �n cRþ11 cA � � 11 cR � � 1þcA½A� KAþcR½R� KR � �n cAcRþ11 cR � � 1þcA½A� KA � �n cAþ11 cA � � 1þcR½R� KR � �n cRþ11 cA � � 11 cR � � : Compared to (14), in this complex representation of the expression rates the activation threshold is not as clear. So we need to do a little more delicate mathematical analysis. In fact, if we impose the threshold equation Fregð½A�;½R�;ffAgcA;fRgcRgÞ ¼ 1;ð20Þ it can be shown that, if n= 3 and c A ,c R �1, this threshold is determined by an unique increasing function f, verifying Fregð½A�;KRfð½A�=KAÞ;ffAgcA;fRgcRgÞ ¼ 1:ð21Þ It separates the concentrations [A]-[R] into those which provoke net activation, when [R]/K R <f([A]/K A ), and those which provoke net repression, when [R]/K R >f([A]/K A )) (see Section E in S1 File for definition and analysis of the function f). Note that the threshold of a BEWARE operator with partial cooperativity is not, in general, a straight line although it shows a linear asymptotic behaviour for large concentrations (see Panels A, C and E in S2 and S3 Figs). Note also that the same analysis implies that the function fis independent of the TF affinities: K A and K R . We can understand the effects of partial cooperativity clearly in certain limit regimes such as: • Cooperativity only between repressors, that is c A = 1 and c R >1. The expression rates have been proved to be monotonic decreasing with respect to repressors cooperativity, that is, the more cooperativity the less expression because cooperativity increases repression effectivity Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 21 / 25 (see Lemma F 2) in S1 File). Then, obviously, less cooperativity implies less repression provoking wider NAC and CER. • Cooperativity only between activators, that is c A >1 and c R = 1. The counterpart results in this other case shows that expression rates are increasing with cooperativity between activators since it increases activation effectivity (see Lemma F 1) in S1 File). Then a reduction in cooperativity between activators reduces NAC and CER. These result have been summarised in Table 1 cells 2b) and 2c) and represented in Panel D in S2 and S3 Figs. When both TFs cooperate simultaneously we have compared the thresholds given by functions fwith the linear relation (15), obtaining a criteria comparison in terms of the variables that follow � a2¼3cAcR a1 r1f2þð1þaÞðr1ÞcA2ar þða1Þðrþ1ÞcRg;ð22Þ � a3¼cAcR a1 ðr1Þ2fð1þaþa2ÞcA 2ðr1Þ23cAðr1Þða2r1Þ  ð1þrþr2Þða1Þ2cR 2þ3ða1Þðar21ÞcRg:ð23Þ Depending on the positive or negative sign of these values it can be proven that the change from partial cooperativity to non cooperativity causes: • If � a2>0,� a3>0: we see an decrement in the net activated cellular range with respect to the non cooperative case. • If � a2<0,� a3<0: we see a increment in the net activated cellular range with respect to the non cooperative case. • In all other cases: An increase or decrease of the net activated cellular range can occur depending on the total amount of TFs considered (h), the binding affinities, K A ,K R as well as the activation/repression intensities. A detailed explanation can be found in Section E in S1 File. Another interesting aspect is how these thresholds depend on the affinities of the TFs in presence of partial cooperativity. Let us consider again gene2 with high affinity binding sites and gene1 with proportionally low affinity binding sites, that is, Kg2 R¼dKg1 Rand Kg2 A¼dKg1 Abeing 0<d<1:ð24Þ Let us now consider the function f. It defines the threshold of a BEWARE operator with partial cooperativity, determined by the cooperativity constants a,r,c A and c R and expression (21). This function and the corresponding affinities determine the activator concentration thresholds for gene1 and gene2,½A�g1 th and ½A�g2 th , using these expressions ½A�g1 th þKg1 Rfð½A�g1 th =Kg1 AÞ ¼ hðxÞand ½A�g2 th þKg2 Rfð½A�g2 th =Kg2 AÞ ¼ hðxÞ;ð25Þ where (17) has been considered. Then, under hypothesis (24), the order of ½A�g2 th and ½A�g1 th can be determined. It can be proven that when only activators cooperate between them these limit values verify ½A�g2 th <½A�g1 th Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 22 / 25 implying that gene2 has a wider expression range than gene1 because of (19). On the other hand, when only repressors cooperate the inverse relation its true ½A�g2 th >½A�g1 th implying that gene2 has a narrower expression range than gene1 (see Lemma J in S1 File for details). In Panels D and F in Fig 3 both these situations have been depicted when partial cooperativity is either for activators or for repressors. More detailed representations about the thresholds and these range determination process can be found in Panels A and B in S2 and S3 Figs. Previous results can be generalised to the case of both species cooperating simultaneously and monotonicity has be related with conditions of concavity and convexity of the threshold function f. This condition can be understood in terms of the Greater-Than-Additive and LessThan-Additive effects described in transcriptional activation [23]. The concavity of fis the same as observing Greater-Than-Additive effects at the basal level. That is, when fis concave, any convex combinations of two different pairs of concentrations, [A], [R], that give basal expression, will produce net activated cells. On the other hand, when fis convex Less-ThanAdditive effects can be observed. This result, and our interpretation, are in contrast with the analysis in [11] for the case of cooperative repressors. The assertions in Table 1, cells 4b) and 5c) can be verified easily from statement 3) Lemma F, in S1 File. Supporting information S1 File. Supplementary material. This file contains all the mathematical details of the theoretical analysis. Furthermore, it contains the discussion of the results that can be deduced for the single gradient case. (PDF) S1 Fig. Transcriptional logic in the presence of opposing gradients and non/total cooperativity. (PDF) S2 Fig. Transcriptional logic in the presence of opposing gradients and partial cooperativity between activators. (PDF) S3 Fig. Transcriptional logic in the presence of opposing gradients and partial cooperativity between repressors. (PDF) S4 Fig. Transcription logic in the case of single gradients. (PDF) Acknowledgments O. S. would like to thank professor J. Garcı ´a-Ojalvo for pointing out reference [22] and R. Bielawski for language editing and proofreading. The authors also would like to thank the reviewers for their contributions. Author Contributions Conceptualization: Manuel Cambo ´n, O ´scar Sa ´nchez. Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 23 / 25 Data curation: Manuel Cambo ´n, O ´scar Sa ´nchez. Formal analysis: Manuel Cambo ´n, O ´scar Sa ´nchez. Investigation: Manuel Cambo ´n, O ´scar Sa ´nchez. Methodology: Manuel Cambo ´n, O ´scar Sa ´nchez. Resources: Manuel Cambo ´n, O ´scar Sa ´nchez. Validation: Manuel Cambo ´n, O ´scar Sa ´nchez. Visualization: Manuel Cambo ´n, O ´scar Sa ´nchez. Writing – original draft: Manuel Cambo ´n, O ´scar Sa ´nchez. Writing – review & editing: Manuel Cambo ´n, O ´scar Sa ´nchez. References 1. Tabata T, Takei Y. Morphogenes, their identification and regulation. Development. 2004; 131: 703–712. https://doi.org/10.1242/dev.01043 PMID: 14757636 2. Torroja C, Gorfinkiel N, Guerrero I. Mechanisms of Hedgehog Gradient Formation and Interpretation. Dev Neurobiol. 2005; 64: 334–356. https://doi.org/10.1002/neu.20168 3. Hallikas O, Palin K, Sinjushina N, Rautiainen R, Partanen J, Ukkonen E, et al. Genome-wide Prediction of Mammalian Enhancers Based on Analysis of Transcription-Factor Binding Affinity, Cell. 2006; 124 (1): 47–59. https://doi.org/10.1016/j.cell.2005.10.042 PMID: 16413481 4. Lorberbaum DS, Ramos AI, Peterson KA, Carpenter BS, Parker DS, De S, et al. An ancient yet flexible cis-regulatory architecture allows localized Hedgehog tuning by patched/Ptch1. eLife. 2016; 5: e13550. https://doi.org/10.7554/eLife.13550 PMID: 27146892 5. Ramos AI, Barolo S. Low-affinity transcription factor binding sites shape morphogen responses and enhancer evolution. Philos T R Soc Lon B. 2013; 368(1632): 20130017. 6. Gurdziel K, Lorberbaum DS, Udager AM, Song JY, Richards N, Parker DS, et al. Identification and Validation of Novel Hedgehog-Responsive Enhancers Predicted by Computational Analysis of Ci/Gli Binding Site Density. PLoS ONE. 2015; 10(12): e0145225. https://doi.org/10.1371/journal.pone.0145225 PMID: 26710299 7. Alexandre C, Jacinto A, Ingham PW. Transcriptional activation of hedgehog target genes in Drosophila is mediated directly by the cubitus interruptus protein, a member of the GLI family of zinc finger DNAbinding proteins. Gene Dev. 1996; 10:2003–2013. https://doi.org/10.1101/gad.10.16.2003 PMID: 8769644 8. Mu¨ller B, Basler K. The repressor and activator forms of Cubitus interruptus control Hedgehog target genes trough common generic Gli-binding sites. Development. 2000; 127: 2999–3007. PMID: 10862738 9. Aguilar-Hidalgo D, Becerra-Alonso D, Garcia-Morales D, Casares F. Towards a study of gene regulatory constraints to morphological evolution of the Drosophila ocellar region. Dev Genes Evol. 2016; 226: 221–233. https://doi.org/10.1007/s00427-016-0541-8 PMID: 27038024 10. Parker DS, White MA, Ramos AI, Cohen BA, Barolo S. The cis-Regulatory Logic of Hedgehog Gradient Responses: Key Roles for Gli Binding Affinity, Competition, and Cooperativity. Sci Signal. 2011; 4: 1– 16. https://doi.org/10.1126/scisignal.2002077 11. White MA, Parker DS, Barolo S, Cohen BA. A model of spatially restricted transcription in opposing gradients of activators and repressors. Mol Syst Biol. 2012; 8: 614. https://doi.org/10.1038/msb.2012.48 12. Ptashne M, Gann A. Transcrition activation by recruitment. Nature. 1997; 386(10): 569–577. https://doi. org/10.1038/386569a0 PMID: 9121580 13. Ptashne M. Regulation of transcription: from lambda to eukaryotes. Trends Biochem Sci. 2005; 30(6): 275–279. https://doi.org/10.1016/j.tibs.2005.04.003 PMID: 15950866 14. Piepenburg O, Vorbru¨ggen G, Ja ¨ckle H. Drosophila segment borders result from unilateral repression of hedgehog activity by wingless signaling. Mol Cell. 2000; 6:203–209. https://doi.org/10.1016/S10972765(05)00011-0 PMID: 10949042 15. Ackers GK, Johnson AD, Shea MA.Quantitative model for gene regulation by λphage repressor. P Natl Acad Sci USA.1982; 79: 1129–1133. https://doi.org/10.1073/pnas.79.4.1129 Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 24 / 25 16. Shea M, Ackers GK. The O R Control System of Bacteriophage Lambda. A Physical-Chemical Model for Gene Regulation. J Mol Biol. 1985; 181: 211–230. https://doi.org/10.1016/0022-2836(85)90086-5 PMID: 3157005 17. Aguilar-Hidalgo D, Dominguez-Cejudo MA, Amore G, Brockmann A, Lemos MC, Co ´rdoba A, et al. A Hh-driven gene network controls specification, pattern and size of the Drosophila simple eyes. Development. 2013; 140: 82–92. https://doi.org/10.1242/dev.082172 PMID: 23154412 18. Aguilar-Hidalgo D, Lemos MC, Co ´rdoba A. Core regulatory network motif underlies the ocellar complex patterning in Drosophila melanogaster. Physica D. 2015; 295: 91–102. https://doi.org/10.1016/j.physd. 2014.12.006 19. Gilman A, Arkin AP. Genetic “Code”: Representation and Dynamical Models of Genetic Components and Networks. Annu Rev Genom Hum G. 2002; 3: 341–369. https://doi.org/10.1146/annurev.genom.3. 030502.111004 20. Ay A, Arnosti DN. Mathematical modeling of gene expression: a guide for the perplexed biologist. Crit Rev Biochem Mol. 2011; 46(2): 137–151. https://doi.org/10.3109/10409238.2011.556597 21. Frank TD, Cavadas MAS, Nguyen LK, Cheong A. Non-linear Dynamics in Transcriptional Regulation: Biological Logic Gates. In: Carballido-Landeira J, Escribano B, editors. Nonlinear Dynamics in Biological Systems, SEMA SIMAI Springer Series. 2016; 7: 43–62. 22. Bintu L, Buchler NE, Garcı ´a HG, Gerland U, Hwa T, Kondev J, et al. Transcriptional regulation by the numbers: models. Curr Opin Genet Dev. 2005; 15: 116–124. https://doi.org/10.1016/j.gde.2005.02.007 PMID: 15797194 23. Frank TD, Carmody AM, Kholodenko BN. Versatility of Cooperative Transcriptional Activation: A Thermodynamical Modeling Analysis for Greater-Than-Additive and Less-Than-Aditive Effects. PLoS ONE. 2012; 7(4): e34439. https://doi.org/10.1371/journal.pone.0034439 PMID: 22506020 24. Zhou X, Su Z. tCal: transcriptional probability calculator using thermodynamic model. Bioinformatics. 2008; 24(22): 2639–2640. https://doi.org/10.1093/bioinformatics/btn494 PMID: 18796477 25. Rogers KW, Schier AF. Morphogen Gradients: From Generation to Interpretation, Annu Rev Cell Dev Bi. 2011; 27: 377–407. 26. Junker JP, Peterson KA, Nishi Y, Mao J, McMahon AP, van Oudenaarden A. A Predictive Model of Bifunctional Transcription Factor Signaling during Embrionic Tissue Patterning. Dev Cell. 2014; 31: 448–460. https://doi.org/10.1016/j.devcel.2014.10.017 PMID: 25458012 27. Balaskas N, Ribeiro A, Panovska J, Dessaud E, Sasai N, Page KM, et al.Gene regulatory logic for reading the Sonic Hedgehog signaling gradient in the vertebrate neural tube. Cell. 2012; 148: 273–284 https://doi.org/10.1016/j.cell.2011.10.047 PMID: 22265416 28. Kicheva A, Bollenbach T, Ribeiro A, Valle HP, Lovell-Badge R, Episkopou V, et al. Coordination of progenitor specification and growth in mouse and chick spinal cord. Science. 2014; 345: 1254927. https:// doi.org/10.1126/science.1254927 PMID: 25258086 29. Saha K, Schaffer DV. Signal dynamics in Sonic hedgehog tissue patterning. Development. 2006; 133: 889–900. https://doi.org/10.1242/dev.02254 PMID: 16452094 30. Verbeni M, Sa ´nchez O, Mollica E, Siegl–Cachedenier I, Carlenton A, Guerrero I, et al. Morphogenetic action through flux–limited spreading. Phys Life Rev. 2013; 10: 457–475. https://doi.org/10.1016/j.plrev. 2013.06.004 PMID: 23831049 31. Lai K, Robertson MJ, Schaffer DV. The Sonic Hedgehog Signaling System as a Bistable Genetic Switch. Biophys J. 2004; 86: 2748–2757. https://doi.org/10.1016/S0006-3495(04)74328-3 PMID: 15111393 32. Bisswanger H. Enzyme Kinetics: Principles and methods. 2nd ed. Weinheim: WILEY-VCH; 2008. Analysis of transcriptional logic governing Hh target genes PLOS ONE | https://doi.org/10.1371/journal.pone.0209349 January 7, 2019 25 / 25