Next Article in Journal
Mathematical Modeling of Recursive Drug Delivery with Diffusion, Equilibrium, and Convection Coupling
Previous Article in Journal
Generalized Thermoelastic Interaction in a Half-Space under a Nonlocal Thermoelastic Model
Previous Article in Special Issue
Explainable Machine Learning for Longitudinal Multi-Omic Microbiome
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Thermodynamic Modelling of Transcriptional Control: A Sensitivity Analysis

by
Manuel Cambón
and
Óscar Sánchez
*,†
Applied Mathematics Department, University of Granada, E18071 Granada, Spain
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2022, 10(13), 2169; https://doi.org/10.3390/math10132169
Submission received: 29 April 2022 / Revised: 11 June 2022 / Accepted: 18 June 2022 / Published: 22 June 2022
(This article belongs to the Special Issue Models and Methods in Bioinformatics: Theory and Applications)

Abstract

:
Modelling is a tool used to decipher the biochemical mechanisms involved in transcriptional control. Experimental evidence in genetics is usually supported by theoretical models in order to evaluate the effects of all the possible interactions that can occur in these complicated processes. Models derived from the thermodynamic method are critical in this labour because they are able to take into account multiple mechanisms operating simultaneously at the molecular micro-scale and relate them to transcriptional initiation at the tissular macro-scale. This work is devoted to adapting computational techniques to this context in order to theoretically evaluate the role played by several biochemical mechanisms. The interest of this theoretical analysis relies on the fact that it can be contrasted against those biological experiments where the response to perturbations in the transcriptional machinery environment is evaluated in terms of genetically activated/repressed regions. The theoretical reproduction of these experiments leads to a sensitivity analysis whose results are expressed in terms of the elasticity of a threshold function determining those activated/repressed regions. The study of this elasticity function in thermodynamic models already proposed in the literature reveals that certain modelling approaches can alter the balance between the biochemical mechanisms considered, and this can cause false/misleading outcomes. The reevaluation of classical thermodynamic models gives us a more accurate and complete picture of the interactions involved in gene regulation and transcriptional control, which enables more specific predictions. This sensitivity approach provides a definite advantage in the interpretation of a wide range of genetic experimental results.

1. Introduction

Fundamental biological processes, such as morphogenesis and tumor development, are the consequence of complex signalling pathways that provoke cell responses [1,2,3,4]. This work focuses on the modelling of one of the key mechanisms mediating these responses: the transcriptional gene regulation by protein binding to the cis-regulatory region. The RNA polymerases (RNAP) bind to a DNA sequence called a promoter, and its transcriptional activity is directly promoted or repressed by other DNA-binding proteins called transcription factors (TFs) [1,4]. The linkage of both RNAP in the promoter and the TFs in the transcription factors binding sites (TFBSs) or enhancers is a complex process that depends on several factors: different families of binding proteins (characterised in general as activators, repressors and RNAP), the relative concentration of these proteins in each cell of the tissue, natural tendency of the proteins to be bound to the cis-regulatory region (affinity), and the rules that drive their binding process, among others. This process in itself can be compounded, where the binding can involve competition for the same TFBSs, or assistance of other already bound elements (binding cooperativity) [5,6,7,8]. Transcriptional control results in the gene (protein) production rates being controlled by the bound transcription factors which in turn use several biochemical mechanisms [9,10,11]. These rates can also depend on post-transcriptional regulation mechanisms [12,13,14].
This process is quite widespread and takes place in many prototypical systems, for instance the morphogens Hedgehog (Hh) and Sonic Hedgehog (Shh). In particular, we will use as a reference the experimental evidence in Drosophila’s wing imaginal disc obtained in reference [15]. In the imaginal disc, composed of the Anterior and Posterior compartments, 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 [16,17]. The morphogens spreading give rise to a gradient in the receptor cells which causes opposing gradients of activator and repressor forms of the transcription factor Cubitus (Ci). Under opposing activators/repressors gradients, we get a simple picture of gene expression where the tissue can be divided into regions of gene activation and repression according to the abundance of activators and repressors.
However, it is known that the same signal of Hh produces different spatial expression of these 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. In order to understand how the same signal and the same TFs can produce different responses, the authors of [15] performed specific alterations of TFBS of these genes located in the vicinity of their corresponding promoters. By electrophoretic mobility shift assays, it was found that Ci binding sites in the ptc TFBS have considerably higher affinity than dpp sites. Indeed, the same authors constructed transgenic fly lines that allow them to compare the transcriptional activity of reporter genes: (a) containing different variants of these sites modifying their affinity and (b) reducing from 3 to 1 the number of TFGS in order to detect cooperativity effects.
The aim of this work is to develop accurate theoretical tools to achieve a deeper understanding of these experiments. Indeed, the numerous variables mentioned above make it difficult solely by using experimental procedures. For this reason, experimental discussions in literature are frequently supported by theoretical models in order to decipher how these mechanisms interact [15,18]. The modelling of transcriptional processes has been tackled from different mathematical perspectives [19,20,21,22,23], such as Bayesian [24,25,26,27,28] and Boolean [29], among others. In this work, we will focus on models based on the statistical thermodynamic equilibrium approach [30,31,32,33,34]. This methodology, introduced in the pioneering works of Ackers–Shea and coworkers [30,31], is also known as the BEWARE method (Binding Equilibrium Weighted Average Rate Expression) because transcription, and thus expression, is considered to be proportional to the probability of transcript initiation [35]. This modelling considers transcript initiation an average of all the possible micro-state configurations where the system (proteins/binding sites) can be present.
As already mentioned in [36], the thermodynamic models follow one of two different biological control processes: the “recruitment” or “stimulated” transcription models [10,11,37]. In the recruitment approach, all the configurations with a bound RNAP have the same transcriptional efficiency, but the TFs control the RNA polymerase recruitment by TFs-RNAP cooperative/anti-cooperative interactions. That is, the TFs are able to alter the RNA polymerase affinity for the promoter both positively (activators) and negatively (repressors) [32,38,39]. Alternatively, in the stimulated transcription approach, the RNA polymerase binding affinity is assumed to be fixed, but the translational strength of any configuration is modulated in terms of the bound activators/inhibitors [18,40,41]. In this work, we also consider Hill-type models since they can be justified from thermodynamic models [42,43,44,45,46,47]. See [48] for a comprehensive introduction to these kinds of modelling and [19,36,49,50] and references therein for a review of case-study scenarios implemented in different areas. Although recruitment is considered to be the main way of controlling gene transcription by Ci/Gli factors [51,52], Hh/Shh pathways have been modelled in either using the recruitment [15,53,54,55,56,57,58], the stimulated [18,41,59,60] or Hill approaches [46]. One of the main goals of this work is to show that modelling is sensitive to the choice of one of these approaches.
It is important to note that the averaging procedure performed in the thermodynamic methodology reveals one of its main drawbacks: the complexity of the deduced mathematical expressions when the model takes into consideration the wide variety of biochemical mechanisms that can be involved in the original process (competition of multiple ligands for multiple TFBSs [61], affinities, roles played by each ligand, signal strength, and cooperativity [62,63,64], among others). The large number of micro-states produce entangled mathematical expressions, so the only way to extract valid information from the models is to go through a multiple-parameter non-trivial calibration process [19,65].
In this regard, in [53], the authors were able to analyse the variations of the sizes of the activated/repressed regions determined by a recruitment thermodynamic model (involving different biochemical mechanisms) under analogous perturbations to the experiments performed in [15]. The authors establish there that, if cooperativity between some of the TFs exists, this would suppose a competitive advantage for such TFs in the linking process. In that case, a theoretical emulation of the experiments developed in [15] predicted the alteration of such a competitive advantage which explains variations in net activated/repressed tissular regions, in concordance with the results in [15]. However, it was not clear if the results obtained in [53], or even the results in [15], could be different using alternative BEWARE operators (stimulated, Hill).
In order to answer these questions, in the present work, we have developed a brand new tool based on sensitivity analysis [66,67]. This allows a deep theoretical analysis of the role of affinity, cooperativity and enhancers in all the BEWARE operators considered. This tool also allows for extracting relevant biological information in a straightforward way. Indeed, we have used it to detect specific parameter relations in the models that lead to unexpected results. That is, the models predictions are sensitive to the election of one of these versions of the BEWARE operators.
The paper is structured in the following sections: In Section 2, the guidelines for the deduction of the BEWARE operators with respect to the stimulated, recruitment and Hill approaches (see details in Supplementary Materials Sections S1 and S2) will be reviewed. Although some particular cases of these expressions are well known in the basic literature, the general simplified expressions deduced in Supplementary Materials Section S1 (valid for an undetermined number of TFs) are new and necessary for the subsequent development. In certain circumstances, one can unify all the activators/repressors to get a global activator variable and a global repressor variable as has been justified in Supplementary Materials Section S3. In Section 2.1, the activation/repression thresholds that can subsequently be deduced from the BEWARE operators are introduced. The approach employed in this work to justify the existence of these thresholds, based on the analysis of ‘inverse logic’ in Supplementary Materials Section S4, is new and quite convenient in order to analyse all the considered models. By using them, an activated/repressed tissular region in the presence of (global) activator/repressor opposing gradients can be established. The sensitivity analysis of the activation/repression thresholds with respect to the analogous variations to the experiments performed in [15] are proved to be determined by the elasticity function in Section 2.2. A comparative analysis for recruitment, stimulated and Hill thresholds will be performed based on their elasticity estimates presented in Section 2.3. In addition, finally, in Section 3, we explore the origin of the differences obtained in Section 2.3. In Section 4, we summarise all the conclusions we have made from the performed analysis.

2. Materials and Methods

In this work, we have computed brand-new expressions of the BEWARE operators integrating a variety of already proposed models. We have worked in a general framework, where we assumed that the gene expression was controlled by any number of cooperative transcription factors. However, for the sake of simplicity, we will present the results in this manuscript for only two transcription factors, A and R (see Supplementary Materials Section S1 for more general expressions). The goal of the statistical thermodynamic model is to describe the synthesis rate of a protein P in terms of the activator, repressor and RNA polymerase concentrations ( [ A ] , [ R ] , [ R N A P ] ), and the cooperativity interactions between them ( C ). That is,
d [ P ] d t = BEWARE m ( [ A ] , [ R ] , [ R N A P ] ; C ) β [ P ]
where BEWARE m ( ) is the function specifying the dependence on the TFs, and β [ P ] is a degradation contribution [68]. Thus, stationary states for Equation (1) are essentially determined by the BEWARE operator, which is subjected to the biochemical mechanisms involved in all these binding processes (see Figure 1). The model takes into account these inputs:
  • Variable number of TFBSs (enhancers): In general, a regulatory machinery is governed by specific cis-regulatory regions, bound by transcription factors and RNAP. In our model, we will work with n transcription factor enhancers, and one binding site for the RNAP (promoter).
  • Binding affinities: The binding process of both transcription factors and RNA polymerase is defined by the binding energy (affinity). This affinity will be described by dissociation constants K A , K R and K R P , for the Activator, Repressor and RNAP proteins, respectively. Please note that in our description the higher the value of K i , the lower the binding affinity of the protein i for their corresponding cis-regulatory region, with i = A , R , R P .
In (1), the subindex m is used to denote what version of model is used:
  • m = r stands for “Recruitment” model,
  • m = s stands for “Stimulated” model,
  • m = H r or m = H s stands for “Hill” model deduced either from the recruitment or the stimulated model.
These versions conform to different model hypotheses:
  • Recruitment model: the synthesis of P depends on the total probability of finding RNA polymerase in the promoter. This probability is obtained from a combination of all the enhancer micro-states where the promoter is occupied by the RNAP. Each micro-state takes into account the number of activators and repressors, j A and j R , and their TFs-RNAP cooperative/anti-cooperative effect modifies the RNA polymerase binding affinity by the promoter
    K ¯ R P = K R P a j A r j R ,
    with a 1 , 0 r 1 , and K R P - K ¯ R P being the “TF empty-occupied” binding affinity of the RNA polymerase. That is, activators/repressors facilitate/impede RNAP binding, and therefore transcription.
  • Stimulated model: the synthesis of P depends on a weighted combination of probabilities of finding the enhancers occupied. These weights are determined for each configuration in terms of the constants, r b a s , ν m a x ( n ) , and r ˜ 1 being r b a s , r b a s + ν m a x ( n ) , r ˜ n r b a s the basal, maximal and minimal transcription rates corresponding to either empty or filled with n activators/repressors. The superscript in ν m a x ( n ) denotes the dependence of this maximal rate on the number of enhancers n in the same manner that the minimal rate does. This dependence is a point not considered in standard modelling of the stimulated BEWARE operator that we will discuss in Supplementary Materials Section S1.
  • Hill versions of both Recruitment and Stimulated operator have been obtained with any number of transcription factors (see specific results in Supplementary Materials Section S2). In order to do this, we have adopted the “extreme cooperativity approach” [48] according to several cooperativity regimes explained in the next paragraphs.
Furthermore, in Equation (1), the set C arranges the transcription factors in subsets of proteins. That is, depending on how the transcription factors interact (from now on, ‘cooperate’) between them, we will aggregate these in different sets:
  • Total cooperativity: C { [ A ] , [ R ] } c , where both activators and repressors cooperate with a cooperativity constant c 1 .
  • Partial cooperativity: C { { [ A ] } c A , { [ R ] } c R } , where the activators cooperate only between other activators with a cooperativity constant c A 1 , and the repressors cooperate only with other repressors with a cooperativity constant c R 1 .
This means that, in the presence of j A - j R already bound activators-repressors, the affinity of a new cooperating i protein for an empty enhancer would be K i / c j A + j R in the case of total cooperativity and K i / c i j i in the case of partial cooperativity, where i = A , R . That is, as soon as the cooperativity coefficients are greater than 1, the already bound proteins assist in the binding of any other future binding protein that cooperates with the already bound protein.
It is important to note that, by using different versions of modelling and different cooperativities, we get different explicit BEWARE operators:
BEWARE r ( [ A ] , [ R ] , [ R N A P ] ; C ) = C B 1 + K R P [ R N A P ] F r e g ( [ A ] , [ R ] ; C ) ,
BEWARE s ( [ A ] , [ R ] , [ R N A P ] ; C ) = r b a s 1 + K R P [ R N A P ] B a s a l ( [ A ] , [ R ] ; C ) + ν m a x ( n ) 1 + K R P [ R N A P ] P r o m o t e r ( [ A ] , [ R ] ; C ) ,
where F r e g ( ) , B a s a l ( ) and P r o m o t e r ( ) are functions that only depend on the concentrations of transcription factors, and C B , r b a s and ν m a x ( n ) are scale factors. Please note also that Equations (3) and (4) will have a total of four explicit expressions depending on what kind of cooperativity C applies between the TFs: Two for the total cooperative case, two for the partial cooperative case. Another four Hill models can be deduced by the extreme cooperative approach (see Supplementary Materials Sections S1–S3 for their derivation). That is,
Recruitment ( B E W A R E r ( · , · , · ; C ) ) E x t .   C o o p .   A p p r . Hill ( B E W A R E H r ( · , · , · ; C ) ) Stimulated ( B E W A R E s ( · , · , · ; C ) ) E x t .   C o o p .   A p p r . Hill ( B E W A R E H s ( · , · , · ; C ) )
where C refers to total or partial cooperativity.
Binding cooperativity is not the only way the proteins can interact. With the stimulated approach, TFs can take into account the situation where two adjacent activators or repressors function more effectively together than separately. This sort of functional cooperativity has been considered in some models in the literature [40]. In addition, indeed, this will be reflected in the differences found in this work between stimulated and recruitment models as can be seen in Section 3.
Usually, Expressions (3) and (4) are so complex that it is impossible to perform (non-numeric) mathematical analysis, even in the case of just two transcription factors. However, in Supplementary Materials Section S1.5, we have been able to reformulate these expressions, obtaining simpler polynomial expressions that can be treated theoretically. In those cases when this is not possible, the numerical approach can always be be considered.
Please note that, in this section, we have shown the basic notation and concepts needed in order to follow the main work, in the case of two transcription factors. The two-transcription factor case is a relevant biological case which is involved, for instance, in the control of the target genes of the Hedgehog morphogen in Drosophila. However, there are other biological systems, where more than two transcription factors are involved in the control of the same genes, for instance, in control of the Shh target genes in vertebrates executed by the Gli proteins family. In previous literature [18], the functional grouping of several TFs was proposed. The expressions obtained allow us to reconsider the same question. In Supplementary Materials Section S3, there is a discussion on what the conditions of the system should be to be treated, in general, as a system governed by what we call global Activation/Repression variables. This is a possible reduction in the system’s degrees of freedom that resembles the two-transcription factors case. Indeed, these extra hypotheses are not so demanding that in some works they have been already considered [41,59,60]. The results stated in our next section are only valid when these global Activation/Repression variables can be assumed.
Although we can simplify the stimulated and recruitment BEWARE operators originally deduced, these simplifications do not in any way modify them since they have only been rewritten in a more accessible way. This is not the case of the Hill approach, as we will see in our conclusions, because we have noticed that applying the “extreme cooperativity hypothesis”, which is used in this approach, seriously alters the sensitivity of the stimulated and recruitment BEWARE operators.

2.1. Activation/Repression Threshold

Most of the genetic experiments presented in previous literature basically analyse the gene expression that results in the segregation of a protein P. The experiments work with either wild-type expression and/or expressions measured in mutants where some of the biochemical conditions have been modified. Examples are: Shh/Hh target genes in [15,18,46,51,69,70,71], other Drosophila’s target genes [72,73] and prototypical biological systems, such as the λ phage [74,75]. Of all these experimental approaches, of particular interest are those that compare transcription rates against basal levels. Correctly defining the basal transcription level is quite important because it depends on which part of the transcriptional control system you are considering. Since in our case we are focusing on the transcriptional effects due to signalling in the specific module of n enhancers, the basal level is the expression of P observed when these enhancers are disabled from receiving that signalling [15,70]. These basal levels can still be signal-dependent since they could collect the signalling effects coming from other modules of enhancers or TFs.
The deduced models allow us to predict when, and hopefully where, cells can be relatively activated or repressed with respect to the basal (net activated-repressed), in the presence of two opposing signals (activators vs. repressors). Our study reveals that this can be easily done by using a threshold level between the transcription factors. In Supplementary Materials Section S5, we have proved that each BEWARE operator defines a curve (function) in the [ A ] [ R ] plane which separates that plane in two regions. If the BEWARE is evaluated on concentrations of TFs that are below this curve, then it will predict an expression of P higher than the basal rate (and lower for concentrations above the curve). We have called these concentrations of transcription factors “Activation/Repression regions”, and they are depicted in Figure 2C. It is important to note that the thresholds (i.e., the interphase between the Activation/Repression regions) depend strongly on the biochemical mechanisms involved in the transcriptional binding process by means of the BEWARE operator used in their determination. We have, in fact, defined the threshold function with n enhancers, f m , l ( [ A ] ; n ) , as the function where the corresponding BEWARE operator fulfils this equation:
BEWARE m ( [ A ] , f m , l ( [ A ] ; n ) , [ R N A P ] ; C ) = BEWARE m ( 0 , 0 , [ R N A P ] ; C )
where C = { [ A ] , [ R ] } c in the case of total cooperativity ( l = t ) and C = { { [ A ] } c A , { [ R ] } c R } in the case of partial cooperativity ( l = p ). That is, the thresholds are the BEWARE operator level curves corresponding to the absence of signalling in the enhancers module, the basal level.
This curve can be used to extrapolate several pieces of transcriptional information from the molecular to the tissular level. f m , l predicts if a cell is activated or repressed by simply checking if the concentration of repressors inside of the cell is higher or lower than the threshold corresponding to its activators concentration—or, in other words, a cell with levels [ A ] and [ R ] of activators and repressors will be
  • relatively activated if [ R ] < f m , l ( [ A ] , n ) ;
  • or relatively repressed if [ R ] > f m , l ( [ A ] , n ) .
Hence, if we know the distribution of TF concentration across the developing tissue, we can deduce the position of cells that will be activated or repressed. The information given by the activation and repression regions, added to the knowledge of the distribution of the transcription factors across the tissue, allows the analysis of how the activation profile is distributed spatially. We have shown in Proposition 2, Supplementary Materials Section S5, that f m , l ( [ A ] ; n ) is a strictly increasing function. This assures that the transcriptional control of an opposing-gradient is enough to establish two simple (and connected) regions in the tissue, where the cells are either activated or repressed (see Figure 2). This pattern of activated-or-repressed cells is essential in the development of several biological systems—for instance, Drosophila’s wing development by the target genes dpp and ptc [17,53], among others. The study of the tissular activation/repression patterns has been extensively studied by experimental methods, using modifications of several biochemical properties of the development system. Variations such as modification of TFBSs affinities, or even the number of enhancers, have been shown to have important effects in the formation of these patterns [15]. Motivated by this line of experiments, we have used the threshold definition in order to theoretically reproduce these variations. Our study has been able to extract specific information of how these biochemical variations alter the regions of activated cells (see Section 2.2 and Section 2.3 for more details). This proves the potential of the threshold function as an analytic tool in the modelling of transcription processes. The information, given in terms of concentrations of transcription factors inside of the cell, now can be decoded at a tissular level outside the cell via the analysis of f m , l ( [ A ] ; n ) . We will explore the importance of this detail in the next section, comparing with the experimental evidence obtained in [15].
It is important to mention that the monotonicity of f m , l does not necessarily correspond to the monotonicity of the BEWARE operator with respect to the TFs. Our analysis predicts concentration regions where the BEWARE operators involving total cooperation can exhibit behaviours of “inverse logic” such as increases in the transcription rate where there is an increase in the concentration of repressors (or decreases where there is an increase in the concentration of activators). We refer to this effect as the “pull-effect”, and it takes place when total cooperativity is very strong. Indeed, the analysis performed allows us to describe in great detail when these effects can be found (see Supplementary Materials Section S4 for mathematical analysis and graphical explanation of these effects).

2.2. Sensitivity Analysis of the Threshold Functions: Elasticity

In order to test out the change in the activation/repression regions, we have analysed the behaviour of the threshold function f m , l under the following perturbations:
  • proportional reduction in affinity for the enhancers:
    K A η K A and K R η K R being η 1 ,
  • decrease in the number of available enhancers:
    n n 1 ,
Equation (6) represents the fact that an increase in the perturbation parameter η corresponds to lower affinities between the TFs and the enhancers, remembering that K * are dissociation constants. Our study predicts that the response obtained by these perturbations are closely related in the case of stimulated and recruitment operators. Although both alterations interfere with the action of activators and repressors in the same way, it is surprising that, in all these models, the response to perturbations (6) and (7) is qualitatively predicted by the elasticity of the function f m , l ,
ϵ m , l ( [ A ] ; n ) = [ A ] f m , l ( [ A ] ; n ) f m , l ( [ A ] ; n ) Δ f m , l ( [ A ] ; n ) / f m , l ( [ A ] ; n ) Δ [ A ] / [ A ] ,
where Δ f m , l ( [ A ] ; n ) = f m , l ( [ A ] + Δ [ A ] ; n ) f m , l ( [ A ] ; n ) . ϵ is a quantity usually used in Economics in order to measure a system’s responses to proportional perturbations [76]. It is also known as a condition number in Numerical Analysis [77]. This index has also been introduced in biological contexts, for instance in Ecology [66].
In Supplementary Materials Lemma S4 and Corollary S1, Section S6, we have proven the direct relationship between elasticity and the threshold response to perturbations (6) and (7). These results are resumed in the following expression:
s i g n δ f m , l δ η ( [ A ] ; n ) η = 1 = s i g n f m , l ( [ A ] ; n 1 ) f m , l ( [ A ] ; n ) = s i g n 1 ϵ m , l ( [ A ] ; n ) .
This is applicable to recruitment and stimulated operators ( m = r , s ) in their total and partial cooperativity versions ( l = t , p ). Please note that both identities in (9) show that the Activation/Repression thresholds under signalling interferences (6) and (7) react in the same qualitative manner. For example, the thresholds:
  • will decrease in the elastic regime, that is, if ϵ m , l ( [ A ] ; n ) > 1 ;
  • will not be modified in the unit elastic regime, that is, if ϵ m , l ( [ A ] ; n ) = 1 ;
  • will increase in the inelastic regime, that is, if ϵ m , l ( [ A ] ; n ) < 1 .
We can interpret this result in terms of a loss of competitive advantage between the transcription factors. Binding cooperativity mechanisms between TFs are the clearest example of this competitive advantage. If TFs cooperate in their binding process, this cooperation constitutes an advantage for such TFs, and this advantage is clearly amplified in the presence of
  • high affinity enhancers because the first binding required for cooperativity is more likely to occur,
  • a high number of enhancers because they allow improvement of TF’s affinity by cooperativity.
Thus, perturbations (6) and (7) are clearly limiting the advantage given by cooperativity. This is more obvious in the case of (7) since, in the limit case, n = 1 , cooperativity can not operate at all. As we will see in next section, binding cooperativity is not the only advantage we can detect by means of the elasticity.
The inelastic case is interpreted as a situation where repressors “lose their advantage” over the activators, since a global decrease in affinity or number of enhancers gives rise to an increase in the threshold, as seen in (9), and consequently an increase of the activation region (see Figure 2). Remember that, given a fixed activator concentration, an increase of the threshold function implies that an increase of the repressor concentration is needed in order to get the same transcription rate as the basal. That is to say, the cells will enter in the activated state (i.e., transcription rates higher than the basal) “more easily” in the new situation because the repressors’ advantage is not as effective as before. The same goes for the elastic case. In this case, the activators seem to lose the benefit of the advantage. Here, (9) shows that the thresholds decrease, so the activation region decreases because of perturbations (6) and (7). The unit elastic regime can be seen as a stable situation, where none of the perturbations modify the threshold (i.e., both activators and repressors have the same transcriptional advantage and perturbations affect both in the same way). Please note also that the same reasoning can be applied to proportional decrease in activator and repressor concentrations because of definition (8). In Figure 3, these qualitative behaviours have been illustrated using a stimulated BEWARE operator where three sets of parameters have been chosen to represent all these situations.
With the threshold definition (5), obviously f m , l ( [ A ] = 0 ; n ) = 0 holds true. Thus, the elasticity coefficient has a very simple geometrical interpretation: the comparison between the slope of the tangent line at the point [ A ] , f m , l ( [ A ] ; n ) and the slope of the secant line intersecting the threshold curve at the origin ( 0 , f m , l ( 0 ; n ) ) and at the point ( [ A ] , f m , l ( [ A ] ; n ) ) . That is,
ϵ m , l ( [ A ] ; n ) > 1 f m , l ( [ A ] ; n ) > f m , l ( [ A ] ; n ) f m , l ( 0 ; n ) [ A ] 0 .
This expression also tells us that, if f m , l is convex, then it is elastic. In the same way, the concavity of f m , l implies being inelastic. Nevertheless, as we can see in Figure 3, this geometrical interpretation of ϵ can not always be easily recognised at first glance.
Although Hill models have been deduced in this work from the stimulated and recruitment BEWARE operators by an extreme cooperativity hypothesis, these types of models do not inherit relationships (9). Thus, this conceptual streamlining loses the original relationship between perturbations (6) and (7) and the elasticity function. Indeed, in Supplementary Materials Section S7, it is proven that Hill threshold functions exist, but they are always straight lines whose slope can change with the number of enhancers.

2.3. Calculations: Recruitment, Stimulated and Hill BEWARE Operators: Comparative Analysis

This section is devoted to comparing the effects of the perturbations (6) and (7) on the activation/repression regions in terms of the thresholds predicted by the different BEWARE operators considered. This study shows the existence of two main factors that are seen to determine the model’s response: the first one is the binding cooperativities involved in the TF binding process and the second is how the transcriptional effects of activators and repressors are modelled.
In order to understand the cooperativity effects between TFs in their binding process, we have divided the partial cooperative case in two extreme scenarios: activator cooperative ( c A > 1 , c R = 1 ) and repressor cooperative ( c A = 1 , c R > 1 ) cases. Then, we have included the null and total cooperative case in the same scenario, since the only difference between the operators is the value of c (i.e., c = 1 for the null cooperative case and c > 1 for the total cooperative case). This allows us to identify more easily cases where the activators or repressors should lose “competitive advantages” (discussed in Section 2.2) or not. Particular cases with partial cooperativity between activators and repressors should be estimated numerically.
We analysed stimulated and recruitment BEWARE operators with each kind of cooperativity, obtaining interesting (but unexpected) results which are summarised in Table 1 (detailed proofs can be found in Supplementary Materials Sections S6.1 and S6.2).

3. Results

From the modelling point of view, the expected behaviour in relation to perturbations (6) and (7) should be:
  • Null/Total cooperative case: The activators and repressors lose no competitive advantage in the the null case or lose exactly the same amount of competitive advantage in the total cooperativity case. Hence, the threshold function should not vary (unit elastic case);
  • Activator cooperative case: The activators lose that competitive advantage over the repressors. Hence, the threshold function should decrease (elastic case).
  • Repressor cooperative case: The repressors now lose the competitive advantage over the activators. Hence, the threshold function should increase (inelastic case).
Estimations exhibited in Section 2.3 show that the Recruitment BEWARE operator behaves as expected, that is, the elasticity of the thresholds determined by these operators are proven to be ϵ r , p > 1 in the Activator Cooperative case, ϵ r , p < 1 in the Repressor Cooperative case, and ϵ r , t = 1 in the case of null/total cooperativity.
In the case of the Stimulated BEWARE operators, the analysis of the elasticity variable (at n = 2 ) is more entangled as can be seen in Table 1. Here, the elasticity value is related not only to the cooperativities but also to other parameters that determine each micro-state transcription rate. In fact, we can observe that, regardless of the (binding) cooperativity considered, we can get elastic, inelastic or unit-elastic situations depending on certain relationships between the parameters involved in the modelling.
Let us explain this conclusion revealed by elasticity. Remembering that the basal transcriptional level, corresponding to the transcriptional rate of an empty micro-state, in this approach is t r , = r b a s . Here, represents an enhancer with non bound TFs. The transcriptional rates associated with micro-states with a single bound activator/repressor, as can be seen in Supplementary Materials Section S1, are t r A , = r b a s + e ˜ ν m a x ( 2 ) = r b a s + ν m a x ( 1 ) and t r R , = r ˜ r b a s , respectively. If the linkage of a second TF of the same family occurs, then the new transcription rates become t r A , A = r b a s + ν m a x ( 2 ) = r b a s + ν m a x ( 1 ) / e ˜ and t r R , R = r ˜ 2 r b a s .
Let us observe that the ratio in the variation of the transcriptional rate due to the binding of a repressor is constant regardless of the existence of other already bound repressor, that is,
t r R , t r , = t r R , t r R , R = r ˜ .
However, in the case of activators, the analogous rates depend on the values of e ˜ , ν m a x ( 1 ) and r b a s
t r A , t r , = 1 + ν m a x ( 1 ) r b a s a n d t r A , t r A , A = r b a s + ν m a x ( 1 ) e ˜ r b a s + ν m a x ( 1 )
and can vary from the first to the second binding protein. Then, the existence of the second enhancer implies
(i)
a transcriptional advantage for repressors if t r A , t r , > t r A , t r A , A , since in that case comparatively the binding of a second activator is less effective than the binding of a second repressor;
(ii)
a transcriptional advantage for activators if t r A , t r , < t r A , t r A , A , since in that case comparatively the binding of a second activator is more effective than the binding of a second repressor;
(iii)
no advantage for any TF when t r A , t r , = t r A , t r A , A because in that case the binding of a second TF is equally as effective as the first bound TF.
It is easy to check that (i), (ii) and (iii) directly correspond to the elastic, inelastic or unit-elastic cases under null/total cooperativity determined in Table 1. That is, in the absence of binding cooperativity, the advantage that elasticity demonstrates is related to the possibility that the functioning of two adjacent activators can be more/less effective together than separately. Indeed, this “functional (anti-)cooperativity” mechanism is not new in literature; it was already proposed in [40]. However, with the elasticity, we are able to analyse some specific cases of the models proposed in [40], and the values that the elasticity function takes reveal inconsistencies with the modelling guidelines proposed in that work.
In Table S4, we introduced some identifications/choices that can be made in order to fit the model proposed in [40] into the expression of a stimulated BEWARE operator such as the ones employed in this work. This allowed us to apply, in some specific cases, the analytical tools we have developed to the model proposed in [40]. The results of this analysis (Supplementary Materials Section S6.2) are summarised in the middle row of Table 1, that is, the estimates of the elasticity corresponding to the thresholds deduced from stimulated BEWARE operators with two enhancers ( n = 2 ). This analysis shows that there are acceptable parameter values for the model that give us inelastic threshold functions that are fully compatible with the non-existence of “functional cooperativity” ( ε A = ε R = 1 according to notation in [40]). In Figure 3, the inelastic threshold function (continuous purple line) has been obtained using some of these parameters (see values in Table S4, Supplementary Materials Section S8).
In the same figure, we can see how perturbed thresholds vary according to the elasticity being less than one. On the other hand, this model was deduced without considering binding or functional cooperativity.
Thus, our analysis shows that the model in [40], deduced from non “functional” and non “binding” affinity assumptions, can have an inelastic threshold. In addition, it reacts as it has been proven to do under perturbations (6) and (7). This can be seen in Figure 3. Our argument supports the notion that the asymmetric approach that has been adopted in stimulated operators for modelling the activators and repressors functionalities is causing this effect. The elasticity can also be estimated in the presence of partial binding cooperativity, and we have seen that, in both cases, elasticity is able to balance the effects of both cooperativities.
This unexpected effect needs to be taken into account in analyses such as those developed in [53] because, if used, it can seriously alter the conclusions. This example is one reason why we were interested in performing sensitivity analysis in this work.
Regarding the Hill versions of the Recruitment and Stimulated operators, we can say that ϵ h = 1 for any kind of cooperativity and operator. However, it is important to note that identities in (9) are not valid in this framework. That is, the extreme cooperativity approach used to get the Hill type models is incompatible with the information that elasticity provides.

4. Discussion

This work explores the applicability of classical thermodynamic modelling to the interpretation of experimental results observed at a tissular level in terms of which biochemical mechanisms are involved in transcriptional control. The thermodynamic approach allows us to represent the transcriptional control exerted by opposing transcription factors (TFs) over a particular set of TFBS (enhancers), considering alternative enhancer configurations and biochemical mechanisms. Using these models, we are able to show that it is possible to theoretically predict activated/repressed tissular regions, that is, cells which are relatively activated/repressed.
One of the possible uses of this type of modelling is to theoretically reproduce experimental procedures which compare the effects of given enhancer perturbations on those activated/repressed celular regions. Specifically, we have considered two enhancer module perturbations performed in [15]: one reduces the number of enhancers on the module and the other modifies the enhancers’ affinity. In the theoretical framework, both of these experiments are the equivalent of performing a sensitivity analysis because we are comparing the effects of the analog enhancer perturbations on the activation/repression regions. Testing the experimental results against this kind of theoretical predictions can give insights into the mechanisms involved in the transcriptional control of the considered gene. This has already been done in paradigmatic systems such as Hh target genes [53]. However, the reliability of this comparison depends on how accurately the biochemical mechanisms are represented by the thermodynamic model. In this work, we have analysed this question in models involving two different TF transcriptional control mechanisms which have been well established in previous literature: the recruitment and stimulated mechanisms.
The first conclusion made from our analysis is that, despite the fact that the perturbations proposed in [15] are quite different, the variations of the activation/repression regions due to both perturbations are qualitatively the same. In both stimulated and recruitment approaches, we have proven that the theoretical response of the activation/repression regions depends on the same variable: the elasticity of the threshold function. At first, it seems surprising that the same variable would determine the cause–effect behavioural tendency of such different perturbations. Nevertheless, we must keep in mind that both perturbations are interfering with the transcriptional control. More specifically, they interfere with those mechanisms in which some bound transcription factors assist other transcription factors, usually called cooperativity mechanisms. The elasticity variable is significant in this framework because it is able to distill the effects of the multiple biochemical hypotheses considered.
Although the enhancer perturbations affect all TFs in the same way, the activation/repression variations occur as a consequence of the loss of some functional advantage that some of the TFs had originally. In this point, we have determined the advantages involved in recruitment and stimulated models.
Through analysis of the recruitment operators, we find that, in this case, their elasticity is mainly determined by the binding cooperativity relations between TFs. However, in the case of the stimulated approach, in addition to the binding cooperativity, functional cooperativity is included in the modelling. That is, TF proteins are also allowed to cooperate or anti-cooperate in their transcriptional control function. Indeed, by using elasticity, we have detected that the modelling can (artificially and accidentally) include this functional cooperativity effect. Care must be taken with the parameters to avoid this unexpected inclusion of functional cooperativity effects. This is a particular example where elasticity is capable of balancing the effects of different mechanisms and funnel them into a single variable.
Although Hill type models can also be deduced from the stimulated/recruitment BEWARE operators, we have proven that they alter the response to this sensitivity analysis.
It is important to remark that we can also discuss the applicability of stimulated or Hill models in the interpretation of the results in [15]. That is, we know by Table 1 that:
  • The Hill model fails to reproduce the experimental evidence. This is because the strong cooperativity assumption leads to an elastic threshold, which is not compatible with the results obtained in [15].
  • The stimulated model could also fail, depending on the selection of the model parameters.
The last point notes another interesting use of the elasticity function. Given the experimental evidence, it could be used to restrict the Stimulated model parameters space, using the elasticity function as a helper in the parameter fitting process.
The possibility of performing sensitivity analysis gives new perspectives and opens new opportunities for the use of thermodynamic modelling, especially in the study of the mechanisms involved in signal interpretation. In addition, as new experimental procedures become available, the reevaluation of classical models becomes possible. This is not only relevant per se, but also because these basic models are intrinsically involved in the interpretation of a wide range of phenomena such as in epigenetic control [78], protein dispersion [60] or transcriptional bursting [79,80,81].
.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/math10132169/s1. See [82,83,84,85,86].

Author Contributions

Conceptualization, M.C. and Ó.S.; Investigation, M.C. and Ó.S.; Software, M.C.; Writing—original draft, M.C. and Ó.S.; Writing—review & editing, Ó.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the MINECO-Feder (Spanish Government) research Grant Nos. FPI2015/074837 and RTI2018-098850-B-100, and the Consejería de Economía, Innovación, Ciencia y Empleo, Junta de Andalucía (Andalucía Government) project PY18-RT-2422, and Junta de Andalucía, Prog. Operat. FEDER Andalucía: A-FQM-311-UGR18 & B-FQM-580-UGR20.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Ó.S. would like to thank R. Biewlanski for language editing and proofreading.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations and nomenclature are used in this manuscript:
BEWARE m Binding Equilibrium Weighted Average Rate Expression. The
subindex m is used to denote what version of the model is used
(i.e., m = r denotes Recruitment model, m = s denotes Stimulated
model, m = H r , m = H s denote Hill models).
C Set of transcription factors. Depending on the cooperativity
interactions between the TFs, this set can be arranged on subsets of
cooperativity, characterized by the cooperativity constant c i .
c i Cooperativity constant of a protein i. Depending on the model, this
constant can be the same for all proteins (total cooperativity), in
which case c i = c , or extreme (Hill), in which case c i = c .
CiCubitus Interruptus. Transcription factor involved in the regulation
of Hedgehog target genes.
cis (-regulatory regions)Regions of non-coding DNA which regulate the transcription of
neighboring genes.
DppDecapentaplegic. Protein synthesized from the transcription of the
gen dpp, one of the target genes of Hedgehog.
ϵ m , l Elasticity function of f m , l . It measures the system response to
proportional perturbations in the model parameters.
f m , l Threshold function. It defines the concentration of TFs needed in
order to get a basal transcription in the a BEWARE m model. The
subindex l denotes what kind of cooperativity is applied between
the TFs (i.e., l = t denotes Total cooperativity and l = p denotes
partial cooperativity).
HhHedgehog. Morphogen involved in the development of
Drosophila melanogaster.
j i Occupation number of a protein i. It denotes the number of
cis-regulatory regions that are bound by the protein i.
K i Dissociation constant of a protein i. It is related to
the binding affinity of the protein i by its inverse (i.e., the larger K i
is, the lower the binding affinity of the protein i is.)
mRNARNA messenger. Single-stranded segment of RNA that corresponds
to the genetic sequence of a gene.
nNumber of TFBSs (enhancers).
RNAPRNA Polymerase (Pol II). Protein that binds the DNA in a specific
cis-regulatory region (promoter) and starts genetic transcription.
ShhSonic Hedgehog. Morphogen member of the Hh family in vertebrates.
TFTranscription Factor. Protein that binds the DNA in a specific
cis-regulatory region (enhancers; TFBSs) and regulates genetic
transcription.
TFBSsTranscription factor binding sites. Specific DNA cis-regulatory region
that transcription factors bind to in order to regulate the
genetic transcription.
PtcPatched. Transcription factor involved in the regulation of Hedgehog
target genes.

References

  1. Bradshaw, R.A.; Dennis, E.A. Handbook of Cell Signaling; Elsevier Science and Technology: Amsterdam, The Netherlands, 2009; Volume 3. [Google Scholar]
  2. Kicheva, A.; Cohen, M.; Briscoe, M. Developmental pattern formation: Insights from physics and biology. Science 2012, 338, 210–212. [Google Scholar] [CrossRef] [PubMed]
  3. Griffiths, P.; Stotz, K. Genetics and Philosophy. An Introduction; Cambridge University Press: Cambridge, UK, 2013. [Google Scholar]
  4. Lee, T.I.; Young, R.A. Transcriptional Regulation and Its Misregulation in Disease. Cell 2013, 152, 1237–1251. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Perutz, M.F. Mechanisms of cooperativity and allosteric regulation in proteins. Q. Rev. Biophys. 1989, 22, 139–237. [Google Scholar] [CrossRef] [PubMed]
  6. De Las Rivas, J.; Fontanillo, C. Protein-protein interactions essentials: Key concepts to building and analyzing interactome networks. PLoS Comput. Biol. 2010, 6, e1000807. [Google Scholar] [CrossRef] [Green Version]
  7. Ben-Naim, A. Cooperativity in binding of proteins to DNA. J. Chem. Phys. 1997, 107, 10242–10252. [Google Scholar] [CrossRef]
  8. Ben-Naim, A. Cooperativity in binding of proteins to DNA. II. Binding of bacteriophage lambda repressor to the left and right operators. J. Chem. Phys. 1998, 108, 6937–6946. [Google Scholar] [CrossRef]
  9. Ptashne, M.; Gann, A. Transcription activation by recruitment. Nature 1997, 386, 569–577. [Google Scholar] [CrossRef]
  10. Ptashne, M.; Gann, A. Genes and Signals; Cold Spring Harbor Laboratory Press: Woodbury, NY, USA, 2001. [Google Scholar]
  11. Ptashne, M. Regulation of transcription: From lambda to eukaryotes. Trends Biochem. Sci. 2005, 30, 275–279. [Google Scholar] [CrossRef]
  12. Tomaszewska, A.; van Mourik, S.; Blom, J.; Westerhoff, H.V.; Carlberg, C.; Bruggeman, F.J. Tracing the molecular basis of transcriptional dynamics in noisy data by using an experiment-based mathematical model. Nucleic Acids Res. 2015, 43, 153–161. [Google Scholar]
  13. McClure, W. Rate-limiting steps in RNA chain initiation. Proc. Natl. Acad. Sci. USA 1980, 77, 5634–5638. [Google Scholar] [CrossRef] [Green Version]
  14. Na, D.; Lee, S.; Lee, D. Mathematical modelling of translation initiation for the estimation of its efficiency to computationally design mRNA sequences with desired expression levels in prokaryotes. BMC Syst. Biol. 2010, 4, 71. [Google Scholar] [CrossRef] [Green Version]
  15. Parker, D.S.; White, M.A.; Ramos, A.I.; Cohen, B.A.; Barolo, S. The cis-Regulatory Logic of Hedgehog Gradient Responses: Key Roles for Gli Binding Affinity, Competition, and Cooperativity. Sci. Signal. 2011, 4, ra38. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Torroja, C.; Gorfinkiel, N.; Guerrero, I. Mechanisms of Hedgehog Gradient Formation and Interpretation. J. Neurobiol. 2005, 64, 334–356. [Google Scholar] [CrossRef] [PubMed]
  17. Tabata, T.; Takei, Y. Morphogenes, their identification and regulation. Development 2004, 131, 703–712. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Junker, J.P.; Peterson, K.A.; Nishi, Y.; Mao, J.; McMahon, A.P.; van Oudenaarden, A. A Predictive Model of Bifunctional Transcription Factor Signalling during Embryonic Tissue Patterning. Dev. Cell 2014, 31, 448–460. [Google Scholar] [CrossRef] [Green Version]
  19. Ay, A.; Arnosti, D.N. Mathematical modeling of gene expression: A guide for the perplexed biologist. Crit. Rev. Biochem. Mol. Biol. 2011, 46, 137–151. [Google Scholar] [CrossRef] [Green Version]
  20. de Jong, H. Modeling and simulation of genetic regulatory systems: A literature review. J. Comput. Biol. 2002, 9, 67–103. [Google Scholar] [CrossRef] [Green Version]
  21. García-Ojalvo, J. Physical approaches to the dynamics of genetic circuits: A tutorial. Contemp. Phys. 2011, 52, 439–464. [Google Scholar] [CrossRef]
  22. Chen, W.W.; Niepel, M.; Sorger, P.K. Classic and contemporary approaches to modelling biochemical reactions. Genes Dev. 2010, 17, 1861–1875. [Google Scholar] [CrossRef] [Green Version]
  23. Filkov, V. Identifying gene regulatory networks from gene expression data. In Handbook of Computational Molecular Biology; CRC Press: London, UK, 2005; pp. 708–736. [Google Scholar]
  24. Friedman, N.; Linial, M.; Nachman, I.; Peer, D. Using Bayesian networks to analyse expression data. J. Comput. Biol. 2000, 7, 601–620. [Google Scholar] [CrossRef]
  25. Friedman, N. Inferring cellular networks using probabilistic graphical models. Science 2004, 303, 799–805. [Google Scholar] [CrossRef] [PubMed]
  26. Markowetz, F.; Spang, R. Inferring cellular networks-a review. BMC Bioinform. 2007, 6, 5. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Elowitz, M.B.; Levine, A.J.; Siggia, E.D.; Swain, P.S. Stochastic gene expression in a single cell. Science 2002, 297, 1183–1186. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Milias-Argeitis, A. Elucidation of Genetic Interactions in the Yeast GATA-Factor Network Using Bayesian Model Selection. PLoS Comput. Biol. 2016, 12, e1004784. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Wang, R.S.; Saadatpour, A.; Albert, R. Boolean modelling in systems biology: An overview of methodology and applications. Phys. Biol. 2012, 9, 055001. [Google Scholar] [CrossRef] [Green Version]
  30. Ackers, G.K.; Johnson, A.D.; Shea, M.A. Quantitative model for gene regulation by λ phage repressor. Proc. Natl. Acad. Sci. USA 1982, 79, 1129–1133. [Google Scholar] [CrossRef] [Green Version]
  31. Shea, M.; Ackers, G.K. The OR Control System of Bacteriophage Lambda. A Physical-Chemical Model for Gene Regulation. J. Mol. Biol. 1985, 181, 211–230. [Google Scholar] [CrossRef]
  32. Bintu, L.; Buchler, N.E.; García, H.G.; Gerland, U.; Hwa, T.; Kondev, J.; Phillips, R. Transcriptional regulation by the numbers: Models. Curr. Opin. Genet. Dev. 2005, 15, 116–124. [Google Scholar] [CrossRef] [Green Version]
  33. Buchler, N.E.; Gerland, U.; Hwa, T. On schemes of combinatorial transcription logic. Proc. Natl. Acad. Sci. USA 2003, 100, 5136–5141. [Google Scholar] [CrossRef] [Green Version]
  34. Cohen, M.; Page, K.M.; Perez-Carrasco, R.; Barnes, C.P.; Briscoe, J. A theoretical framework for the regulation of Shh morphogen-controlled gene expression. Development 2014, 141, 3868–3878. [Google Scholar] [CrossRef] [Green Version]
  35. Gilman, A.; Arkin, A.P. Genetic “Code”: Representation and Dynamical Models of Genetic Components and Networks. Annu. Rev. Genom. Hum. Genet. 2002, 3, 341–369. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Segal, E.; Widom, J. From DNA Sequence to Transcriptional Behaviour: A Quantitative Approach. Nat. Rev. Gen. 2009, 10, 443–456. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Kininis, M.; Isaacs, G.D.; Core, L.J.; Hah, N.; Kraus, W.L. Postrecruitment regulation of RNA polymerase II directs rapid signalling responses at the promoters of estrogen target genes. Mol. Cel. Biol. 2009, 29, 1123–1133. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Fakhouri, W.D.; Ay, A.; Sayal, R.; Dresch, J.; Dayringer, E.; Arnosti, D.N. Deciphering a transcriptional regulatory code: Modeling short-range repression in the Drosophila embryo. Mol. Syst. Biol. 2010, 6, 34. [Google Scholar] [CrossRef] [PubMed]
  39. He, X.; Samee, M.A.; Blatti, C.; Sinha, S. Thermodynamics-based models of transcriptional regulation by enhancers: The roles of synergistic activation, cooperative binding and short-range repression. PLoS Comput. Biol. 2010, 6, e1000935. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  40. Keller, A.D. Model Genetic Circuits Encoding Autoregulatory Transcription Factors. J. Theor. Biol. 1995, 172, 169–185. [Google Scholar] [CrossRef]
  41. Lai, K.; Robertson, M.J.; Schaffer, D.V. The Sonic Hedgehog Signalling System as a Bistable Genetic Switch. Biophys. J. 2004, 86, 2748–2757. [Google Scholar] [CrossRef] [Green Version]
  42. Hill, A. The combinations of haemoglobin with oxygen and with carbon monoxide. J. Physiol. 1910, 40, 4–7. [Google Scholar] [CrossRef] [Green Version]
  43. Hill, T. Cooperativity Theory in Biochemistry; Springer: New York, NY, USA, 1985. [Google Scholar]
  44. von Dassow, G.; Meir, E.; Munro, E.M.; Odell, G.M. The segment polarity network is a robust developmental module. Nature 2000, 406, 188–192. [Google Scholar] [CrossRef]
  45. Cerone, L.; Neufeld, Z. Differential Gene Expression Regulated by Oscillatory Transcription Factors. PLoS ONE 2012, 7, e30283. [Google Scholar] [CrossRef]
  46. Aguilar-Hidalgo, D.; Dominguez-Cejudo, M.A.; Amore, G.; Brockmann, A.; Lemos, M.C.; Córdoba, A.; Casares, F. A Hh-driven gene network controls specification, pattern and size of the Drosophila simple eyes. Development 2013, 140, 82–92. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Garcia-Morales, D.; Navarro, T.; Iannini, A.; Pereira, P.S.; Míguez, D.G.; Casares, F. Dynamic Hh signalling can generate temporal information during tissue patterning. Development 2019, 146, dev176933. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Sherman, M.S.; Cohen, B.A. Thermodynamic State Ensemble Models of cis-Regulation. PLoS Comput. Biol. 2012, 8, e1002407. [Google Scholar]
  49. Segal, E.; Raveh-Sadka, T.; Schroeder, M.; Unnerstall, U.; Gaul, U. Predicting expression patterns from regulatory sequence in Drosophila segmentation. Nature 2008, 451, 535–540. [Google Scholar] [CrossRef] [PubMed]
  50. Frank, T.D.; Cavadas, M.A.S.; Nguyen, L.K.; Cheong, A. Nonlinear Dynamics in Transcriptional Regulation: Biological Logic Gates. In Nonlinear Dynamics in Biological Systems; SEMA SIMAI Springer Series; Springer International Publishing: Cham, Switzerland, 2016; Volume 7, pp. 43–62. [Google Scholar]
  51. Mü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. [Google Scholar] [CrossRef]
  52. Niewiadomski, P.; Niedziółka, S.M.; Markiewicz, L.; Uśpieński, T.; Baran, B.; Chojnowska, K. Gli Proteins: Regulation in Development and Cancer. Cells 2019, 8, 147. [Google Scholar] [CrossRef] [Green Version]
  53. Cambón, M.; Sánchez, O. Analysis of the transcriptional logic governing differential spatial expression in Hh target genes. PLoS ONE 2019, 14, e0209349. [Google Scholar] [CrossRef] [Green Version]
  54. Gertz, J.; Siggia, E.D.; Cohen, B.A. Analysis of combinatorial cis-regulation in synthetic and genomic promoters. Nature 2009, 457, 215–218. [Google Scholar] [CrossRef] [Green Version]
  55. Gertz, J.; Cohen, B.A. Environment-specific combinatorial cis-regulation in synthetic promoters. Mol. Syst. Biol. 2009, 5, 244. [Google Scholar] [CrossRef]
  56. Audun, B.; Metzler, R.; Sneppen, K. Sensitivity of OR in phage λ. Biophys. J. 2004, 86, 58–66. [Google Scholar]
  57. Law, S.M.; Bellomy, G.R.; Schlax, P.J.; Record, M.T. In vivo thermodynamic analysis of repression with and without looping in lac constructs. J. Mol. Biol. 1993, 230, 161–173. [Google Scholar] [CrossRef] [PubMed]
  58. Zinzen, R.P.; Senger, K.; Levine, K.; Papatsenko, D. Computational models for neurogenic gene expression in the Drosophila embryo. Curr. Biol. 2006, 16, 1358–1365. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Saha, K.; Schaffer, D.V. Signal dynamics in Sonic hedgehog tissue patterning. Development 2006, 133, 889–900. [Google Scholar] [CrossRef] [Green Version]
  60. Verbeni, M.; Sánchez, O.; Mollica, E.; Siegl–Cachedenier, I.; Carlenton, A.; Guerrero, I.; Ruiz i Altaba, A.; Soler, J. Morphogenetic action through flux–limited spreading. Phys. Life Rev. 2013, 10, 457–475. [Google Scholar] [CrossRef] [PubMed]
  61. Granek, J.A.; Clarke, N.D. Explicit equilibrium modeling of transcription-factor binding and gene regulation. Genome Biol. 2005, 6, R87. [Google Scholar] [CrossRef] [Green Version]
  62. Simpson-Brose, M.; Treisman, J.; Desplan, C. Synergy between the hunchback and bicoid morphogens is required for anterior patterning in Drosophila. Cell 1994, 78, 855–865. [Google Scholar] [CrossRef]
  63. Szymanski, P.; Levine, M. Multiple modes of dorsal-bHLH transcriptional synergy in the Drosophila embryo. EMBO J. 1995, 14, 2229–2238. [Google Scholar] [CrossRef]
  64. Shearwin, K.E.; Perkins, A.J.; Burr, T.; Hochschild, A.; Egan, J.B. Cooperativity in the long-range gene regulation by the λ cI repressor. Genes Dev. 2004, 18, 344–354. [Google Scholar]
  65. Zhou, X.; Su, Z. tCal: Transcriptional probability calculator using thermodynamic model. Bioinformatics 2008, 24, 2639–2640. [Google Scholar] [CrossRef] [Green Version]
  66. Caswell, H. Matrix Population Models: Construction, Analysis, and Interpretation; Sinauer Associate, Inc.: Sunderland, MA, USA, 2001. [Google Scholar]
  67. Turányi, T.; Tomlin, A.S. Analysis of Kinetic Reaction Mechanisms; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
  68. Rabani, M.; Levin, J.Z.; Fan, L.; Adiconis, X.; Raychowdhury, R.; Garber, M.; Gnirke, A.; Nusbaum, C.; Hacohen, N.; Friedman, N.; et al. Metabolic labeling of RNA uncovers principles of RNA production and degradation dynamics in mammalian cells. Nat. Biotechnol. 2011, 29, 436–442. [Google Scholar] [CrossRef] [Green Version]
  69. Nahmad, M.; Stathopoulos, A. Dynamic Interpretation of Hedgehog Signaling in the Drosophila Wing Disc. PLoS Biol. 2009, 7, e1000202. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  70. Lorberbaum, D.S.; Ramos, A.I.; Peterson, K.A.; Carpenter, B.S.; Parker, D.S.; De, S.; Hillers, L.E.; Blake, V.M.; Nishi, Y.; McFarlane, M.R. An ancient yet flexible cis-regulatory architecture allows localized Hedgehog tuning by patched/Ptch1. eLife 2016, 5, e13550. [Google Scholar] [CrossRef] [PubMed]
  71. Ramos, A.I.; Barolo, S. Low-affinity transcription factor binding sites shape morphogen responses and enhancer evolution. Philos. T. R. Soc. Lon. B 2013, 368, 20130017. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Arnosti, D.N.; Gray, S.; Barolo, S.; Zhou, J.; Levine, M. The gap protein knirps mediates both quenching and direct repression in the Drosophila embryo. EMBO J. 1996, 15, 3659–3666. [Google Scholar] [CrossRef]
  73. Janssens, H.; Hou, S.; Jaeger, J.; Kim, A.R.; Myasnikova, E.; Sharp, D.; Reinitz, J. Quantitative and predictive model of transcriptional control of the Drosophila melanogaster even-skipped gene. Nat. Genet. 2006, 38, 1159–1165. [Google Scholar] [CrossRef]
  74. Arkin, A.; Ross, J.; McAdams, H.H. Stochastic kinetic analysis of developmental pathway bifurcation in phage lambda-infected Escherichia coli cells. Genetics 1998, 149, 1633–1648. [Google Scholar] [CrossRef]
  75. Von Hippel, P.H.; Revzin, A.; Gross, C.A.; Wang, A.C. Non-specific DNA binding of genome regulating proteins as a biological control mechanism: 1. The lac operon: Equilibrium aspects. Proc. Natl. Acad. Sci. USA 1974, 71, 4808–4812. [Google Scholar] [CrossRef] [Green Version]
  76. Sloman, J.; Wride, A. Economics, 7th ed.; Pearson Education Limited: London, UK, 2009. [Google Scholar]
  77. Gautschi, W. Numerical Analysis, 2nd ed.; Springer: New York, NY, USA, 2012. [Google Scholar]
  78. Folguera-Blasco, N.; Perez-Carrasco, R.; Cuyas, E.; Menendez, J.A.; Alarcon, T. A multiscale model of epigenetic heterogeneity-driven cell fate decisicion-making. PLoS Comput. Biol. 2019, 15, e1006592. [Google Scholar] [CrossRef] [Green Version]
  79. Bokes, P.; Lin, Y.T.; Singh, A. High Cooperativity in Negative Feedback can Amplify Noisy Gene Expression. Bull. Math. Biol. 2018, 80, 1871–1899. [Google Scholar] [CrossRef]
  80. Fukaya, T.; Lim, B.; Levine, M. Enhancer Control of Transcriptional Bursting. Cell 2016, 166, 358–368. [Google Scholar] [CrossRef] [Green Version]
  81. Yokoshi, M.; Kawasaki, K.; Cambón, M.; Fukaya, T. Enhancer Control of Transcriptional Bursting. NAR 2022, 50, 92–107. [Google Scholar] [CrossRef]
  82. Érdi, P.; Tóth, J. Mathematical Models of Chemical Reactions: Theory and Applications of Deterministic and Stochastic Models; Manchester University Press: Manchester, UK, 1989. [Google Scholar]
  83. Bisswanger, H. Encyme Kinetics: Principles and Methods; Wiley-VCH: Weinheim, Germany, 2008. [Google Scholar]
  84. Weiss, J.N. The Hill equation revisited: Uses and misuses. FASEB J. 1997, 11, 835–841. [Google Scholar] [CrossRef] [PubMed]
  85. Frank, T.D.; Carmody, A.M.; Kholodenko, B.N. Versatility of Cooperative Transcriptional Activation: A Thermodynamical Modeling Analysis for Greater-Than-Additive and Less-Than-Additive Effects. PLoS ONE 2012, 7, e34439. [Google Scholar] [CrossRef] [PubMed]
  86. Dilão, R. The regulation of gene expression in eukaryotes: Bistability and oscillations in repressilator models. J. Theor. Biol. 2014, 340, 199–208. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Representation of the biochemical mechanisms involved in transcriptional control. Competitive binding process of the TFs (red and purple circles) to n identical enhancers (green). RNA Polymerase (brown oval) binds to the promoter (brown rectangle). Affinities are indicated by black arrows. Alternative biochemical mechanisms are represented with discontinuous arrows: total vs. partial binding cooperativity between TFs (red) and recruitment vs. stimulated approaches for transcriptional control exerted by the TFs (blue). See assumption H2, HS4 and HR4 in Supplementary Materials Section S1 for details.
Figure 1. Representation of the biochemical mechanisms involved in transcriptional control. Competitive binding process of the TFs (red and purple circles) to n identical enhancers (green). RNA Polymerase (brown oval) binds to the promoter (brown rectangle). Affinities are indicated by black arrows. Alternative biochemical mechanisms are represented with discontinuous arrows: total vs. partial binding cooperativity between TFs (red) and recruitment vs. stimulated approaches for transcriptional control exerted by the TFs (blue). See assumption H2, HS4 and HR4 in Supplementary Materials Section S1 for details.
Mathematics 10 02169 g001
Figure 2. Activation-repression threshold and net activated/repressed tissular regions. (A) BEWARE operator representation. Blue, green and red correspond to values of the BEWARE operator, lower than, close to or higher than the basal level. (B) Triangles represent the tissular activator/repressor gradients governing gene transcription. Lower circles show net activated/repressed tissular regions determined by the upper gradients and the BEWARE operator in (A). The limiting position, x 0 corresponds to cells that have TF concentrations determined by the white circled in figure (C); (C) representation of (A,B) in the [A]–[R] plane. The blue region is concentrations leading to BEWARE values below the basal level. The red region is concentrations leading to BEWARE values above the basal level. The green curve is the activation/repression threshold, that is, concentrations leading to the basal level. The grey curve shows the TF concentrations represented by the triangles in (B). The intersection between the green and grey curves, i.e., white circle, determines tissular regions where cells are net activated or repressed according to the BEWARE operator. Parameter values can be found in Supplementary Materials Table S3.
Figure 2. Activation-repression threshold and net activated/repressed tissular regions. (A) BEWARE operator representation. Blue, green and red correspond to values of the BEWARE operator, lower than, close to or higher than the basal level. (B) Triangles represent the tissular activator/repressor gradients governing gene transcription. Lower circles show net activated/repressed tissular regions determined by the upper gradients and the BEWARE operator in (A). The limiting position, x 0 corresponds to cells that have TF concentrations determined by the white circled in figure (C); (C) representation of (A,B) in the [A]–[R] plane. The blue region is concentrations leading to BEWARE values below the basal level. The red region is concentrations leading to BEWARE values above the basal level. The green curve is the activation/repression threshold, that is, concentrations leading to the basal level. The grey curve shows the TF concentrations represented by the triangles in (B). The intersection between the green and grey curves, i.e., white circle, determines tissular regions where cells are net activated or repressed according to the BEWARE operator. Parameter values can be found in Supplementary Materials Table S3.
Mathematics 10 02169 g002
Figure 3. Activation-repression threshold variations under TFs signalling interferences (6) and (7). Thresholds determined by a stimulated BEWARE operator exhibiting different elasticities depending on the values of the model parameters. Orange, green and purple continuous lines are used to represent thresholds with elasticity greater than, equal to or less than one, respectively. In (A) perturbation, (6) has been applied with η = 10 . In (B) perturbation, (7) is where the number of enhancers changes from n = 2 to n = 1 . The thresholds after perturbations are depicted in both cases as dotted lines. Orange and purple arrows show the threshold variations in the elastic and inelastic cases, since in the unit elastic case (green threshold) there is no change, in accordance with Equation (9). Parameter values can be found in Table S4, Supplementary Materials Section S8 and the estimation of the elasticities for these values are in Table 1.
Figure 3. Activation-repression threshold variations under TFs signalling interferences (6) and (7). Thresholds determined by a stimulated BEWARE operator exhibiting different elasticities depending on the values of the model parameters. Orange, green and purple continuous lines are used to represent thresholds with elasticity greater than, equal to or less than one, respectively. In (A) perturbation, (6) has been applied with η = 10 . In (B) perturbation, (7) is where the number of enhancers changes from n = 2 to n = 1 . The thresholds after perturbations are depicted in both cases as dotted lines. Orange and purple arrows show the threshold variations in the elastic and inelastic cases, since in the unit elastic case (green threshold) there is no change, in accordance with Equation (9). Parameter values can be found in Table S4, Supplementary Materials Section S8 and the estimation of the elasticities for these values are in Table 1.
Mathematics 10 02169 g003
Table 1. Analytical estimations of elasticity for thresholds deduced from the BEWARE operators in the global activator/repressor framework. The values t 1 , h 1 , t 2 , h 2 appearing in the case of the stimulated operator with two enhancers are defined in terms of the rest of the model parameters (see Section S6.2 in Supplementary Materials, for explicit definitions).
Table 1. Analytical estimations of elasticity for thresholds deduced from the BEWARE operators in the global activator/repressor framework. The values t 1 , h 1 , t 2 , h 2 appearing in the case of the stimulated operator with two enhancers are defined in terms of the rest of the model parameters (see Section S6.2 in Supplementary Materials, for explicit definitions).
Act Coop.Null/Total Coop.Rep Coop.
Recr. ϵ > 1 ϵ = 1 ϵ < 1
Stim . ( n = 2 ) ϵ < 1 if e ˜ > t 2 ϵ 1 if e ˜ < t 2 & [ A ] h 2 ϵ > 1 if e ˜ < t 2 & [ A ] > h 2 ϵ < 1 if e ˜ > t 1 ϵ = 1 if e ˜ = t 1 ϵ > 1 if e ˜ < t 1 ϵ < 1 if e ˜ > c R t 1 & [ A ] > h 1 ϵ 1 if e ˜ > c R t 1 & [ A ] h 1 ϵ > 1 if e ˜ c R t 1
Hill ϵ = 1 ϵ = 1 ϵ = 1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Cambón, M.; Sánchez, Ó. Thermodynamic Modelling of Transcriptional Control: A Sensitivity Analysis. Mathematics 2022, 10, 2169. https://doi.org/10.3390/math10132169

AMA Style

Cambón M, Sánchez Ó. Thermodynamic Modelling of Transcriptional Control: A Sensitivity Analysis. Mathematics. 2022; 10(13):2169. https://doi.org/10.3390/math10132169

Chicago/Turabian Style

Cambón, Manuel, and Óscar Sánchez. 2022. "Thermodynamic Modelling of Transcriptional Control: A Sensitivity Analysis" Mathematics 10, no. 13: 2169. https://doi.org/10.3390/math10132169

APA Style

Cambón, M., & Sánchez, Ó. (2022). Thermodynamic Modelling of Transcriptional Control: A Sensitivity Analysis. Mathematics, 10(13), 2169. https://doi.org/10.3390/math10132169

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop