Exploring the mono-/bistability range of positively autoregulated signaling systems in the presence of competing transcription factor binding sites

Binding of transcription factor (TF) proteins to regulatory DNA sites is key to accurate control of gene expression in response to environmental stimuli. Theoretical modeling of transcription regulation is often focused on a limited set of genes of interest, while binding of the TF to other genomic sites is seldom considered. The total number of TF binding sites (TFBSs) affects the availability of TF protein molecules and sequestration of a TF by TFBSs can promote bistability. For many signaling systems where a graded response is desirable for continuous control over the input range, biochemical parameters of the regulatory proteins need be tuned to avoid bistability. Here we analyze the mono-/bistable parameter range for positively autoregulated two-component systems (TCSs) in the presence of different numbers of competing TFBSs. TCS signaling, one of the major bacterial signaling strategies, couples signal perception with output responses via protein phosphorylation. For bistability, competition for TF proteins by TFBSs lowers the requirement for high fold change of the autoregulated transcription but demands high phosphorylation activities of TCS proteins. We show that bistability can be avoided with a low phosphorylation capacity of TCSs, a high TF affinity for the autoregulated promoter or a low fold change in signaling protein levels upon induction. These may represent general design rules for TCSs to ensure uniform graded responses. Examining the mono-/bistability parameter range allows qualitative prediction of steady-state responses, which are experimentally validated in the E. coli CusRS system.

Introduction Cells rely on accurate control of signaling systems to adapt to ever-changing environments. Appropriate responses to external stimuli often involve transcriptional regulation of gene expression, mediated by various interactions of transcription factors (TFs) with proteins, RNA and regulatory DNA sequences. A set of recurring network motifs, defined as "wiring" or connectivity patterns of TF interactions, can perform specific information-processing functions and constitute the basic building blocks of gene regulatory networks [1]. Architectures of these motifs, as well as the underlying biochemical properties that enable specific cellular functions, have been extensively studied [2][3][4]. These studies often consider only a limited set of transcription factor binding sites (TFBSs) that directly relate to a defined regulon although there is increasing evidence that many seemingly "non-functional" TFBSs exist for many TFs [5,6] and can compete with the autoregulated promoter to impact how individual motifs function [7][8][9]. Here we examine how the relative abundance of TFBSs affects the steady-state response and places constraints on the biochemical parameter range of positively autoregulated signaling systems.
Positive autoregulation is one of the most common motifs in gene regulatory networks. It has often been associated with bistability, i.e., existence of two stable steady states [1,3,10], which can facilitate cell differentiation and bet-hedging under challenging conditions [11]. Approximately 20% of TFs in E. coli activate their own expression according to RegulonDB [12]. More than 10 of these auto-activated E. coli TFs belong to the family of two-component systems (TCSs), one of the major signaling systems in prokaryotes [13,14]. A typical TCS involves a sensor histidine kinase (HK) that responds to environmental stimuli and regulates the phosphorylation level of its cognate response regulator (RR). HKs in prototypical TCSs are usually bifunctional with autokinase/phosphotransferase activities to increase RR phosphorylation and phosphatase activity to decrease phosphorylation ( Fig 1A). Most RRs contain a DNA-binding domain and the phosphorylated RR is the active form that regulates transcription of response genes, many times including genes encoding the HK and RR themselves. Positive autoregulation is more prevalent than negative autoregulation in TCSs, yet bistable responses are not frequently seen while graded responses have been often observed [10,13,15]. The archetype TCSs with bifunctional HKs may have evolved to avoid bistability. It has been argued that bistability, with two cell populations having distinct responses to identical inputs, is often undesirable for many signaling systems [4]. A monostable, uniform and graded response may be favorable for many systems because it allows for a continuity of levels of gene expression or even activation of different sets of genes in response to different signal strengths.
Achieving either the mono-or bistable response demands specific ranges of biochemical parameters, and the mono-/bistability ranges have been explored previously [4,16]. These studies focused on autoregulated pathways, while effects of TF binding to other genomic TFBSs were not considered. This is applicable for many TFs with TF protein abundance in great excess to the number of TFBSs. However, a considerable fraction of TFs has low TF/ TFBS ratios with a large number of TFBSs and high binding demand. In E. coli, the median ratio of TF/TFBS is~10 with many TFs having small ratios ( Fig 1B). Nearly 20% of TFs have a TF/TFBS ratio smaller than 2, insufficient for protein molecules to occupy every DNA binding site, assuming a 2:1 binding stoichiometry that is used by most bacterial TFs. About one third of E. coli TFs have a TF/TFBS ratio smaller than 4, which is far from great excess. The curated TFBS data [12] are traditionally biased toward regulatory sites in proximity to promoters. With recent advances in global TF binding studies revealing many TFBSs in gene coding regions without apparent functions [5,6,17], the fraction of TFs with small TF/TFBS ratios will be even higher and a large fraction of TFs will be considered not to be in great excess to TFBSs [18]. Binding of TFs to these competing TFBSs can sequester the active TF molecules from the autoregulated TF promoter, and protein sequestration may cause ultrasensitivity and bistability [18][19][20][21]. Thus the number of TFBSs can impact the mono-/bistability range and further place constraints on biochemical parameters of the positively autoregulated signaling systems.
In this study, we focused on TCSs and incorporated a DNA binding competition model with autoregulation and RR phosphorylation/TF activation models [4,16,22] to examine the bistability-promoting effect from TFBS competition. TF autoregulation requires a high fold change of transcription for bistability, however, TFBS competition can lower the fold change requirement. We show that bistability is sensitive to the relative binding affinities of competing TFBSs to the TFBS of the autoregulated promoter and the phosphorylation capacity of TCSs. Bistability can be largely avoided with a low fold change, a high TF affinity for the autoregulated promoter and low RR phosphorylation because they all weaken the overall positive feedback. Our model allows qualitative prediction of the mono-/bistability of TCS responses and such prediction was experimentally validated in the E. coli CusRS system.

Model
The general model of a positively autoregulated signaling circuit is illustrated in Fig 1A. Based on biochemical functions and reaction timescales (see details in Methods and S1 Text), the system is decoupled into three modules: TF activation, competitive DNA binding and transcriptional autoregulation. The activation module defines TF activation that can vary for different signaling events such as binding of small molecules, protein-protein interaction or post-translational modification. Here we focus on phosphorylation of RR TFs in TCSs and the terms 'phosphorylation' and 'activation' are used interchangeably in this study. Output of the activation module, the concentration of active/phosphorylated TF molecules R P , is a function of the total concentration of TF R T . Active TF molecules are redistributed among different DNA binding sites via competitive DNA binding, affecting the concentration of free active TF molecules R f . R f is the input of the autoregulation module, determining occupancy of the auto-activated promoter, thus the production rate of total TF levels R T . As shown below, A, B and F represent individual functions of three modules: where k dil is the protein degradation or growth dilution rate. Steady states correspond to intersections of the TF dilution curve and TF production function F (Fig 1C). Bistability is associated with an ultrasensitive response function that has multiple intersection points. Existence of two stable steady states requires existence of an unstable steady state. As demonstrated previously [3,16,23], a necessary condition for bistability, or existence of an unstable steady state, is that the slope of the response function on the log-log dimension is larger than 1. Large values of the slope correspond to ultrasensitive or steep transitions in the response. Such slope has been termed open-loop logarithmic gain (LG) [16]. The overall gain LG 3 is the product of LGs from three modules: where LGs of individual modules are defined as: For simplicity, three modules are considered decoupled so the LG of each module only depends on a limited set of independent parameters. Exploring the parameter space that enables bistability is thus reduced to examining LGs of individual modules that allow the overall gain to be larger than 1. An LG 3 value larger than 1 is a necessary condition for bistability. If the maximum of LG 3 is smaller than 1, the system will be monostable. Details of calculation of LG are included in Methods and S1 Text. Below we examine the autoregulation, DNA binding and phosphorylation modules in order, and discuss effects of individual parameter values on the monostability range.

Monostability requires low transcriptional fold change in autoregulation
The module of transcription autoregulation (Fig 2A) has been extensively investigated [4,16] and transcriptional fold change of the autoregulated promoter greatly impacts the bistability range. Steady-state total TF levels R T range between the basal TF expression level R b and the fully induced level fR b (Fig 2B) where f is the fold change of transcription. Autoregulated transcription depends on promoter occupancy, which is determined by the concentration of free active TF R f , binding affinity K auto and binding cooperativity h: in which all concentration parameters are normalized to dimensionless values by dividing by K auto to reduce complexity, As described previously [4,16], the logarithmic gain LG F can be derived from the above equations and the maximum is: It is apparent that high values of the fold change f or cooperativity h can result in high values of maximal LG F that exceed 1 and promote bistability (Fig 2C). Considering that many RRs bind DNA as a dimer with the cooperativity h = 2, the bistability range is defined by the fold change f. When f is smaller than 9, the maximal LG F will be always smaller than 1, thus monostability favors low values of fold change. Bistability can potentially occur when f is larger than 9. When f ! 1, the maximal LG F is 2. For systems with high values of f, bistability can be prevented if high LG F values can be offset by coupled negative feedback, which reduces LG of the autoregulation module [24], or low LGs of the other two modules.

Monostability range is affected by binding competition between the autoregulated promoter and other TFBSs
It has been shown previously that TFBS competition can lead to bistability in positively autoregulated circuits [18,19]. Here, we quantitatively evaluate how individual parameters of DNA binding competition impact the mono-/bistability range. To model the DNA binding module, we consider an extremely simplified scheme with all other DNA sites being identical with a cooperativity at 2 and a binding affinity K (Fig 3A). All concentration and affinity parameters are normalized by K auto to ensure consistency with the autoregulation module.
The amount of competing TFBSs D � and the relative TF affinity K � determine the extent of competition that the autoregulated promoter faces (Fig 3 and S1 Fig). When there are no other TFBSs (D � = 0), the model is similar to previous analyses that assume a great excess of TF molecules to the number of TFBSs. The concentration of free unbound phosphorylated TF (R f ) is approximately equal to the concentration of phosphorylated TF R P , and the logarithmic gain

PLOS COMPUTATIONAL BIOLOGY
Competing transcription factor binding sites impact the mono-/bistability range LG B of the binding module is always 1 (Fig 3B and S1 Fig). When competing TFBSs are present, binding competition has opposite effects on LG B depending on the relative values of R f � and K � : LG B : The presence of other TFBSs reduces LG B below 1 when R f � is less than K � . This can be understood as these competing TFBSs serving as a "sink" to buffer against accumulation of free active TF molecules. When R f � exceeds K � , the "sink" is gradually getting filled with more than half of TFBSs bound, the accumulation rate of free active TF increases and LG B becomes larger than 1. Once a large fraction of TFBSs is bound (S1 Text), R f � becomes ultrasensitive to R P � and LG B reaches the maximum (Fig 3B and 3C and S1 Fig). The maximum of LG B increases with a larger amount of TFBSs D � (S1 Fig) or lower values of K � , i.e., stronger relative affinity of other TFBSs (Fig 3C), and can reach well above 1.
However, a large maximum of LG B does not necessarily result in a large overall gain or bistability; LG B and LG F are not always synergistic. The combined LG, i.e., the product of LGs from the autoregulation and binding modules, LG F LG B , also depends on whether the maximums of the two modules are achieved at similar R f � levels. Eq (4) indicates the peak LG F value occurring at R f � = f -1/2h , and the parameter space having LG B >1 (R f � >K � ) needs to overlap with the peak LG F region to boost the combined LG. This requires K � values less than or near f -1/2h . Because f -1/2h is smaller than 1 for positively autoregulated systems (f >1), low K � values correspond to relatively strong affinities for competing TFBSs. As shown in Fig 3D, high amounts of competing TFBSs with intermediate TF affinity (D � = 3, K � = 1) yield a similar maximum of LG B as intermediate amounts of TFBSs with strong affinity (D � = 1, K � = 1/3), but the LG B peak of the former (purple dashed line) occurs at a higher R f � level and is not well aligned with the LG F peak of the autoregulation module (black dotted line). Thus, the combined LG remains below 1 while the latter, with high affinity (brown solid line), shows an increase of LG above 1. For TFBSs with a relatively strong affinity, e.g., K � = 0.5, DNA binding competition can promote bistability even if the fold change f is smaller than 9 and a limited parameter space is monostable with the maximal LG smaller than 1 (striped in Fig 3E). With a relatively weak affinity (K � = 2), the bistability range shrinks ( Fig 3E). Even with a large number of competing TFBSs, the system can lessen the bistability-promoting effect of TF sequestration with a high K � value. This is also true with TFBSs having binding cooperativity different from 2, although the mono-/bistability range can vary greatly with different cooperativities (S1 Fig).
For auto-activated signaling systems, a high K � value may be favored to avoid bistability caused by TFBS binding competition and TF sequestration. Because K � is the relative binding affinity between competing TFBSs and the autoregulated promoter, a high K � value can be achieved if the affinity for the autoregulated promoter K auto is among the strongest of all TFBSs. We assessed the relative affinity of autoregulated promoters by inspecting information content of individual binding sites or binding peak intensities in chromatin immunoprecipitation (ChIP) experiments (Fig 4). Information content of individual binding sequences reflects similarity to the overall binding matrix [25] and high information contents often correlate with high affinities [26]. For each of the three selected RRs with well-curated TFBSs, binding site information content of the autoregulated promoter (magenta diamond) is among the top 15% of all TFBSs (Fig 4A). Similar patterns are also observed in ChIP binding intensities for multiple auto-activated E. coli TFs curated in the prokaryotic ChIP database, proChIPdb [27]. For individual TFs, binding peak intensities for positively autoregulated promoters are among the highest of all binding regions across the chromosome (Fig 4B), suggesting that strong affinities for the auto-activated promoters may be selected for. The capability to lessen TFBS competition and avoid bistability with a strong affinity for the auto-activated promoter may be one of the reasons for such selection.

Phosphorylation/dephosphorylation in a typical TCS limit the bistability range
To investigate how signal-regulated activation or phosphorylation of RRs impact the bistability range, we adopted a phosphorylation model described previously [22,28] based on a typical TCS that contains a bifunctional HK ( Fig 5A). The steady-state phosphorylation level of an RR can be approximated with a simple function (see details in Methods and S1 Text) dependent only on the total RR concentration R T , and composite parameters Cp and Ct: where all parameters with an asterisk are normalized by K auto to align with other modules. Signal-regulated activation of RRs is often modeled as an increase of Cp � , leading to an increase of RR phosphorylation level that can be easily approximated with Eq (6). Such approximation relies on the assumption that the RR is in great excess to the HK. We assessed the effects of different RR/HK ratios on RR phosphorylation using the full scheme with typical TCS parameters (S2 Fig). Phosphorylation output with lower RR/HK ratios is not much different from that derived from Eq (6) except when the RR/HK ratio reaches 1. Based on proteomic studies in E. coli [29], a large fraction of TCSs have an RR/HK ratio much higher than 1 (S2 Fig), so the effects of RR/HK ratios on phosphorylation are expected to be small and are not considered here.
As noted previously [22,28], Cp � is the maximal phosphorylation level of the RR. Thus, we term Cp � as the phosphorylation capacity. The RR phosphorylation level R P � increases with the total RR concentration R T � and saturates at Cp � at high R T � levels ( Fig 5B). Correspondingly, the slope of R P � , LG A , is high at low R T � levels when R P � is far from saturation, and  gradually decreases to zero when approaching phosphorylation saturation.
LG A can be derived: LG LG A : The phosphorylation module alone does not promote bistability because LG A will never exceed 1. Instead, phosphorylation and dephosphorylation reactions can limit the bistability range if LG A is small. For example, in the absence of significant TFBS competition (LG B = 1), the LG maximum of the autoregulation module is 2 for an autoregulated TF promoter with a cooperativity of 2 (see Eq (4)). The system can be monostable even with a large fold change f if LG A is smaller than 0.5. The range of LG A smaller than 0.5 is defined by the sum of Cp � and Ct � . Higher values of Cp � result in wider ranges of high LG A values, and thus have a less limiting effect on bistability; vice versa, lower values of Cp � result in systems more likely to be monostable. In general, high Cp � leads to high phosphorylation levels (Fig 5B and 5C) while Ct � is negatively correlated with phosphorylation (S3 Fig). It has been shown that stimuli modulate the autokinase rate k k or phosphatase rate k p or both, which would all be reflected in changes of the phosphorylation capacity Cp � . A high stimulus could boost phosphorylation by increasing the phosphorylation capacity Cp � . Effects of Cp � and Ct � on LG A are illustrated in Fig 5D and  S3 Fig. Parameter ranges that result in a high percentage of RR phosphorylation (red shaded in Fig 5C) give high LG A values close to 1 (white in Fig 5D), while ranges with low LG A values usually have low phosphorylation percentages (gray shaded in Fig 5C and 5D). A high phosphorylation percentage indicates a small difference between R T � and R P � , and thus gives high LG A values according to Eq (7). Biochemically, the phosphatase activity of the bifunctional HK has a negative effect on phosphorylation, causing the phosphorylation level not to rise similarly to R T � and offsetting the positive feedback from transcription. A high percentage of phosphorylation usually suggests low phosphatase activity, thus, less negative impacts on phosphorylation and stronger positive feedback, which promotes bistability.
The full parameter space for mono-/bistability is assessed by multiplying LGs from three modules to examine the overall LG 3 (Fig 6). For bistability with LG 3 >1, the peak LG of the autoregulation and binding modules must align with the region where LG A is relatively large

PLOS COMPUTATIONAL BIOLOGY
Competing transcription factor binding sites impact the mono-/bistability range ( Fig 6A). A large amount of competing TFBSs (orange lines in Fig 6A) promotes bistability, but requires a high level of phosphorylated TF molecules to occupy the DNA sites, thus a high phosphorylation capacity Cp � . As shown in Fig 6B, a high Cp � together with a high fold change f (f = 25) allows a large bistability range. For a system with a low value of fold change (f = 5), bistability can occur via TF sequestration in a limited parameter range (cyan). It demands a high TFBS amount D � , a high TF level R T � and a fairly small TF/TFBS ratio near 2 where the TF is not in great excess. A low phosphorylation capacity Cp � greatly reduces the bistability range (green). Fig 6C summarizes which parameter values of f and Cp � lead to monostability with the maximum of LG 3 smaller than 1. In the absence of other TFBSs (gray shaded), high values of f (f>9) enable bistability at high Cp � values while the system is monostable regardless of f at low Cp � values. The presence of extra TFBSs changes the bistability range that allows for lower f values but requires higher Cp � .
Feedback control in TCSs can involve more complicated schemes with multiple coupled feedbacks [13]. Coupled negative feedback (CNF) has been shown to contribute to robustness, accelerate the response and shape the temporal dynamics [30-32]. Coupled negative feedback can also impact the mono-/bistability range. We consider a simple negative autoregulation scheme with the phosphorylated RR independently binding to a repression site within its own promoter ( Fig 6D). The LG of the coupled regulation module, LG C , is the sum of the LGs of positive and negative autoregulation. With the LG of negative autoregulation always negative (S4 Fig), the LG C is always smaller than the LG F of positive feedback alone. Thus, coupled negative autoregulation always decreases the overall LG and reduces the bistability range (purple line in Fig 6C and S4D Fig). Negative feedback can also occur through an RR activating production of proteins that inhibit RR phosphorylation (dashed line in Fig 6D), such as the E. coli RR PhoP, which activates expression of MgrB, repressing the HK PhoQ [31]. In such a scheme, coupled negative feedback functions to reduce Cp � at the steady state, which also limits the bistability range.

Mono-/bistable response output in the presence of TFBS competition
The mono-bistability range shown in Fig 6C highlights the parameter space that can have a maximal LG 3 larger than 1 and allows bistability. However, a maximal LG 3 larger than 1 is only a necessary condition for bistability. Whether the maximal LG 3 can be achieved at steady states depends on other parameters, such as R b � . Steady states of an autoregulated system correspond to intersections of the RR production function F and the RR dilution curve (gray line in Fig 7A). R b � does not impact LG 3 but will determine the range of RR production levels and affect the steady-state intersections. A sample parameter set with R b � = 2.5 (red line) gives multiple intersections with one unstable and two stable steady states. In contrast, only one intersection or steady state exists for R b � = 0.75 (brown line) with the identical parameter set allowing LG 3 >1. Deriving steady-state levels of R T � and R f � at different signal strengths, i.e., C p � , allows examination of the mono-/bistable responses. Transcriptional responses depend on the steady-state levels of R f � and are modeled as occupancy of RR-regulated promoters with the indicated affinity (Fig 7B and S5 Fig). In the presence of strong TFBS competition (D � = 4, K � = 0.5), the system displays a bistable response (red line in Fig 7B) at high C p � values. Weakening the competition with K � = 2 gives a monostable response (purple) to signals. In general, binding competition reduces the signaling capacity, i.e., "information capacity", of This is experimentally explored in the E. coli CusSR system using increasing numbers of decoy binding sites.

Experimental examination of the mono-/bistability in the CusSR system
Our model reveals how individual parameters shape the mono-/bistability range of TCSs. With accumulation of multi-omics data that allow estimation of some parameter values, it becomes possible to make qualitative predictions about the mono-/bistable response for specific systems. To evaluate the applicability of our model, we experimentally examined the output response of the E. coli CusRS system (Fig 8A). CusS is an HK that senses environmental copper and CusR is an RR transcription factor regulating several genes, including cusCFBA, that are responsible for exporting and detoxifying copper ions [34]. Environmental Cu 2+ concentrations alter the phosphorylation capacity and modulate outputs.
Cellular CusR levels range from~170 to more than 500 molecules per cell in different measurements [29,35,36] while about a dozen CusR binding sites across the genome have been identified [37]. This suggests that CusR is likely in excess to the binding sites and TFBS competition may not play a significant role in the wild-type system (LG B �1). The fold change f has been observed to be small with only a modest increase in CusR protein level upon copper exposure [36]. We estimated the value of f to be~5 based on the PcusR-yfp transcription reporter (S6 Fig). Considering that the autoregulated cusR promoter contains a single binding site [34], the binding cooperativity h will not likely exceed 2 and LG F will be smaller than 1 given f <9. The overall gain will be smaller than 1, thus the wild-type (WT) system will be monostable. The system can become bistable by altering the relative abundance of TFBSs to boost binding competition and increase LG B . This can be achieved by introducing additional CusR binding sites. Exact prediction of DNA binding competition is difficult because of the complexity of CusR bindng modes and the lack of biochemical characterization (see details in S1 Text). However, a number of TFBSs at the same order of magnitude as the number of CusR molecules is likely required to promote bistability.
Output response from a chromosomal PcusC reporter expressing a bright fluorescent protein, mGreenLantern [38] was measured using flow cytometry (S7 Fig). As predicted, in the

PLOS COMPUTATIONAL BIOLOGY
Competing transcription factor binding sites impact the mono-/bistability range absence of extra TFBSs, the WT CusRS system displayed a monostable and graded response to different Cu 2+ concentrations (black colored in Fig 7B and 7C). TFBS competition was exerted by introducing plasmids carrying different numbers of decoy CusR binding sites. Bimodal responses were observed with addition of 96, 120 and 192 decoy binding sites (Fig 7B and S7  Fig). Lowering Cu 2+ concentrations converted the bimodal output to monostable (central and right panels in Fig 7C), e.g., a single uniform distribution of cell response was apparent at intermediate Cu 2+ concentrations with 96 decoy sites. This is consistent with the prediction of our model that lowering the phosphorylation capacity can reduce the overall LG and limit the bistability range.

Discussion
In this study, we investigated how the relative abundance of TFBSs to TF protein molecules impacts the steady-state responses of positively autoregulated TCSs. The particular focus here is on how TFBS competition constrains biochemical parameters in TCS phosphorylation and autoregulation to ensure either mono-or bistable responses. With increasing numbers of TFs discovered to be expressed in limited amounts relative to their binding sites, TF sequestration by competing TFBSs appears to be common for many TFs and can have a wide range of effects on response output, dynamics, noise and timing of gene expression [8,9,18,19,[39][40][41][42]. Such effects are tightly correlated with the abundance of TFBSs and TF proteins, which are subject to considerable variations in cells, e.g., copy number variations in chromosome replication and/or changes in protein levels under different growth conditions. Many signaling systems need to operate under a variety of conditions with different relative abundances of TFBSs.
Here we characterized the mono-/bistability range and showed that a low fold change between basal and maximal transcription rates of the autoregulated promoter, high TF affinity for the

PLOS COMPUTATIONAL BIOLOGY
Competing transcription factor binding sites impact the mono-/bistability range autoregulated promoter and low phosphorylation can ensure monostability for positively autoregulated TCSs.
It has long been known that a large fold change (f >9) is required for bistability [4,16]. TF sequestration by TFBSs enhances ultrasensitivity and positive feedback, lowering the fold change required for bistability. In many autoregulated TCSs, transcription of the TCS genes is under control of at least two promoters, a constitutive one for basal expression and an autoregulated one for maximal induction [13]. The fold change corresponds to the ratio of the strengths of the two promoters. Modulation of either promoter by mutation or stimuli can alter the fold change and the population distribution of the response. In the BvgA TCS, where a monostable and graded response is needed for precise control of the temporal order of gene expression [43], a promoter mutation that lowers the basal expression and increases the fold change has been observed to cause phenotypic bistability [44]. Assessment of bistability requires knowledge of the basal and maximal transcription levels to derive the fold change value. Even though the TF level under high stimulus does not always reach the maximal transcription level, the fold change value can be roughly approximated by the relative change of TF levels in response to high stimulus, which may be readily available in multi-omics datasets or easily measurable. Estimation of the fold change facilitates an initial assessment of the bistability, while the ultimate steady-state response will depend on the relative abundance of TFBSs and phosphorylation status. As demonstrated here for CusR, low fold change coupled with sufficient TF proteins for the TFBSs predicts monostable responses. Many TCSs appears to have high fold change values, e.g., 21  Protein sequestration can cause ultrasensitivity and bistability. In TCSs, sequestration caused by formation of a non-functional "dead-end" complex between the HK and the RR, has been suggested to promote bistability and has been discussed in detail previously [10,50]. The "dead-end" complex provides a mechanism for an implicit feedback to enable bistability even in the absence of transcriptional feedback. The sequestration-based bistability studied here is dependent on the extent of competition between the binding site within the autoregulated promoter and other TFBSs.
Our analyses show that bistability can occur when the affinity of other TFBSs is strong relative to the affinity of the autoregulated promoter. In the presence of a large number of TFBSs, monostable responses can be achieved if these competing TFBSs have a relatively weak affinity in comparison to the autoregulated promoter. Therefore, a much stronger TF affinity for the autoregulated promoter relative to that of other TFBSs may be preferred to alleviate competition and promote monostability. Within available datasets that allow evaluation of the relative affinities, TF affinity for the autoreregulated promoter is always among the strongest of all TFBSs for positively autoregulated TFs, while such a trend is not observed in negatively autoregulated TFs. More global binding data are needed to evaluate this observation. Strong affinity for the autoregulated promoter is also advantageous for a fast transcriptional response, and weak affinity can result in an extremely slow response that does not reach the steady state within a physiologically relevant timeframe [4,9,26,30]. Benefits in both monostability and response speed may contribute to selection for strong affinity of TFs for binding sites within positively autoregulated promoters.
For sequestration-based bistability, high logarithmic gain or sensitivity occurs at a point where most TFBSs are bound by phosphorylated TFs, which, for most RRs, requires high phosphorylation levels of RR to occupy the TFBSs. The phosphatase activity of bifunctional HKs plays an essential role in setting the phosphorylation capacity. As noted previously [15,16], upregulation of the bifunctional HK via positive autoregulation also increases the phosphatase activity, offsetting the positive feedback and limiting the bistability range. Any mechanism that lowers the phosphatase activity can potentially promote bistability in positively autoregulated TCSs. Bimodal and heterogeneous responses have been observed for a PhoQ phosphatase mutant [51,52] and PhoB cross-activated by a high amount of non-cognate HKs [53], both of which lack phosphatase activity. Phosphatase activities can protect TCSs from bistability and our analyses suggest that the "protection" range with low logarithmic gain is roughly correlated with low phosphorylation. Many autoregulated TCSs display a maximum percentage of RR phosphorylation around or below 50% [45, [54][55][56]. It is likely that these autoregulated TCSs operate in a range with low RR phosphorylation to ensure monostability.
Lack of full RR phosphorylation may represent the cost for positively autoregulated TCSs to maintain a monostable and graded response. Similarly, other parameter values that facilitate bistability or monostability can also carry costs in different aspects of response. For example, a low fold change f can reduce the propensity for bistability, but it restricts the range of RR protein levels, which can lead to a weak response to signals as shown in S5 Fig. High K � for other TFBSs limits the sequestration-based bistability but may also impact the binding of TFs to these TFBSs, among which are the regulatory targets for transcriptional responses. Satisfying different signaling needs often requires a concerted tuning of multiple parameters or even evolution of novel regulatory schemes [30,31]. Furthermore, TCS architectures can be much more complex than the canonical one analyzed here. Bistability and ultrasensitivity can originate from alternative configurations, such as phosphorelays and split HKs [10,57].
Overall, we investigate how TFBS competition impacts the mono-/bistability parameter space and how positively autoregulated TCSs can overcome bistability-promoting effects to ensure accurate responses. Our model is limited by a deterministic approach. The equilibrium method for modeling TF binding employed here will be inadequate for many TFs that are expressed in low copies and are subject to fluctuations in a noisy system. Stochastic simulation, such as the one described previously for decoy competition [42], may be required to further explore the constraints imposed by competing TFBSs. Furthermore, TFBSs can be promiscuous and may be bound by other TFs or DNA-binding proteins, such as nucleoid associated proteins. Binding competition will be altered when a considerable fraction of the TFBSs are not accessible due to occupancy by other proteins. We also assumed a universal slow degradation/dilution rate for all proteins and protein-DNA complexes. Differences in degradation rates between the two could impact predictions [58]. Nevertheless, our analyses provide a quantitative framework to account for the TFBS competition effect for designing synthetic signaling circuits or predicting the steady-state output for naturally occurring signaling systems.

Model description
The model is based on a deterministic description of TCSs. Such an approach allows simple treatment of various reactions, especially the complex phosphorylation/dephosphorylation reactions of TCSs. The positively autoregulated signaling circuit is separated into three modules (Fig 1A). Reactions of different modules occur at different timescales with DNA binding being the fastest at a timescale of seconds, and transcriptional regulation the slowest, usually longer than 30 min and extending to hours. Phosphorylation and dephosphorylation of many TCSs occur at a timescale in min, between the other two modules [59,60]. Therefore, the three modules are considered as decoupled for simplicity. Further, growth dilution is ignored for the phosphorylation and DNA binding modules because their reaction timescales are much faster than typical bacterial growth rates. Details of decoupling approximation and mass action kinetics are described in S1 Text.

Autoregulation module
Transcription and translation of the RR TF are lumped together to model the production of RR protein molecules, which depends on the binding probability of the promoter. Accounting for both production and growth dilution (see S1 Text for details), the ordinary differential equation (ODE) for total RR level R T is: in which α is the lumped protein production rate constant and k dil is the growth dilution rate. Normalizing all concentration parameters by K auto and letting The logarithmic open loop gain of the autoregulation module can be derived: LG This is used to derive the LG values at different R f � values and combined with LGs from the other two modules to explore the mono-and bistability range. The maximum value of LG F can be solved analytically: In autoregulated TCSs, the RR and the HK are commonly transcribed from the same operon. Expression of HK protein molecules is modeled with a rate constant proportional to that of RR expression so that a constant ratio of RR to HK is maintained. Impact of RR/HK ratios on phosphorylation is discussed in the modeling of the phosphorylation module.

DNA binding module
In the presence of competing TFBSs, TFs can be sequestered by DNA binding, lowering the concentration of free unbound active TF, R f . Because many bacterial TFs bind DNA as dimers and each TFBS usually contains two half-sites, the binding cooperativity is set to 2 for the binding model discussed here. More analyses are provided in S1 Text for TFBSs with alternative binding cooperativity. An extreme case with all TFBSs having the same affinity K is modeled. With the DNA concentration at D, R P is the sum of R f and DNA-bound phosphorylated TF molecules: All concentrations and relevant parameters are normalized by K auto to ensure consistency of variables from different modules: Both sides of Eq (11) are differentiated with respect to R P � to obtain dR � f =dR � P . The LG of the DNA binding module can be calculated as: ð12Þ LG B values are computed with Eq (12) above and multiplied by LG F at individual R f � levels. The maxima of combined LG values within a R f � range are numerically computed to explore how the concentration and affinity of TFBSs affect the monostability range.

Activation/phosphorylation module
Phosphorylation of the RR is modeled as described previously [18,22,28]. Degradation or growth dilution rates of TCS proteins are considered to be very slow in comparison to reaction rates within the phosphorylation module, and thus are ignored (see details in S1 Text). Concentration of the phosphorylated RR TF can be derived as follows: in which Cp and Ct are composite parameters defined in Fig 5A and [RR] is the concentration of the free unphosphorylated RR. From Eq (13), we can conclude P P < Cp, and Cp defines the maximal RR phosphorylation level. Based on mass conservation, When the RR is in great excess to the HK, the last two terms for the HK-RR complexes can be ignored and [RR] ffi R T − R P . Substituting this into Eq (13) and normalizing all concentrations, we obtain: Solving Eq (15) gives: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi This is the RR phosphorylation approximation equation described previously [22]. Logarithmic gain of the phosphorylation module can be obtained from Eq (16) (see S1 Text for details): LG Because R � T � R � P ; Cp � � R � P , the last term of Eq (17) is non-negative, and we can conclude: LG A � 1: It can also be derived from Eq (17) that: LG A < 0:5: To calculate the overall LG, R P � and R f � are numerically computed for each R T � level and corresponding LG values for three modules are multiplied to derive LG 3 .

Parameter ranges
Parameter ranges are estimated based on experimental measurements from well-studied TCSs.

Bacterial strains and growth conditions
E. coli strain BW25113 is the parent strain for all strains constructed for TF activation assays. Similar to methods described previously [26], PcusC and PcusR promoters were amplified by PCR and cloned into pJZG146 that has a promoter-less yfp gene, resulting in pJZG157 and pJZG209 that were used for estimating fold change of transcription. To examine CusR activation in single cells, a PcusC promoter was fused with the mGreenLantern gene [38] to obtain stronger fluorescence than YFP. The fluorescent reporter was integrated into the HK022 phage attachment site in the chromosome of BW25113 using recombination strategies [64], yielding RU2118. Decoy CusR binding sites were made by annealing oligonucleotides (see details in S1 Text) with sequences based on previous analyses [34,65] and inserted into pRG475 [18], a plasmid with a copy number that has been determined previously and that can be further manipulated by arabinose. The resulting plasmids pCusRBS1 (1 site/plasmid), pCusRBS2 (2 sites), pCusRBS3 (3 sites) and pLH10 (0 site) [18] were introduced into RU2118 for flow cytometry analyses. All strains were grown at 37˚C in LB broth or minimal A medium [66] supplemented with 0.4% (w/v) glucose and amino acids (40 μg/ml each).

CusR activation assay
Bacterial cultures grown in minimal A medium overnight were used to inoculate fresh cultures with~1:20 dilution. Fresh cultures in mid-log phase were harvested and resuspended in fresh media containing different concentrations of CuSO 4 (0-500 μM). The starting optical density measured at 600 nm (OD 600 ) was 0.15 and bacterial cultures were transferred to either 96-well plates for fluorescence reporter assays, or flasks for flow cytometry. For reporter assays measured by a Varioskan plate reader (Thermo Scientific), cells were incubated with shaking with fluorescence and OD 600 measured every 5 min. Data were processed as described previously [18]. For flow cytometry analyses, cells were incubated for 2 h and analyzed using a Gallios Flow Cytometer (Beckman Coulter). GreenLantern fluorescence (FL1) and CFP (FL10) were measured for~60000 cells per sample. The binding site plasmids also carry a constitutively expressed cfp gene as a control. Dead cells with no CFP expression were excluded from reporter analyses. Arabinose can alter the plasmid copy number and change CFP fluorescence. Thus, the relative plasmid copy number was calculated by comparing median CFP fluorescence from cells with 0.2% arabinose to cells without arabinose. The number of CusR binding sites was derived by multiplying the plasmid copy number by the number of binding sites on individual plasmids.

Evaluation of TFBS binding strength
To assess binding affinities of TFBSs, individual information content (R i ) was calculated as described [25] for curated sets of TFBSs of PhoB from E. coli [26,67], PhoP from Salmonella [68], and BvgA from Bordetella [69]. Peak intensities from chromatin immunoprecipitation (ChIP) experiments can reflect binding strengths of a TF to various sites. ChIP peak intensities, defined as signal noise ratio (SNR), were obtained for autoregulated TFs from the prokaryotic ChIP database, proChIPdb [27]. After subtracting the median of logarithmic intensities of all identified peaks, the resulting values (Δlog 10 SNR) were used to assess the relative binding strength of individual TFBSs. Definition of positive or negative autoregulation is based on Reg-ulonDB [12]. Three other autoregulated TFs (ArgR, FNR and PuuR) in proChIPdb were not included due to the lack of peak intensity for the autoregulated promoter. Two other TFs with a fairly large number of TFBSs, YfeC and YdhB, were chosen because of the presence of binding sites within their own promoters [70] however, the regulatory role of these sites is unknown.  (16) and (17) with our approximation model discussed in the main text, assuming the RR is in great excess to the HK. Colored lines represent simulated data with the indicated RR/HK ratios using the full model described in S1 Text. Significant deviation from the approximation model is only apparent with an RR/HK ratio of 1. Parameter values are as follows: Cp � , 18.3; Ct � , 2.1; K auto , Transcriptional response is defined as the promoter occupancy that is determined by R f � and the binding affinity K � rep . K � rep will be different for different promoters, such as the autoregulated RR promoter and the RR-regulated promoters that contain sites among the competing TFBSs. (B-D) Impacts of R b � and K � on R f � , and subsequently the promoter occupancy. For systems with identical parameter sets (f = 5, D � = 4 and K � = 2), R b � determines the concentration range of R f � levels (blue or purple shaded areas in B) and the range of promoter occupancy in response to signals. In comparison to a weak TFBS affinity (K � = 2), a strong TFBS affinity (K � = 0.5) always leads to stronger TFBS competition, thus lower R f � and lower promoter occupancy (C and D). For high R b � and low K � rep values, such differences can be small, and both (purple and red lines in C) can reach near full occupancy. (E) Impacts of R b � and f on transcriptional responses. High f can result in bistable responses (pink) while low f gives monostable but extremely weak responses (navy). RU2118 (PcusC-mGreen-Lantern) carrying different decoy plasmids were assayed. All decoy plasmids carry a constitutively expressed CFP used for gating of the flow cytometry data (A) and estimation of plasmid copy number (B). A small population of cells with only background CFP fluorescence was often observed, especially in the presence of high concentrations of toxic CuSO 4 . None of these cells showed GreenLantern (GL) fluorescence. They are considered as dead cells or cells with impaired gene expression, thus were gated out with CFP histograms and excluded from further analyses. (B) Plasmid copy number estimation. The copy number of decoy plasmids has been determined as 96 in previous studies [18]. Addition of 0.2% arabinose (Ara) can inhibit replication of the plasmid origin and reduce the copy number, leading to reduced CFP fluorescence (left panel). The median of CFP fluorescence was used to derive the relative fluorescence. Solid and open circles represent CFP fluorescence of independent samples in the absence and presence of arabinose (middle panel). Arabinose reduced fluorescence to 62%, thus the copy number is estimated to be 0.62×96�60. Plasmid DNA extracted from corresponding cultures showed a similar value of the relative DNA amount (right panel), consistent with the relative copy number estimated from CFP. (C) Dot plots showing GL and CFP fluorescence of individual cells in response to 500 μM CuSO 4 . Decoy numbers are computed based on the plasmids and conditions used as follows: 0, pLH10 (0 sites); 60, pCusRBS1 (1 site, +Ara); 96, pCusRBS1; 120, pCusRBS2 (2 sites, +Ara); 192, pCusRBS2; 288, pCusRBS3. (TIF)