Quantitative models of cis-regulatory activity have the potential to improve our mechanistic understanding of transcriptional regulation. However, the few models available today have been based on simplistic assumptions about the sequences being modeled, or heuristic approximations of the underlying regulatory mechanisms. We have developed a thermodynamics-based model to predict gene expression driven by any DNA sequence, as a function of transcription factor concentrations and their DNA-binding specificities. It uses statistical thermodynamics theory to model not only protein-DNA interaction, but also the effect of DNA-bound activators and repressors on gene expression. In addition, the model incorporates mechanistic features such as synergistic effect of multiple activators, short range repression, and cooperativity in transcription factor-DNA binding, allowing us to systematically evaluate the significance of these features in the context of available expression data. Using this model on segmentation-related enhancers in Drosophila, we find that transcriptional synergy due to simultaneous action of multiple activators helps explain the data beyond what can be explained by cooperative DNA-binding alone. We find clear support for the phenomenon of short-range repression, where repressors do not directly interact with the basal transcriptional machinery. We also find that the binding sites contributing to an enhancer' s function may not be conserved during evolution, and a noticeable fraction of these undergo lineage-specific changes. Our implementation of the model, called GEMSTAT, is the first publicly available program for simultaneously modeling the regulatory activities of a given set of sequences.
We began by examining whether our models agree with existing data on transcriptional gene regulation during Drosophila embryonic development (anterior–posterior axis specification). This involved training our model on 37–44 experimentally characterized CRMs and 6–8 transcription factors. The overall quality of fit as well as predictive ability of our models was remarkably high. Next, we applied different model variants to investigate mechanistic questions. We found that the transcriptional synergy arising from simultaneous contact of activators with the BTM contributes significantly to the accurate specification of expression patterns, and this contribution extends beyond the contribution from mutual interactions (DNA-binding cooperativity) between activators. Shifting attention to repressors, we then found that competition between repressors and activators for binding sites is an insufficient mechanism of repression [29]. We found evidence in favor of a short range repression mechanism for two of the TFs, consolidating experimental evidence that exists for this mechanism. However, our results also raised the possibility that long-range mechanisms (such as direct interaction with the BTM) may also contribute to the repressors' function. We also studied the importance of cooperative DNA-binding (of both activators and repressors) in this system. Our results provide clear evidence of cooperative effects of some TFs but give mixed signals with respect to other TFs. We also used our model to examine a contentious evolutionary issue. Several studies [30]–[32] have reported that TF binding sites undergo rapid turnover (loss and/or gain) during evolution. However, due to the difficulty of establishing true functionality of binding sites in practice (e. g. , binding to a TF does not necessary lead to regulatory function [33]), it is not clear whether such turnover is largely limited to non-functional sites. We investigated this issue using our model in conjunction with evolutionary sequence comparison, and found that lineage–specific losses affect functional sites to a noticeable extent. As mentioned above, a few thermodynamics-based models have been proposed in the past, which we now review briefly. The approach of Reinitz and colleagues exploits physicochemical principles, and includes important mechanistic aspects such as short range repression through quenching [13], [34]. However, the Reinitz model does not consider all possible molecular configurations, a fundamental tenet of the statistical thermodynamic treatment. Also, cooperative DNA-binding by TFs is not included in the model. Segal et al. [6] presented a model based on enumeration of all configurations of bound and unbound TFs. This model uses statistical thermodynamics to model TF-DNA interactions and to compute relative probabilities of configurations, but models the mapping from these configurations to their transcriptional output in a heuristic manner. Also, the Segal model ignores important mechanistic issues such as transcriptional synergy (discussed above) and short range repression. Furthermore, the formulation of transcriptional output in this model makes the computational task intractable. (The authors adopted sampling methods to deal with this issue, thereby sacrificing exactness of the model computation.) The models developed by other researchers make various simplifying assumptions, e. g. , binding of a single activator is strong enough to activate transcription [14], and their implementations are often limited in their generality, e. g. , only sequences with a small number of binding sites are considered [12], or all sites are assumed to have identical binding affinities [14]. See Table 1 for a summary of the strengths and weaknesses of the models discussed above. We have not undertaken a rigorous comparison of our approach versus the above-mentioned approaches, for three reasons. First, none of the previous models have a publicly available implementation that we could use in our setting. Bauer & Bailey' s implementation [11] of the Reinitz model comes closest, but cannot be applied to more than one CRM at a time. Second, while Segal et al. [6] make their data set (and their predictions for this set) available, their method uses a much larger number of free parameters (the position weight matrices of the TFs were estimated from data), precluding a fair comparison. Third, and most importantly, our main goal in this study was to search for mechanistic insights that are revealed by the data, rather than engineering a model with the best possible fit to the data. For the same reason, we have not attempted here to position our work in comparison to machine learning-based models of gene expression [35], [36]. To facilitate future studies by other researchers, we make the source code of our implementation freely available online.
We begin with an overview of the statistical thermodynamic theory of gene expression, following Buchler et al. [7]. This theory has two components, one dealing with the occupancy of TFs at DNA sequences, and the other with the interactions of occupied TFs with the BTM. We first describe the model of TF occupancy. Consider a CRM with n binding sites (e. g. , n = 2 in Figure 1A). A molecular configuration, denoted by σ, specifies which sites are bound and which are free. Thus there are 2n possible configurations. The statistical weight of configuration σ, denoted by W (σ), and which we shall endeavor to compute, gives us the relative probability, P (σ), of the configuration when the system is in equilibrium. In other words, we have P (σ) = W (σ) /Z, where Z is a normalization constant, defined as ∑σ W (σ), and known as the partition function. Calculation of P (σ) would allow us to answer questions like: “What is the relative probability of site S being in the bound state? ” This may be computed by summing P (σ) over all σ in which S is bound, and is also called the fractional occupancy of the site S. The statistical weight W (σ) depends on the number and affinities of the occupied binding sites in the configuration σ, and on interactions between bound TF molecules. We will present details of W (σ) when discussing specific models below. We next describe, at a high level, how the above molecular configurations (σ) affect gene expression. We assume that the gene expression level (on a scale of 0 to 1) is equal to the fractional occupancy of the promoter by the BTM. Each of the configurations σ considered above (specifying bound or unbound TFs) may now correspond to two states, depending on whether BTM is bound or not. The statistical weight of these two states will be given by W (σ) Q (σ) and W (σ) respectively, where W (σ) is the contribution from TF–DNA interactions as explained above, and Q (σ) is the contribution from TF–BTM interactions, present only if the BTM is bound. Q (σ) depends on the bound TFs in the configuration σ, and may be construed as the transcriptional output from the configuration. We now calculate the relative probability of bound BTM as ZON = Σσ W (σ) Q (σ), and that of unbound BTM as ZOFF = Σσ W (σ), to obtain the gene expression level as follows (note that “ON” and “OFF” represent the state of BTM occupancy, which is separate from the occupancy states of binding sites in the CRM sequence, as indicated by σ): (1) Here, we present details of how the W (σ) and Q (σ) terms are specified by the first of our two models. Under this model, DNA-bound transcription factors interact favorably (activators) or unfavorably (repressors) with the BTM, thus affecting the probability of the BTM being bound at the promoter. We call this the Direct Interaction (“DirectInt”) model. For a configuration σ, the statistical weight W (σ) has terms reflecting binding of TFs to their binding sites, and those reflecting interactions between TFs. Let q (S) denote the contribution of a single occupied site S to W (σ). This depends on the concentration of the TF and the strength of the site, and is given by: (2) where (See Text S1 for how Equation (2) is derived.) Note that two unknown constants, one related to TF-DNA binding (K (Smax) ), and the other (ν) a constant of proportionality for TF concentration, appear together as a product, which can be treated as a single free parameter to be estimated from data. The above equation makes the implicit assumption that the binding energy of each position of a site is additive. This assumption has been questioned in several studies [37], [38], but is necessary in our case because there is not enough TF-DNA interaction data to construct accurate models incorporating higher-order interactions. Furthermore, the additivity assumption seems to be a reasonably good approximation for the TFs in the segmentation system [6], [13]. The statistical weight of a configuration σ, in the absence of cooperative binding, is then given by, where σi is an indicator variable (values 0 or 1) for Si being occupied by its TF in the configuration (Figure 1A) [7]. If two bound TFs interact (protein–protein interaction), they make an additional contribution to the statistical weight of the configuration. We denote this contribution by ω (d), where d is the distance between their binding sites (Figure 1B). The dependence of this cooperativity term on the distance is discussed in Text S1. The statistical weight of a configuration, accounting for cooperative binding, is the product of contributions of all occupied sites and all TF-TF interactions implied by that configuration [7]: (3) where ωij (dij) denotes the statistical weight contribution due to interaction between the TFs bound to sites Si and Sj, and dij is the distance between these sites. We assume that cooperative binding is possible only if the bound sites are adjacent in the configuration, i. e. , there is no other bound site in between. We also assume that it is predetermined whether any given pair of transcription factors exhibit cooperative binding or not. The model allows interactions between adjacent binding sites that may be either homotypic (of the same TF) or heterotypic (of different TFs). Next, we describe how we model Q (σ), the statistical weight contribution from TF-BTM interactions. We assume that each TF is either an activator or repressor. A bound activator A interacts with the bound BTM with statistical weight αA>1, while a repressor R interacts with weight αR<1 (Figure 1C). Q (σ) is the product of the α terms corresponding to each bound TF in the configuration. This corresponds to the intuition that a bound activator makes the configuration more energetically favorable (thus, a greater weight) while a bound repressor makes it less favorable. We also assume that each bound TF interacts independently with the BTM, with energy contributions that add up, which is reflected in the statistical weights being multiplicative. Computation of Equation (1) involves summation of an exponential number of configurations. We developed an efficient algorithm based on dynamic programming to carry out the computation (see below and Text S1). We note that Gertz et al. made the same model assumptions [12], but, unlike their method, we provide a general and efficient implementation that is applicable to arbitrary sequences. The DirectInt model presented here largely follows Buchler et al. [7], with the treatment of sequence-specific DNA binding (Equation 2) being borrowed from Berg & von Hippel [9]. However, the approach of Buchler et al. [7], designed for prokaryotic systems, assumed repressors to work by competition with the polymerase, and does not extend to distally located binding sites. In the DirectInt model above, repressor action is independent of the location of binding sites for repressors or activators. However, experimental work has shown that certain repressors act on activators only if they are bound within a “short range”, e. g. , less than 150 bp, of the activator binding site [17]. Such short range repression, also called “quenching” [17], may work by repressors inhibiting DNA-binding of activators [39], possibly by modifying chromatin accessibility. We model this mechanism by assuming that a bound repressor does not directly interact with the BTM, instead, it creates a new possible configuration, one where DNA in its “neighborhood” (defined by a range parameter dR) is inaccessible to binding by any other TF, for example by localized chromatin modification (Figure 1D). A configuration where the neighboring chromatin is inaccessible (Figure 1D, bottom) competes with the configurations where the chromatin is accessible to activators, thus effectively reducing the occupancy of activators. We call this model the short-range repression, or SRR, model. Note that there are more configurations under this model than in the DirectInt model. In any configuration, an activator site may exist in one of two states (bound or unbound) as in DirectInt. In contrast, each repressor site may now exist in one of three states: unbound, “bound-only”, and “bound-effective” (the bound-only state has the repressor bound but not interacting with either the BTM or the neighboring DNA, while in the bound-effective state the bound repressor makes the neighboring DNA inaccessible). Not all possible configurations are allowed, however. We assume that within a certain range of a bound-effective repressor, an activator site is not allowed to be bound (thus implementing the idea of short-range repression). For a legitimate configuration σ, W (σ) in the SRR model is given by Equation (3), multiplied by a repressor-specific constant βR for each bound-effective site of the repressor R (Figure 1D, bottom). The parameter βR may be interpreted as the equilibrium constant of the reaction that changes the chromatin state from accessible to inaccessible. The value of βR controls the strength of the repressor. When it is close to 0, there is no repression effect; when it approaches +∞, the repressor completely shuts down all activator sites in the neighborhood. Thus, in this alternative to the DirectInt model, repression is modeled by augmenting the calculation of W (σ), instead of direct interaction terms (αR) for the repressor in Q (σ). Q (σ) is now a product of the direct interaction terms (αA) for activators alone. We show that even with this new model, it is possible to perform efficient computation of Equation (1) using dynamic programming (see Text S1). We consider the following question: how are the effects of multiple bound activators combined? In both models described above (DirectInt, SRR), their individual statistical weights (αA) were multiplied, in calculating the overall contribution of activator-BTM interactions. This is the “multiplicative effect” model of combined action by multiple activators. It reflects a scenario where the bound activators interact with different parts of the BTM (or different steps of transcription initiation), and the energy terms are added. Veitia [26] shows that this multiplicative effect model results in “transcriptional synergy”, where the activating effect of two binding sites is greater than the sum of their individual effects, even in the absence of cooperative DNA-binding. We next consider an alternative scenario where in any given configuration, at most one activator molecule may interact with the BTM. This is plausible if for example the bound activators must interact with the same subunit of the BTM. In this case, the TF-BTM interaction term is written as Q (σ) = Σαi, where the sum is over bound activators in the configuration. This is the called the “additive effect” model (Figure 1E). In this case, there will be no synergistic activation due to TF-BTM interaction, though some level of synergy may still arise from cooperative DNA-binding by activators. In Text S1, we compare the two mechanisms that may lead to transcriptional synergy: multiplicative effect model, and additive effect model in combination with cooperative DNA binding. The basic insight is that synergistic effect will disappear at high activator concentration under the cooperative binding model (activator binding has already been saturated under this condition, thus cooperative interactions will not be further helpful), but not under the multiplicative model. This difference in model behavior suggests that it is theoretically possible to distinguish two models from the data. To investigate the mechanism of synergistic activation, we implement both “multiplicative effect” and “additive effect” models as special cases of a more general model for combined activator action: a user-defined parameter NMA (positive integer) sets the limit on the maximum number of bound activators that can simultaneously interact with the BTM. We call this the “limited contact” model of activator action (see Text S1 for details). The cases NMA = 1 and NMA = ∞ correspond to the additive and multiplicative effect models respectively. This general model can be combined with cooperative binding of TF molecules, thus allowing us to study the relative contribution of multiplicative activation and cooperative binding. As discussed earlier, the computation of Equation (1) involves summation of an exponential number of configurations. In this section, we describe an efficient algorithm for computing the DirectInt model with multiplicative effect of activation. (The algorithms for other models are based on similar dynamic programming techniques and are presented in Text S1.) Let ZOFF (i) denote the total statistical weight of all configurations of sites up to the site i, with site i being occupied. We obtain the following recurrence, by summing over the position of the occupied site j nearest to site i: (4) where q (i) is the statistical weight of the site i, as defined in Equation (2), ω (i, j) is the interaction between the occupied sites i and j, and Φ (i) is the set of sites to the left of i that do not overlap with i. This recurrence equation is similar to that in [40], [41]. The constant term, +1, corresponds to the case where no site to the left of i is occupied. Under this model, Q (σ) is the product of the transcriptional effects (α terms, as described before) of all occupied TF molecules in the configuration σ. Let f (i) be the factor bound at the site i, and αf (i) be the transcriptional effect of f (i), then we have the following recurrence for ZON: (5) To calculate the values required for Equation (1), we simply take the sum over all possible values of i: and. The time complexity of the algorithm is O (n2), where n is the number of sites in the sequence. However, if cooperative interaction between adjacent sites is not modeled, or the interaction only occurs within a constant range, the time complexity is linear in n. We started with the Drosophila segmentation data set from Segal et al. [6]. This set includes 44 bona fide CRMs with their A/P expression profiles, eight TFs (Bcd, Cad, TorRE, Hb, Gt, Kni, Kr, and Tll) with their concentration profiles and PWM motifs. Each expression profile (or concentration profile) consists of 100 real numbers between 0 and 1 representing the relative expression level of the CRM (or relative concentration of a TF) in positions along the A/P axis, divided into 100 bins from anterior to posterior. One problem with this data set is that not all relevant TFs in the terminal regions (e. g. , Slp1) are included or known [42]. Also, the TorRE (Torso Response Element) motif included in this data set is assumed to correspond to a (yet unknown) TF that has activating role in the terminal regions of the embryo. Recent evidence suggests that on the contrary TorRE may correspond to the Capicua transcription factor, which is a repressor expressed in the trunk region of the embryo and post-transcriptionally degraded at the termini in response to Torso signaling [43]. This casts doubts over the inclusion of TorRE and in general the terminal regions of the expression profiles as part of the data set, especially for evaluating models that distinguish between activator and repressor mechanisms. We thus limited the CRM expression profiles to their portions lying between 20% and 80% egg length. The number of CRMs came down to 37, after excluding those without patterned expression in this spatial range. This final data set included six motifs (Tll and TorRE were excluded), of which five (Cad, Gt, Hb, Kr, Kni) were taken from Noyes et al. [44] and one (Bcd) was obtained from FlyREG [45]. Binding sites were annotated as those with log likelihood ratio (LLR) scores greater than 0. 4 times the LLR score of the optimal site [46]. This threshold is weak enough to include a large number of putative sites for each TF, while keeping the running time low. Parameter training was performed using the Nelder-Mead simplex method and the quasi-Newton method, and restarts were used to deal with potential local optima. Optimization of the correlation coefficient between predicted and known expression values was alternated with optimization of the sum of squared errors (See Text S1 for details). Note that model training is performed separately for each model (DirectInt or SRR, with or without cooperative interaction, etc.). Thus, even though two models may share certain parameters, e. g. , (K (Smax) ν), their values may be different under the two models after training. The running time of the program scales linearly with the number of TFs, and the total length of sequences (for all models except the “limited contact model”, see Text S1). In our dataset, with 6–8 TFs and about 40 CRMs of average length 1450 bp, the parameter training phase took about 3–4 hours of running time on a desktop computer with 2. 2GHz CPU and 2GB memory.
GEMSTAT offers the following choices between various model features: The above choices are accompanied by parameters that may be set manually, and some of which may be left as free parameters to be trained from the data. All model parameters are described in Table S1. The program takes as input the sequence and expression profiles of a set of CRMs, and the PWMs and concentration profiles of a set of TFs. Expression profiles and concentration profiles are specified as vectors of a fixed dimension, allowing it to be easily used to model any regulatory system. (In our application, vector components correspond to positions along the A/P axis of the embryo, but in other applications these could be distinct anatomical domains or temporal points.) The source code is available at http: //veda. cs. uiuc. edu/Seq2Expr/. The data set consists of 37 experimentally characterized CRMs driving patterned expression along the anterior-posterior axis of the blastoderm stage Drosophila embryo (see Methods). We used several different approaches to objectively evaluate a model and compare models. Our first metric of model performance is the correlation between the model predictions and the observations. For each CRM, we calculated the Pearson correlation coefficient (CC) between the predicted and the observed expression profiles (over 60 bins), and computed the average CC over all CRMs. We also recorded the number of CRMs with CC>0. 65. Additionally, we estimated for each CRM the significance of improvement (in CC) due to one model versus another, and combined these estimates into a p-value of improvement over the entire data set, as described in Text S1. We also calculated the average CC under 10-fold cross validation (denoted by “CVCC”), as a test of predictive ability, and for fair performance comparison between models with different numbers of parameters. For any given choice of model, an identical model that uses randomly permuted PWMs was evaluated as negative control. Any observation about model comparison based on correlation coefficients was also confirmed by visual inspection of the predicted expression patterns on all 37 CRMs. We note that there is no consensus yet on the most reasonable way to evaluate predictions of expression models for data sets such as that used here. We chose the correlation coefficient because of its ability to capture the salient pattern along the A/P axis, and we based all of our claims on this measure to keep our analysis objective and unbiased. We began by exploring the effect of cooperative DNA-binding by molecules of the same TF, i. e. , homotypic interactions. (Modeling heterotypic interactions would involve many more free parameters and was not pursued in this study). Segal et al. [6] also studied this effect, but since their model lacks mechanistic details of activation, the effect of cooperative binding may not be distinguishable from simultaneous interaction of TFs with the BTM (the “multiplicative effect”, also called “transcriptional synergy” [8]). As a baseline, we evaluated the DirectInt model that excludes any cooperative binding terms, but allows for transcriptional synergy. The average correlation coefficient (CC) of this model (of 13 free parameters) is 0. 547, with accurate predicted readout (CC>0. 65) on 16 of the 37 CRMs (Table 2). In contrast, 25 independent negative controls yielded a mean average CC of 0. 211 (standard deviation of this mean across the 25 trials was 0. 075). The cross validation correlation coefficient (CVCC) supports the high predictive ability of the model (average CVCC of 0. 4, compared to 0. 02±0. 083 from negative controls). We then included self-cooperativity of each TF separately (only one additional parameter at a time), and computed the average CC and CVCC as before. Each of the TFs showed an improved CVCC over the baseline of no cooperativity across almost every replicate of the cross validation exercise (Table 2, Table S2), while Bcd and Kni showed the most pronounced effects of cooperativity in terms of average CC. When both Bcd and Kni were included as cooperative factors, the average CC improved further over the model with each factor alone. The improvement in going from no cooperativity (average CC = 0. 547) to self-cooperativity for Bcd and Kni (average CC = 0. 587) was highly significant (p-value 1. 3E-6). Examination of the expression predictions on individual CRMs identified 12 CRMs where the cooperativity model was better and two where it was worse. (Two cases are shown in Figure 2 (A, B), and the complete list is in Figure S2.) Our results are consistent with Segal et al. [6], who found self-cooperativity to improve prediction. Moreover, we find this to be the case even in the presence of transcriptional synergy, which if not accounted for could have confounded the effects of cooperative DNA-binding by activators. As a visual aid for interpreting the quantitative evaluations reported above, we present in Figure 3 all of the expression predictions from the above model (with Bcd and Kni self-cooperativity), alongside their respective known expression patterns. A detailed summary of the model' s performance is given in Table 2, along with results from an appropriate negative control. This model was also fit to the entire data set of Segal et al. (44 CRMs, inclusive of terminal bins) and found to have slightly (but not significantly) higher average CC than the published predictions of the Segal model [6], although our model uses fewer free parameters (see Figure S3 for details.) GEMSTAT implements two alternative approaches to combining the effects of multiple activator sites, using the parameter NMA described in Table S1: the additive effect (NMA = 1) and the multiplicative effect (NMA = ∞), as well as approaches that are in between these two extremes. The “multiplicative effect” model allows any number of activator molecules to simultaneously interact with the BTM, which as discussed in Methods, leads to transcriptional synergy, a source of synergistic activation that is distinct from cooperative DNA-binding [8], [26]. We used the two extreme values of NMA to test whether this phenomenon leads to improved agreement with the data, while keeping other aspects of the model fixed (Table 3). The baseline model here was one with NMA = 1 (no synergy) and with no self-cooperative DNA-binding. The average CC from this model (0. 516) improved significantly (to 0. 547; p-value 3. 7E-4) when we introduced synergy due to the multiplicative effect of multiple activators (NMA = ∞). (This change does not involve any additional free parameters.) This was further confirmed by a greatly improved CVCC (0. 295 to 0. 40, see Table 3 and Table S3), as well as by examination of predictions for individual CRMs (Figure 2 (C, D), and for the complete results see Figure S4): the multiplicative effect model showed clear improvements on 6 CRMs and was worse on 3 CRMs. These observations suggest that simultaneous interaction of multiple activators with the BTM is a plausible source of synergistic activation. Cooperative binding was kept out of the model in the above test. We next introduced cooperative binding (only for the two activators) into the model, and examined the contribution of the multiplicative effect. We found that the model with both sources of synergistic activation shows better quality of fits compared to the model with cooperative binding alone, in terms of average CC (from 0. 558 to 0. 581, p-value 7. 3E-11, see Table 3) as well as CVCC (0. 292 to 0. 396). We also confirmed this improvement by examination of individual CRMs: the model using multiplicative effect along with cooperative binding led to better fits for 8 CRMs compared to the model with cooperative binding alone (Figure S5) and was worse in no case. This result suggests that synergistic activation due to multiplicative effect of activators may be over and beyond that due to cooperative binding [26]. In all of the above tests, we had used a “Direct Interaction” model of repressor function, where a bound repressor is assumed to interact directly with the BTM, destabilizing the configuration, and thus curbing the roles of activator sites in the entire CRM. GEMSTAT also allows us to deploy a more “localized” form of repressor action, in the form of the short range repression (SRR) model, where a bound repressor makes the neighboring chromatin (within some range dR) inaccessible. Prior experimental work [20] suggests that the four repressors in our data set – Kr, Hb, Kni, and Gt – act over short distances (∼100–150 bp [14]), and in two of these cases (Kr and Kni) repression depends on the histone deacetylase dCtBP, which suggests a possible mechanistic basis for the short range action [20]. In our tests, we sought to examine if the SRR model implemented in GEMSTAT is realistic enough to capture the repressors' contributions to expression patterns. Starting with a baseline where every repressor was modeled by “Direct Interaction”, we introduced the SRR model for one repressor at a time (with dR = 250bp), and compared the resulting model with the baseline. Although none of the four resulting models (Kr-, Hb-, Kni-, Gt-SRR) showed clear improvement over the baseline, we found strong evidence that for Kr and Hb, the SRR model implemented by GEMSTAT was able to capture the repressive effects of the factors almost to the same extent as the Direct Interaction model, as described next. We first extended our evaluation metric, the average CC, in the following way: we considered the best K CRMs for a model (in terms of CC), and plotted the average CC over these K CRMs, for all values of K (1 … 37). We found the Kr-SRR model to be highly similar (in terms of average CC) to the baseline model throughout the range (Figure 4A). Additionally, for each model and each value of K, we plotted the average CC of the same model under a Kr “knock down” condition, i. e. , where the Kr concentration was set to 0 across the A/P axis. Such a “knock down” plot allows us to visualize the contribution of a TF (Kr in this case) to the model. We found Kr to contribute significantly to both models, although the contribution to the SRR model was not as strong as to the DirectInt model. This may reflect certain limitations of the SRR implementation in GEMSTAT, but the results strongly suggest that the short range nature of Kr action [20] is largely captured by our model. We also examined the performance differences between the models on individual CRMs. We found seven CRMs where the SRR model was as effective as or better than the Direct Interaction model in predicting readout, with a significant contribution from Kr (Figure S7). In five other cases, the Direct Interaction model yielded superior fits (plots not shown). Similar evidence for the effectiveness of the Hb-SRR model is shown in Figures 4B and Figure S8. However, the Gt-SRR model does not seem to elicit significant contribution from Gt, even though this repressor is found to be effective within the DirectInt model (Figure S6A). A similar lack of evidence was encountered for the Kni-SRR model (Figure S6B). Thus, we find that for Kr and Hb, quenching of activator sites within a distance is sufficient to capture the repressive effect of the TFs, supporting the hypothesis that these two TFs act mainly as short-range repressors, confirming what has been reported in earlier studies, which admittedly relied on a small number of CRMs and synthetic enhancers. On the other hand, we did not find strong evidence of short-range repression for Gt and Kni, and even for Kr and Hb the SRR model' s performance was only as good as and not better than the DirectInt model. This is somewhat unexpected; it may be in part due to limitations of our SRR model, but may also be hinting that these TFs use long-range repression mechanisms as well (see Discussion). Repression by competitive binding, as proposed in the literature [7], [29], involves the binding of repressors to sites overlapping activator sites, thereby suppressing their occupancy by activators. This mechanism may be thought of as a special case of the SRR model in GEMSTAT, with the repression range parameter (dR) set to ∼10 bp. At such a small value of dR, a repressor can only make its immediate neighborhood inaccessible, equivalent to inactivation of overlapping activator sites. Having observed above that the Kr and Hb repressors are effectively modeled in the SRR mode, we compared the Kr-SRR and Hb-SRR models at dR = 250 to their respective versions at dR = 10. As shown in Figure 4 (C, D), in both cases the competitive binding model (dR = 10) was significantly worse than the SRR model, both in terms of average CC and in terms of the repressor' s contribution. Finally, we sought to use the GEMSTAT program to probe an important question regarding the function and evolution of transcription factor binding sites. A number of recent studies have reported the “turnover” (evolutionary gain and loss) of binding sites, based on sequence comparison [30], [31], [47] or from ChIP-based experiments [48]. However, it is possible that such lineage-specific loss and gain are largely limited to non-functional sites, i. e. , the false positive matches to PWMs, or sites that are bound by TFs but do not regulate expression [33], [49], [50]. Here, we explored this possibility by asking if sites that change in lineage-specific ways are functional in contributing to the expression patterns. We note that lineage-specific losses may in part be artifacts of alignment errors (i. e. , sites were completely conserved but not deemed so, due to misalignment). However, in practice, the true gain/loss of sites may be hard to distinguish from alignment errors, so we will call both cases as lineage-specific changes here. We predicted sites by demanding that any predicted site be conserved (in the sense of being above threshold) in all species analyzed, and examined how the quality of fit varies as this evolutionary filter was made more stringent by including more species. We found that more conservative evolutionary filters lead to greatly reduced average CC (Figure 5, red). This shows that a noticeable part of the CRMs' functionality is carried by sites (in D. melanogaster) that are not found to be conserved across all phyla. Those sites could, broadly speaking (a) be deeply conserved in the examined phylogeny, but with some lineage specific losses or (b) have arisen specifically in D. melanogaster or a recent ancestor. Next, we modified the evolutionary filter to demand deep (but not necessarily complete) conservation across the phylogeny (see Text S1) and found that above-mentioned loss in quality of fits disappears (Figure 5, blue, compare number of species = 2 vs. 6). Since the new evolutionary filter discards sites of type (b) mentioned above, we inferred that a noticeable part of the CRMs' functionality is carried by sites that are largely conserved but also undergo lineage-specific losses.
Three different mechanisms have been previously hypothesized for repressor action: (i) competition with activators for access to binding sites, (ii) direct interference with BTM recruitment and assembly, and (iii) local interference (“quenching”) with the function of nearby activator sites. The last hypothesis seems to be most likely in the context of the regulatory system we analyze, as suggested by the following observations: first, repressor and activator sites are often found to be close to each other [19]; second, CRMs of the same gene often work independently, i. e. , a repressor site within one module does not stop the function of another module for the same gene [54]; third, some repressors are found to depend on a co-repressor, CtBP, a histone deacetylase that presumably increases the association of nucleosomes to DNA, making it less accessible [20]. However, direct evidence of this so-called short-range repression (SRR) phenomenon is limited to a few CRMs and synthetic enhancers [20]. We implemented models that could investigate all three mechanisms with respect to their agreement with data on a moderate number of CRMs. Note that even though the short range mechanism has been implemented (in other forms) previously [13], [14], it has never been tested within a framework that also implements alternative mechanisms. We report the first direct data-based comparison between alternative hypotheses regarding repression. Our results clearly exclude the hypothesis of competitive binding being the main mechanism of repression, and are consistent with the SRR hypothesis for two of the four repressors studied (Kr and Hb). It is somewhat unexpected that the SRR model does not explain the data as well as the DirectInt model for Gt and Kni. We note that while Gt is believed to be a short range repressor, Nibu et al. [20] leave open the possibility of this protein having long range mechanisms of action as well, in light of the fact that it does not require dCtBP to mediate repression. Similarly, Kr has been found to have long range mechanisms as well [21], [55]. It is also likely that to some extent the inability of the SRR model to match (for Gt and Kni) or exceed (for Kr and Hb) the effectiveness of the DirectInt model arises from shortcomings of our model and evaluation procedure. Our model assumes that once a repressor molecule is bound, it will make its entire neighborhood inaccessible, defined by a range parameter. We would intuitively expect that the repression effect is stronger for closer chromatin regions, and this is not modeled due to our lack of understanding of the exact mechanism by which repressors may change the chromatin structure. Similarly, we do not know exactly how the effects of two repressor molecules are combined in the regions that may be affected by both, and this part is treated in a simplistic manner under our SRR model. The dataset may also limit our ability to study detailed mechanisms: the resolution of expression patterns is low and the dataset lacks informative negative controls (all sequences are wild type CRMs). Finally, our tests are likely to have been weakened by the fact that models are compared on individual CRMs and not entire control regions. It is generally assumed that the short range mode of repressor action is necessary for the functional modularity of CRMs. For example, Kr is a key input to the eve stripe 2 enhancer, but it can adversely affect the expression readout of the adjacent eve stripe 3 enhancer; this interference is avoided presumably because of its short range of action [56], [57]. Thus, the effect of SRR is already manifested in the compactness of CRMs, and if it were possible to compare SRR with the direct interaction model on entire gene control regions, we would likely observe a clear advantage to the former. Despite these limitations, the SRR model along with a detailed activation model allows to ask questions that cannot be addressed with simple non-mechanistic models of CRM function. Another important issue we explored is how multiple activator sites contribute to expression. It is likely that this multiplicity is important for the synergistic activation, where the total effect of multiple sites is larger than the sum of their individual effects. That such synergy is real and important has been shown through in vitro experiments on the effect of the number of sites [51], as well as in vivo experiments on the typically sharp boundaries of gap gene expression domains [2]. Mechanistically, synergy may result either from cooperative DNA binding of multiple activator molecules or from simultaneous interaction of multiple activators with the BTM (Text S1) [26]. Our model implements both mechanisms, and is thus able to examine the effect of each mechanism on readout, both in the absence and in presence of the other mechanism. We found that both mechanisms are involved in setting the precise expression profile; the effect of transcriptional synergy is evident, and complementary to that of cooperative binding. We have not explored in this study some important details on how synergistic interactions with BTM may occur, and these may worth further investigation. For example, we did not make any distinction between different activators. It is plausible that two different activators may interact with BTM simultaneously, contacting different subunits [52], while the two molecules of the same TF may act in an additive fashion, contacting the same subunit. Other possibilities remain to be explored with regard to cooperative DNA-binding as well. One possibility stems from our assumption that only two adjacent bound molecules may interact with each other. Although this assumption has been commonly made in other studies dealing with cooperativity [6], it is based partly on computational considerations and partly on our lack of understanding of the mechanistic details of interactions among TF molecules. On the topic of mechanistic limitations of our models, we note also that in equating gene expression to the fractional occupancy by the BTM, we are ignoring the internal dynamics of transcription initiation and elongation [58]–[60]. We found that for a number of CRMs, the model (mis-) predicts expression outside the CRM' s primary expression domain (s). For instance, the CRM “kni_ (-5) ” drives anterior expression only, but the model additionally predicts modest expression in the central and posterior regions of the embryo (Figure 3). We noted that kni_ (-5) has many binding sites for Cad, which is an activator present in the posterior half of the embryo. Presumably, the model fails to find strong evidence of appropriate repressive influence, and predicts kni_ (-5) to drive expression in the posterior regions, mediated by the putative Cad sites. A similar observation was made with respect to the CRM “eve_stripe5”, which drives expression in the posterior half (in a stripe between bins 60 and 70, see Figure 3). This CRM harbors several high quality putative sites for Bcd, which is an anterior activator, and this is presumably the reason why the model predicts modest anterior expression as well. That such incongruous predicted expression is often seen under multiple models suggests that the errors may not be due to the specifics of the model that we have been varying. Rather, it is possible that we are missing some additional repression mechanism, e. g. , from chromatin modifications, from unknown repressor sites, or mischaracterization of binding affinity. A relevant fact worth noting here is that there is some ambiguity about the appropriate binding profile to use for the important repressor Gt. In the current study, we used the profile estimated from in vitro Bacterial-one-hybrid (B1H) experiments [44], which happens to be quite different from the profile estimated from verified Gt binding sites in DNA footprinting experiments [45]. However, because relatively few sites were verified, the footprinting-based Gt profile is too un-specific to be used for prediction of new sites. We observed that the total number of Gt sites in all CRMs is considerably smaller than most other factors. This may have led to underestimation of the repressive influence of Gt, and a consequent lack of repression (as per the model' s predictions) in the region where Gt is expressed. An important area of future improvements to our approach will be the quality and amount of data. The spatial expression profiles used here were obtained from manual parsing of stained (in situ hybridization) images, and are essentially qualitative. This is one of the reasons why our evaluations were based on correlation between expression patterns rather than more absolute measures of prediction accuracy. More accurate quantifications that are under way [61] should lead to improved analysis. Our approach assumes that the expression profiles of TFs and CRMs were synchronized (from the same developmental time), although this is not entirely true: the temporal resolution of the data set is not high enough to ensure such synchronization, and this is another direction where future, higher resolution data sets will be needed. Moreover, since we do not have data characterizing the dynamic state of chromatin (nucleosome distributions and their chemical modifications), we did not explicitly model the changes of chromatin structure that may be induced by TF association. With more high-quality expression data and ideally more epigenetic data as well, it should be possible to extend our models with additional details and to incorporate theoretical models of chromatin structure [62], [63]. The models presented here are intended to be usable in a variety of regulatory systems in different species. It is true however that a regulatory system would need to be very well understood at a qualitative level and characterized by quantitative measurements at multiple levels, before we can apply such models. We would need the following information to train the models: (1) the expression readouts of a set of promoters or CRMs, (2) a reasonably complete set of TFs involved in the regulatory network, (3) quantification of their concentration profiles, and (4) their binding specificities. At this time, such a data set is often not available, making it difficult to evaluate the generalizability of the models. A promising application of the proposed quantitative models lies in the prediction and characterization of novel CRMs. Once a sequence-expression model is trained, it may be applied genome-wide to predict segments that have the potential to direct the expression patterns of neighboring genes. The model may also be used to predict the effect of individual transcription factor perturbations, leading us to individual TF-CRM interactions. This paradigm requires quantitative measurements of TF levels, a requirement that may be mitigated to some extent by using mRNA levels of TF genes, but ideally by direct protein level measurements. Recent developments in proteomics and in high-throughput assays of post-translational modifications offer great hope in providing the necessary TF activity data [64]. The models offer new ways to approach the study of regulatory sequence evolution. Transcription factor binding sites have been reported to undergo frequent loss and gain, but it is not clear what the functional consequences of these changes are. We saw an example of how the functional context provided by the model may be combined with cross-species sequence comparison to provide new insights into binding site turnover. In general, sequence-expression models allow us to predict the changes in expression pattern that result from any evolutionary change at the sequence level. This interpretative power may be harnessed to investigate how regulatory sequences evolve under different schemes of selection, and begin to answer questions such as “With gene expression under purifying selection, how tolerant is a sequence to the gain and loss of binding sites? ” or “How feasible is it to evolve a novel expression pattern using only simple nucleotide level changes, i. e. , substitutions, insertions and deletions? ” [65]. Quantitative models have a natural relevance in the field of synthetic biology. In order to design gene networks with a well-defined input/output characterization, we need the ability to engineer gene promoters or enhancers that direct specific expression patterns (outputs) in response to the specific levels of the regulators (inputs). This ability in turn requires a tool to predict the expression pattern corresponding to any given sequence. Moreover, to search in a very large sequence space, an efficient sequence-to-expression mapping will be crucial. This will be a place where our dynamic programming-based algorithms make a large difference. In the long run, we expect quantitative models to be able to consider for example the entire intergenic region next to a gene (and not only individual CRMs) and predict the gene' s spatial-temporal expression pattern. The GEMSTAT models are an important preliminary step towards this grand goal.
The development of complex multicellular organisms requires genes to be expressed at specific stages and in specific tissues. Regulatory DNA sequences, often called cis-regulatory modules, drive the desired gene expression patterns by integrating information about the environment in the form of the activities of transcription factors. The rules by which regulatory sequences read this type of information, however, are unclear. In this work, we developed quantitative models based on physicochemical principles that directly map regulatory sequences to the expression profiles they generate. We evaluated these models on the segmentation network of the model organism Drosophila melanogaster. Our models incorporate mechanistic features that attempt to capture how activating and repressing transcription factors work in the segmentation system. By evaluating the importance of these features, we were able to gain insights on the quantitative regulatory rules. We found that two different mechanisms may contribute to cooperative gene activation and that repressors often have a short range of influence in DNA sequences. Combining the quantitative modeling with comparative sequence analysis, we also found that even functional sequences may be lost during evolution.
genetics and genomics/gene expression
computational biology/transcriptional regulation
biophysics/theory and simulation
computational biology/macromolecular sequence analysis
biophysics/transcription and translation
Thermodynamics-Based Models of Transcriptional Regulation by Enhancers: The Roles of Synergistic Activation, Cooperative Binding and Short-Range Repression