Cooperative Binding of Transcription Factors Promotes Bimodal Gene Expression Response

In the present work we extend and analyze the scope of our recently proposed stochastic model for transcriptional regulation, which considers an arbitrarily complex cis-regulatory system using only elementary reactions. Previously, we determined the role of cooperativity on the intrinsic fluctuations of gene expression for activating transcriptional switches, by means of master equation formalism and computer simulation. This model allowed us to distinguish between two cooperative binding mechanisms and, even though the mean expression levels were not affected differently by the acting mechanism, we showed that the associated fluctuations were different. In the present generalized model we include other regulatory functions in addition to those associated to an activator switch. Namely, we introduce repressive regulatory functions and two theoretical mechanisms that account for the biphasic response that some cis-regulatory systems show to the transcription factor concentration. We have also extended our previous master equation formalism in order to include protein production by stochastic translation of mRNA. Furthermore, we examine the graded/binary scenarios in the context of the interaction energy between transcription factors. In this sense, this is the first report to show that the cooperative binding of transcription factors to DNA promotes the “all-or-none” phenomenon observed in eukaryotic systems. In addition, we confirm that gene expression fluctuation levels associated with one of two cooperative binding mechanism never exceed the fluctuation levels of the other.


Introduction
At the transcriptional level gene expression is mainly controlled by the transcription factor (TF) proteins that bind specifically to regulatory binding sites on the DNA [1,2]. TFs influence transcription rates by interacting with other components of the core transcriptional apparatus, including RNA polymerase. Due to the fact that TFs bind to DNA regulatory sites in a stochastic fashion, the transition between states of the cis-regulatory systems (CRS) is a stochastic process. Since the number of TF molecules and the number of regulatory sites are too small, the deterministic assumptions, which are valid in macroscopic systems, fail to describe a mesoscopic system such as this [3]. Therefore, due to the fundamentally random nature of chemical reactions, trajectories of individual cells are noisy and do not follow a smooth deterministic course. It is known that the gene expression response of an individual cell to a regulatory signal may be graded or binary [4][5][6]. In the graded response, the output varies smoothly with the input stimulus, whereas in the binary response, also termed the ''all-or-none'' phenomenon, gene expression response mainly occurs at either low or high levels. In the latter case, the resulting heterogeneous response of an ensemble of cells leads to a bimodal distribution of the protein level. This is a mechanism that can contribute to phenotypic diversity in genetically identical cell populations and is critical for increasing population survival in a fluctuating environment [7]. The bimodal response of gene regulatory networks can arise from closed loops (e.g., a two-gene system whose proteins mutually repress their transcriptional activity) or a single gene (where the gene expression product induces its own expression). These systems present bistability and have been reported previously [8][9][10][11]. Additionally, the ''all-ornone'' gene expression response has also been experimentally observed in some eukaryotic systems that do not involve bistability [4,5,[12][13][14], where gene expression often occurs in stochastic bursts. This suggests that the binary responses observed in inducible gene expression could be explained by fluctuations in the binding of TFs to DNA [6,15].
Contrariwise to prokaryotic RNA polymerases, eukaryotic polymerases require the prior assembly of general TFs at the typical eukaryotic promoter [16,17]. These factors assemble in a particular order, beginning with the binding of TFIID to the TATA box. The ordered assembly provides several stages at which the initiation of transcription can be regulated [18,19]. Thus, eukaryotic TFs can either facilitate or hinder the assembly of the transcriptional complex. Consequently, it is of paramount importance to contemplate the potential diversity of the CRS architecture and functionality when considering the various known mechanisms by which proteins and DNA interact [20]. However, most of the existing stochastic models for gene regulation are based on transitions between two CRS states (active and inactive) [21][22][23][24][25][26]. Despite their simplicity, these models extract valuable information about gene expression fluctuation. For example, they have illustrated that graded responses arise from fast chemical kinetics, whereas slow kinetics lead to a binary output [3]. Nevertheless, simple models may not be suitable for studying the role of different mechanisms that participate in complex transcriptional regulation processes.
Recently we proposed a mathematical model for transcriptional regulation in cooperative activator switches, which considers a CRS with several regulatory binding sites for a single kind of activator molecule [27]. In this study, by means of the master equation approach, we derived analytical expressions for the first two moments of the steady-state probability distribution for mRNAs and identified two cooperative binding mechanisms [27]: (i) the recruitment mechanism (RM) where the interaction between TFs increases the probability of binding another TF to DNA; (ii) the stabilization mechanism (SM), where the interaction between TFs decreases the unbinding rate of TFs from DNA. These mechanisms affect the fluctuation level in different ways, but not the mean response [27]. In the present paper, we demonstrated what we previously suggested by examination of some regions of the parameters space [27]: that the stabilization cooperative binding mechanism always presents a level of fluctuation greater than or equal to the recruitment mechanism.
Furthermore, in this paper we incorporate two novel generalizations to our previous model: (i) the capacity to understand cooperative mechanisms for repressor or biphasic switches, by considering that bound TFs can repress transcriptional complex formation and modulate transcriptional initiation in different ways [28]; (ii) the inclusion of analytical expressions for the first two moments of the steady-state probability distribution for proteins, enabling contrastation of theoretical and experimental data. In addition, this is the first study to show that cooperative binding plays an important part in determining the transition from graded to binary responses. In this sense, we establish the parameter space regions where each cooperative binding mechanism presents a graded or a binary response. Thus, our findings show that, as well as slow kinetics [3], cooperativity plays a key role in determining the transition from graded to binary responses.

A General Framework for Complex CRS Modeling
Here we present a framework for study models with many states and an arbitrary number of transitions between the different states. This extension of our previous model [27] includes the stochastic production of proteins.
In principle, the CRS states can represent nucleosome organization, DNA loops, TFs bound or unbound to regulatory sites, RNA polymerase binding, etc. Figure 1 depicts a particular outline for this type of complex model, considering eight possibles states, denoted by s~1,2, . . . ,N~8, and fourteen allowed transitions. In general, the CRS can make transitions from a given state s to state r with probability t s,r . Some CRS states are able to synthesize mRNAs at a state-dependent rate, a 1,s . Each mRNA generates proteins, at a constant rate a 2 . Thus the state of the system is specified by three stochastic variables: the chemical state of the CRS s, the number of mRNAs m and the number of proteins n. s,m and n are integers, where m,nw0 and s is 1, . . . ,N. The model also assumes both mRNAs and proteins are degraded at rates c 1 and c 2 , respectively.
Since our model assumes transcriptional regulation as a stochastic process, the theory of stochastic processes is required to analyze the resulting heterogeneous response of an ensemble of cells to a particular signal. Like other authors [21][22][23]26,27,29], we used the master equation approach to study the average gene expression response in the steady state. We can write the probability of finding, at any given time t, the system in the state (s,m,n) as a vector P m,n t ð Þ~P 1,m,n t ð Þ,P 2,m,n t ð Þ, . . . ,P N,m,n t ð Þ ð Þ . The time evolution for this probability is governed by the following master equation: where t s,r is the transition probability per time unit from state r to state s. The last term on the right-hand side of Eq. (1) describes the CRS dynamics, while the others correspond to the production and degradation of mRNAs and proteins. Unlike the master equation for the previous model [27], Eq. (1) has a new random variable n corresponding to proteins and two new terms associated with their production and degradation.

The Steady-state Solution
A time-dependent solution of Eq. (1) is very difficult to obtain even in simpler models. Nevertheless, we are mainly interested in the steady-state solution for mRNA and protein mean levels and their fluctuations. By elaborating on the approach developed in [27], we were able to compute the first two moments of these where P m~P n,s P s,m,n is the marginal probability of the system to have produced m mRNAs, regardless of both the CRS state and the number of proteins for that state, while P n~P m,s P s,m,n is the marginal probability of the system to have n proteins, regardless of both the CRS state and the number of mRNAs for that state. The fluctuations are measured through the corresponding variances, related to the second moments, The summation limits were suppressed for the sake of readability. From now on, every sum over mRNAs or proteins will run from m,n~0 to m,n~?, while the sum over CRS states will be from s~1 to s~N.
Following [27], the moments of jth order can be written in terms of their associated partial moments. Note that the partial moments of order zero are the marginal probabilities for the operator to be in state From Eq. (1) we can derive a set of ordinary differential equations for the time evolution of the partial moments for any j. As there is no feedback, the equations for the partial moments factorize into independent sets of linear equations, which can easily be solved. For j~0,1, and 2 they are From these we can readily find first-order differential equations governing the time evolution of the first moments and variances _ m m m m~{c 1 mz X s a 1,s P s _ n n n n~{c 2 nza 2 m ( ð11Þ _ s s 2 m~{ 2c 1 m 2 zc 1 mz2c 1 m 2 z z X s ½2a 1,s m s za 1,s (1{2m)P s _ s s 2 n~{ 2c 2 n 2 zc 2 nz2a 2 mnza 2 m Equations (11) immediately reduce, in their steady states, to.
where the above Ã denotes the steady-state solution for the random variable. The steady-state solution for the probability vector P Ã corresponds to the normalized eigenvector related to the zero eigenvalue of the CRS transition matrix, TP Ã~0 . From Eqs (12) for the steady-state variances we find where from the last differential equation for the j~2 partial moments, Eqs (10), the second order moment in its steady state mn Ã can be related to the steady-state first-order partial moments of m and n by mn Ã~a Expressions (13)- (18) are general in the sense that they are valid for any CRS, whatever the details of their dynamics. The expressions for mRNA have been previously reported in [27].
Here we incorporate the expressions of mean and standard deviation for proteins, which will allow contrasting models with experiments as many times as experimentalists assess protein levels. The expression for protein fluctuations predicted here does not differ only in an offset from the mRNA fluctuation, as was found in a previous study that considers a many-state CRS [29]. This model assumes that each mRNA generates a burst of proteins, whose size is geometrically distributed. Indeed, our resulting expressions for steady-state fluctuations of mRNAs and proteins, expressed in the form of normalized variance, conform to the general equation described previously by Paulsson [24], but with a more complicated term for the activation-inactivation transitions. Thus, our results expand upon previous studies that were either limited to the modeling of promoter state transitions as a two-state on/off switch [3,21,24,26] or which excluded translation when more than two promoter states were modeled [27,30].

Modeling Genetic Switches
The expressions of the previous subsection are independent of the specific form of the CRS transition matrixT T and of the number of states N. In this section, we will specify the CRS states and the form of the transition matrixT T associated with a particular CRS that is suitable for modeling the transcriptional regulation of switches. For the sake of simplicity, we will consider only eight states (N~8) as sketched in Fig. 1. In order to study the cooperative regulation our model includes three regulatory binding sites for the same TF (N~3), but the generalization to an arbitrary number of sites is straightforward. As in [27], the states s~1,2,3,4 represent states with zero, one, two, and three binding sites occupied by TFs, respectively. The states s §5 correspond to transcriptional preinitiation complex formation, where all components required for transcription are assembled in the CRS. For simplicity, we consider that TFs do not bind or unbind after the formation of the preinitiation complex; the allowed transitions between the CRS states are indicated by arrows in Fig. 1. Once the core transcriptional apparatus is formed, the synthesis of one mRNA copy begins at rate a 1,s . Each mRNA generates proteins at a constant rate a 2 . Our model also assumes that both mRNAs and proteins are linearly degraded at rates c 1 and c 2 , respectively. In the model we can distinguish four regulatory layers. Layer I corresponds to CRS dynamics of TF binding to/unbinding from DNA, layer II corresponds to preinitiation complex formation, layer III corresponds to mRNA production/degradation, while layer IV corresponds to protein production/degradation.
In order to obtain the explicit expressions of the steady-state solutions in terms of the parameters of the system, we need to specify the CRS transition matrixT T. The TFs can bind to regulatory sites with a probability proportional to TF concentration c, following the law of mass action for elementary reactions. Thus, the transition probabilities t 12~c k 12 ,t 23~c k 23 and t 34~c k 34 , while transition rates t ij to and from other states of the operator are denoted simply as k ij . In this case the transition matrixT T can be written as^T whose associated steady-state solutions of the partial probabilities P s Ã involved in Eq. (13) were calculated in [27]. The explicit expression for levels of mRNAs in the steady state is Closed expressions for the variances can also be obtained for Nƒ3 but are too long to be reported here. As in this paper we deal with steady states only, hereafter we write m to denote m Ã . The working hypothesis in our model is that TFs bound to DNA alter the probability of transcriptional complex formation. Consequently, states s~1, . . . ,4 are characterized by different kinetics for the formation of the preinitiation complex. For simplicity, we consider that the sites are functionally identical. The last assumption implies that the model does not distinguish among states with the same number of TFs bound to the regulatory binding sites. Thus, in our model, the states of CRS are more related to the occupancy number rather than to the binding status of each site. This additional simplification reduces the number of states accessible to the CRS and allows us to explore the role of cooperative binding in the noise expression without considering a combinatorial number of states. In this model, with several states able to transcribe, it will be useful to define the transcriptional efficiencies e m for each occupational number i as the rate of mRNA production when there are i TFs bound to DNA, i.e., e m i ð Þ~a 1,iz5 Ã k iz1,iz5 =k iz5,iz1 for i~0,1,2,3. As in the model the regulatory sites are assumed to be functionally identical, we can introduce a relationship between TF binding/unbinding when there is no interaction between the TFs. Thus, if the probability per time unit that a single TF molecule binds to a regulatory site is p, we have k o s,sz1~( N{sz1)p, with s~1,2,3, and u indicates that there is no interaction between TFs. Similarly, unbinding rates are given by k o sz1,s~s q, where q is the probability per time unit that a single TF molecule unbinds from an occupied site.
A further relationship in layer I can be obtained from the principle of detailed balance, which establishes a relationship between the kinetics and the thermodynamic properties of the system [31]. Thus, we will assume that the probability for a TF molecule to bind to a given regulatory site arises from: (i) the free energy of binding a TF to the specific site DG DNA , (ii) the free energy of interaction between TF molecules bound to adjacent sites DG I . Thus, when there is no TF interaction, we have where k o s,sz1 represents the transition rate from state s to state sz1 when there is no interaction between TFs (k o sz1,s represents the rate of reverse transition) and where R is the gas constant and T is the absolute temperature. In general, the TF molecules interact with each other, i.e., DG I ƒ0. If we now assume that each new bound TF interacts with all TFs already bound to the DNA sites, and furthermore, that this energy is the same for all of them, we have where e~e { DG I RT represents the intensity of the interaction between TFs and a s represents the number of interactions, which, because of our assumption, will be a s~s {1 with s~1,2, . . . ,N.
Relationship (22) leaves an extra degree of freedom, because the interaction between TFs can increase the binding rate k s,sz1 , increasing the ability for the recruitment of new TF for DNA binding, or it can diminish the unbinding rate k sz1,s , increasing the stability of the TF bound to DNA. The first case was denoted as the RM, while the second case was denoted as the SM [27]. In order to understand the effect of these cooperativity binding mechanisms on the regulatory response and their associated fluctuations, we will first consider these mechanisms separately. Thus, using relation (22) and the relations for binding/unbinding rates, we obtain for the first mechanism, while for the second mechanism we have Additionally to the two cooperativity binding mechanisms mentioned above, introduced for the first time in [27], we will here consider the case where both mechanisms are acting simultaneously. In this case, we can write the free energy of interaction as DG I~D G RM zDG SM , where DG RM corresponds to the free energy that increases the ability for new TF recruitment for DNA, while DG SM corresponds to the portion of the free energy that diminishes the unbinding rates k sz1,s . Thus, in this more general scenario, we can write the kinetic constants of layer I as where e RM~e RT , noting that e~e RM |e SM . These thermodynamic relationships allow us to write the kinetic parameters of layer I in terms of three parameter p, q and e.
In the next section we will study the transcriptional response of CRS when the TF concentration c is increased. The mean response can be characterized by three parameters: (i) the saturation value (known as V max ), which is defined as lim c?? m(c); (ii) the half-maximum concentration (denoted here by K d ), which is defined as the concentration c at which m(c)~V max =2 (i.e., K d is a root of the polynomial of degree N); (iii) the steepness n H , which is defined as When these definitions are applied to a Hill function H(c)~V max c n H =(c n H zK n H d ), one can determine the three parameters (V max ,K d ,n H ) related to H(c). The above definition allows us to characterize the sigmoidal response given by Eq. (20) analytically, avoiding a nonlinear fitting procedure. Additionally, we characterize the fluctuation around the mean transcript number by the value of standard deviation, given by Eqs. (12) and (15), at c~K d , which is denoted by s max .

Activator, Repressor and Biphasic Switches
In [27] we reported two different cooperative binding mechanisms for activator switches. Here we expand the proposed model to include different types of switches by appropriately setting kinetic rates in layer II and/or in layer III. For example, a repressor switch is obtained if the transcriptional efficiencies e m i ð Þ decrease monotonically with the occupancy number i. This means that, in the example of Fig. 1, a 1 The kinetic parameter values used here are listed in Table 1. For typical experimental conditions, e~6 corresponds to DG I1 1 kcal/mol. This value is similar to the interaction energy between two l-repressor molecules [32] and a bit higher than the free energy associated with the cooperative binding of E2 proteins (DG I *0:7 kcal/mol.) [33]. The binding and unbinding rates of TFs are consistent with the measured values for the lac repressor [34], and for E2 [35], when the TF concentration is given in nM and mM, respectively. Other parameters are assigned plausible but arbitrary values, due to the absence of kinetic information with regard to the other state transitions.
Here we consider activator and repressor switches where the kinetic rates for preinitiation complex formation increase or decrease linearly with the occupancy number, respectively. For the kinetic parameters used in this case the repressor (K d~0 :98) is less sensitive than the activator (K d~0 :25). We also observe that the peak of s m associated with the repressor is slightly smaller than that associated with the activator. Figure 3 illustrates two examples of biphasic switches when the transcriptional efficiency e m i ð Þ does not depend monotonically on the occupancy number i. In Fig. 3A the modulation of the transcriptional efficiency occurs in layer II, while in Fig. 3B the biphasic response to the TF is obtained by the modulation of the rates a 1,s (i.e., layer III). In both cases mean responses are biphasic. Again the mean response does not depend on the cooperative binding mechanism which is acting. In Fig. 3A the fluctuation level, estimated by the standard deviation s m associated with the SM has peaks near the two values of the concentrations where the response is half the maximum, while in the RM case s m presents only one peak. The fluctuation level around the second halfmaximum concentration depends strongly on the acting cooperative binding mechanism. In order to observe the effect of the second type of transcriptional efficiency modulation, we keep the same overall transcription rates by setting k s,sz4~1 :5 for all s, a 1,5~a1,8~0 :01, and a 1,6~a1,7~2 :0. In this case, depicted in Fig. 3B, the mean response decreases and the standard deviation has only one peak with higher amplitude than in the previous case in which the modulation is acting over the kinetic rates related to layer II. We also compare the fluctuation levels associated with these two types of biphasic switches for different transcriptional efficiencies. In this sense, we compute the coefficient of variation CV , as noise measurement (defined as CV~s=m) at the TF concentration c where SmT reaches the maximum, as a function of the overall transcriptional efficiency e m . In the case of a biphasic switch with modulation of the layer II kinetics, different values of e m are obtained by increasing the rates k s,sz4 and keeping a 1,sz4 constant. For a biphasic switch originated by the modulation of layer III kinetics, this is done by increasing the rates a 1,sz4 keeping the kinetic rates k s,sz4 constant. When comparing the respective cooperative binding mechanisms at different transcriptional efficiencies (Fig. 4), we found that the biphasic switch with the latter modulation is always noisier than that where the biphasic response occurs due to the kinetics of layer II.
In order to study how the response of CRS (i.e., the mean and the fluctuations of the mRNA level) depends on the cooperativity parameter e and on the unbinding rate q, we computed three parameters to characterize the mean response and one to characterize the fluctuation around the mean (see subsection  Vertical gray dashed lines are at c~K d . See Table 1  Modeling genetic switches). The binding rate p only affects the dissociation constant K d as reported in [27]. Figure 5 illustrates the activator response behavior of s max , i.e. s m K d ð Þ, as a function of the unbinding rate q (Fig. 5A) and as a function of e (Fig. 5B). In this case it is observed that s max corresponding to the RM does not exceed that associated with the SM. s max decreases sigmoidally with q. The saturation value at low q and the half-maximal q-value increase with e as can be seen in the inset of Fig. 5A. In fact, for the noncooperative case, s max is lower than for the cooperative cases at low and intermediate values of q, but equal at high values of q. On the other hand, s max increases sigmoidally with e (Fig. 5B); again, s max corresponding to the RM does not exceed that associated with the SM, but both saturate to the same value at high values of e. The curves of Fig. 5A and 5B were computed by evaluating the analytic expression for s max , while the symbols were obtained by simulation using the Gillespie method [36]. Figure 5 also illustrates the behavior of the dissociation constant K d and the steepness n H vs. the unbinding rate q (panel C) and e (panel D). As expected, the sensitivity decreases with the unbinding rate but increases with e. On the other hand, the steepness n H depends only on e and not on q or p (data not shown). As expected, n H (blue line) increases with and saturates at 3 at a high value of . A CRS with N~2 activation sites saturates at 2 (data not shown). Thus, in the limit DG I ??, we can recover the Hill function from the expression of the mean response (Eq. 22). A similar behavior is observed for a repressor switch, Fig. 6, with the exception of the steepness n H as a function of e (Fig. 6D). In this case n H decreases with and saturates at 23 at high interaction energy, as expected for a negative regulator. Further differences between Fig. 5 and Fig. 6 are the sensitivity and the fluctuation level. For these parameter values the repressor is less sensitive and noisier than the activator, as we noted in Fig. 2.

Comparing the Cooperative Binding Mechanisms
Results presented up to this point suggest that the SM is associated with a level of noise greater than, or at least equal to, the RM. Now we are interested in determining whether this behavior is a general feature of these mechanisms or if a different scenario can be expected in some regions of the parameter space. In order to address this question we computed the difference between the variances of SM and RM. For the sake of simplicity we considered a switch with two binding sites. Such simplification is sufficient to consider the effects of the binding cooperative mechanisms and to reduce the number of CRS states to six allowing an analytical approach. That is, referring to Fig. 1, we set  Table 1. Panel B parameters are k s,sz4~1 :5 a 1,5~a1,8~0 :01 and a 1,6~a1,7~2 :0, while the rest of the parameters correspond to those in Table 1. doi:10.1371/journal.pone.0044812.g003 k 34~k43~k48~k84~0 , so as to keep only CRS states with s~1,2,3,5,6 and 7. Noticing that the mean values m m Ã do not depend on the acting mechanism, such difference is given by where m m sÃ are solutions to Eq. (15), with the corresponding transition matrixT T, given by Eq. (19) for RM, and by Eq. (20) for SM. If we consider a switch with a 1,s~0 for s~1,2,3 as discussed previously, then the difference between the variances of SM and RM can be written as where the denominator f 0 is a parameter dependent factor, which is the sum of positive terms and consequently it is always positive definite. The difference between the variances of SM and RM can be written as The above expression for Ds 2 is positive for all the parameter space whenever ew1, thus supporting the presumption that follows from our numerical results, i.e., that, for a switch, such as the one depicted in Fig. 1, but with two sites, the fluctuation level associated to the RM never exceeds the fluctuation level associated with the SM.
In live organism, it is more plausible than these cooperative binding mechanisms act simultaneously rather than in an excluding manner, as illustrated for a clearer interpretation. In this context, we also consider some cases where both RM and SM are acting together. Figure 7 depicts the standard deviation s m for an activator switch with the same as Fig. 2, but each fluctuation curve corresponds to different contributions from each mecha- nism. The solid light-gray line corresponds solely to the SM, the dashed light-gray line corresponds to a contribution of 75% from SM and 25% from RM, the dotted black line corresponds to equal contributions from each mechanism, the dashed dark-gray line corresponds to a contribution of 25% from SM and 75% from RM, and the solid dark-gray line corresponds solely to the RM. From this plot, we can observe that fluctuations have a proportional dependence on the mechanism contributions.

Graded and Binary Responses
While the influence of the different mechanisms in which cooperativity can affect the gene expression response is evident from Figs. 5 and 6, their effects are even more dramatic when the steady-state distribution function is studied. In Fig. 8 we compare the recruitment and stabilization mechanisms along several kinetic rates that render the same mean response function using the same kinetics as in the previous activator case (see Table 1), but with e~12. By multiplying all parameters related to a particular regulatory layer by a factor, we alter the fluctuation level, but not the mean response that depends on ratios rather than on individual kinetic rates. The first panel of row A, A1, is the time series of the mRNA number in one cell generated by stochastic simulations using the parameter values of the RM case. The associated histogram (panel A2) shows the number of times a cell shows a given number of mRNAs measured every 10 minutes over a population of 20000 cells. Panels A4 and A5 are the time series and histogram, respectively, of the mRNA number generated by stochastic simulations for the SM case. All time series and histograms in Fig. 8 were obtained using c~K d . For comparison, in panel A3 of Fig. 8 we depict the noise strength Q~s 2 m =m as a function of m r~m =V max , corresponding to the RM case and the SM case. Interestingly, the mechanism acting by stabilization which reduces the unbinding rates, presents a bimodal distribution (panel A5), while the cooperative recruitment mechanism with the same parameters p, q and e, is associated with a unimodal distribution (panel A2). In the panels corresponding to row B the kinetic rates of layer I were amplified by a factor of 10. In this case, the fluctuation levels diminish considerably for both mechanisms. The opposite occurs when the TF binding/unbinding rates decrease (panels of row C). For slow binding/unbinding rates our model predicts that the level of fluctuation increases and the histogram associated with the RM case becomes bimodal. In panels corresponding to row D, the kinetic rates of layer II were amplified by a factor of 10, which does not have much influence over the histograms. Nevertheless, slower rates in this layer promote a higher level of fluctuation, as shown in the panels of row E, in a similar way to slow rates in layer I depicted in the panels of row C. In the panels corresponding to row F, the kinetic rates of layer III were amplified by a factor of 10. At high production and degradation rates, the fluctuation level of the system produces mRNAs in a burst fashion for both mechanisms, which are not very distinguishable in this regime. Contrariwise, when the kinetic rates of this layer decrease by the same factor, the histograms are narrow and unimodal in both cases. We have also noted that the scaling of layer III has opposite consequences to the scaling of layers I and II. The time series and histograms were obtained at m r~0 :5. Nevertheless, the panels in column 3 show that differences between the two cooperative mechanisms are in general greater at low m r (^0:2) and disappear when m r reaches one. Figure 8 illustrates that both binding cooperative mechanisms are able to produce both unimodal and bimodal distributions depending on the kinetic parameters. This feature depends on the relationship between kinetic rates in layers I and III. In fact, bimodal behavior appears when the kinetics of layer I is slower than the kinetics of layer III. This was also observed in simpler models [3]. For a given kinetic relationship between these layers, there exists a region in the parameter space (q,e) related to bimodal distribution. Parameter p does not affect this distribution feature (data not shown). Figure 9 shows the bimodal regime for both mechanisms is in region I (low q), while the unimodal regime is in region II (high q). The region denoted by I-II corresponds to a region in which the unimodal regime of the RM and the bimodal regime of the SM coexist. Interestingly, the interface between these two regions depends on the acting cooperative binding mechanism and when the RM is acting, lower values of q are required to get a unimodal distribution. In Figs. 9B-9E we can see the dynamic behavior and histogram for q{e values indicated by a star in the phase diagram (e~10, q~0:6, c~K d~0 :13 and other values are the same as in Table 1 for the activator). Thus, bimodality can also be a consequence of the cooperative binding mechanism. But at sufficiently lower unbinding rates q, cooperative binding is not necessary to reach bimodal response.

Discussion
Despite the rich variety of gene regulatory mechanisms acting at the transcriptional level [16,17,20,28], most models consider only one or two states for the CRS. These models approximate the transcriptional control by using a regulatory expression function (Hill function in [11,22,37,38] or an ad-hoc function to fit the model to the experimental data in [29,30]).
We have shown that cooperative regulatory function can be derived from a model based on the law of mass action for elementary reactions [27], which allows understanding the consequences of TF cooperative interactions from first principles. For example, we have shown that response steepness depends on the energy involved in the interaction between TFs. However, this analysis was restricted to activator switches. Consequently, in this study we have generalized our previous approach in order to model repressor and biphasic switches. In our model, the basic components of the CRS can be arranged in different ways to modulate gene expression in response to a given signal. For example, a repressor molecule bound to DNA can block further assembly by interacting with general factors of the transcriptional complex [28,39]. This aspect can be modeled in our approach, representing a repressor switch with cooperative response. Many features of the response associated with this switch are analogous to others reported previously for cooperative activator switches [27]. For example, as expected, the Hill function with integer exponent is recovered for infinity interaction energy [40,41]. Furthermore, we show that, for switches such as the one depicted in Fig. 1, fluctuation levels associated with the recruitment cooperative binding mechanism never exceed those associated with the stabilization mechanism.
We also show two types of switches related to a biphasic response, namely, the CRS that allows full activation when the regulatory TF occurs within a narrow concentration band. The biphasic response has been reported underlying a variety of mechanisms [42][43][44]. For example, Kruppel in Drosophila acts as an activator at low levels but dimerizes at high concentrations and acts as a repressor in the same binding site [45]. Recently, it was observed that E3f1 had a biphasic response to MYC [44]. Yet another mechanism known as transcriptional interference [46] was reported to respond biphasically [47]. Our model illustrates that biphasic responses can also arise from two other mechanisms: (i) when an intermediate occupancy number of binding sites promotes the formation of the transcriptional complex, while inhibition occurs at low and high binding site occupancy numbers; (ii) when the transcriptional complex has a poor ability for RNAPol recruitment or activation and a consequent low rate of mRNA synthesis at low and high occupancy numbers. The former mechanism appears to be associated with a lower fluctuation level than the latter.
It is commonly accepted that systems that present bistability (i.e., two stable steady states under the same external conditions) are associated with a bimodal response. In this sense, some mathematical models provide examples for that [23,48] Fig. 2A but each fluctuation curve corresponds to different contributions from each mechanism. The solid light-gray line corresponds solely to the SM (i.e., e SM~6 and e RM~0 ), the dashed light-gray line corresponds to 75% and 25% contributions from SM and RM, respectively (i.e., e SM^~3 :83 and e RM^1 :57), the dotted black line corresponds to equal contributions from each mechanism (i.e., e SM~eRM~ffi ffi ffi 6 p ), the dashed dark-gray line corresponds to 25% and 75% contributions from SM and RM, respectively (i.e., e SM^~1 :57 and e RM^3 :83), and the solid dark-gray line corresponds solely to the RM (i.e., e RM~6 and e SM~0 ). doi:10.1371/journal.pone.0044812.g007 sufficiently strong regulation can also constitute a mechanism for bimodality [49]. More recently and from a perspective of population balance, it has been shown that bistability is neither sufficient nor necessary for bimodal distributions in a population [50].
On the other hand, the all-or-none phenomenon has been observed in inducible gene expression and has been attributed to a purely stochastic origin. Several stochastic models of gene expression suggest that fluctuations in the binding/unbinding of TFs to/from DNA can explain both graded and binary responses to inducing stimuli [3,[51][52][53][54]. Pirone and Elston showed that the slow transitions are responsible for binary responses, whereas fast transitions produce graded responses [51]. Even though their model contemplates several regulatory binding sites, they do not consider the effects of cooperative binding on the inducible response. In the context of cooperativity, Sanchez et al. developed a repressor model that includes two regulatory sites [55]. In their model, cooperativity acts by decreasing the unbinding rate and is equivalent to our SM case. These authors found that induced responses change from long-tailed to bimodal distribution when the cooperative factor increases (see Fig. 3C in [55]). When SM is acting, simulation results from our model are in agreement with this previous observation which could be expected because, in this case, cooperativity is slowing CRS transitions. Notably, our model suggests that bimodal distributions are also promoted by the cooperative RM when cooperativity is reflected in binding rate increases, which in turn accelerate CRS transitions. To our knowledge, this has not been previously reported and adds new insight to the origin of bimodality and the effects of cooperative binding on gene expression. In particular, our finding that cooperative binding promotes bimodal distributions could explain the bimodal response observed in a stably integrated NF-AT construct in clones of the Jurkat T-cell line [56]. NF-AT molecules bind cooperatively to DNA as has been reported in [57,58] and the construct employs three tandem copies of the NF-AT-binding site. The phase diagram obtained for our model shows that bimodal distribution can be obtained for high interaction energy between TFs even for high unbinding rates in SM. The unimodal and bimodal phases in the q,e-space are delimited by a cooperative binding mechanism dependent curve. Thus, there is a region in the space parameter (q,e) where SM shows a bimodal response while RM is associated with an unimodal regime.
Summarizing, our results with regard to the stochastic model for gene expression suggest that the gene expression regulatory architecture is measurably reflected in its associated mean response and intrinsic noise profiles.  1 and 4). Histograms that show the number of cells with a given number of mRNAs are also shown (columns 2 and 5). Strength noise Q as a function of m r for both mechanisms (column 3). The above features were studied at different kinetic rates. Row A corresponds to the same kinetics shown in Fig. 2A. Rows B and C, all kinetic rates of layer I were multiplied by a factor of 10 and 0:1, respectively. Rows D and E, all kinetic rates of layer II were multiplied by factors of 10 and 0:1, respectively. Rows F and G, all kinetic rates of layer III were multiplied by a factor of 10 and 0:1, respectively. All simulations were performed using Table 1 parameter values, except for e~12 and c~0:14. doi:10.1371/journal.pone.0044812.g008