In order to understand the physiological effect of a third component (TC) on the function of a prototypical TCS, we built models of TCS with and without that TC and compared the dynamical behavior of those models. shows a schematic representation of the three models used in our analysis. These models are mathematically described by using a mass action system of ordinary differential equations (ODE) 
. The resulting ODE systems for each of the three alternative models can be analyzed and compared numerically by running appropriate simulations on a computer.
Analyzed Two Component Systems modules.
Models and Comparisons
The network model that we use to describe the prototypical TCS in our analysis is that defined in Igoshin et al.
, which is based on earlier work 
. In Model A, shown in , the SK can autophosphorylate and/or autodephosphorylate in response to an external signal. Both phosphorylated and unphosphorylated forms of SK are allowed to bind RR with similar affinities, as reported in 
. Binding of unphosphorylated SK and RR is reversible and forms a dead-end complex (SKRR). Phosphorylated SK (SKP) can transfer its phosphate to the RR. The phosphorylated RR (RRP) will modulate the biological levels and activity of relevant proteins.
This network for the prototypical TCS was modified to study the effect of a TC binding to either the SK or the RR. The changes in the network are also shown in . Model B represents a TCS where a third component binds to the SK (TCSK
), inactivating it. Model C represents a TCS where a third component binds to the phosphorylated RR (TCRR
) and stabilizes this phosphorylated form. In prototypical TCS modules with bifunctional sensors
, the unphosphorylated SK can destabilize the phosphorylated form of the RR and it increases the dephosphorylation rate of RRP (k8
>0 in ). In prototypical TCS modules with monofunctional sensors
, the unphosphorylated SK has no effect upon the dephosphorylation rate of RRP (k8
0 in ). The model includes a phosphatase that dephosphorylates RRP independently of the SK. This is done for generality. In the cases where no such phosphatase exists, this set of reactions can be replaced by a single reaction where the unstable RRP phosphate bond hydrolyzes over time. An appropriate choice of parameter values will make the results of the analysis similar to those described for the full model.
In this study we analyze the potential effect of a TC in the physiological behavior of TCS modules with bifunctional and monofunctional sensors independently. If the TC has no effect on the physiological behavior of the TCS, then the presence of TC in particular instances of TCS should be understood as an evolutionary accident. If the TC has an effect on the physiological behavior of the TCS, this could provide a rationale for the selection of a TCS design that includes a TC. To perform the analysis, we compare the dynamical behavior of Model A to that of Models B or C, independently. This comparison is done in two ways.
First, Models A and B (or C) are compared ensuring that the parameter values of all processes that are common are the same in the two models. This guarantees that whatever differences are found are only due to the addition of the TC. This comparison is equivalent to comparing an organism where the TCS interacts with a TC to another where the TC has been deleted from the genome. This situation could occur, for example during the creation of a new biological circuit by genetic manipulation in a biotechnological context. Thus, this type of comparison is relevant for understanding the differences in behavior of biological circuits created using synthetic biology techniques.
Second, we also perform a mathematically controlled comparison between Models A and B (or C). This is a well established method for evaluating the irreducible effect of a change in the design of a biological circuit on the dynamic behavior of the network 
. In this comparison, in addition to ensuring that Models A and B (or C) have the same values for corresponding parameters of all processes that are common, we use the differences between the designs as degrees of freedom that evolution can use as a substrate to minimize differences between the dynamic behavior of the two systems. If the alternative designs can be made equivalent by using these degrees of freedom, then one may argue that they cannot be distinguished by natural selection. If, after making the systems as equivalent as possible, there are still irreducible differences in the physiological behavior between designs, then one may expect one of them to be preferably selected when its functionality provides better adaptive advantage. In the models under comparison, the difference is the deletion of a protein from the module between Model B (or C) and Model A. In this situation, the protein burden caused by Model A is lower than that caused by its alternative designs. Hence, we allow that the system changes the total concentrations of the remaining proteins (SK and/or RR). The details for this comparison are given in the methods
section. This comparison is thus relevant for understanding the differences in the dynamic behavior that are intrinsic to the differences in design between Models A and B (or C), and to those alone, in evolutionary terms.
Effect of a third component on TCS signal amplification and bistability
Signal amplification is an important physiological property of TCS. TCS with appropriate signal amplification can provide evolutionary advantages to organisms harboring them. Thus, understanding how signal amplification is affected by adding a TC to a TCS would help in predicting under which conditions to expect such a design to be selected. shows that all models can achieve the same signal amplification, whether the environmental signal modulates the autophosphorylation (k1) or the autodephosphorylation (k2) of the SK. This can be seen because the difference between the amount of RRP (phosphorylated RR) when k1 is low (k2 is high) and when k1 is high (k2 is low) can be similar for all models. Nevertheless, Model B responds at higher signal intensities and Model C responds at lower signal intensities than Model A, when the stimulus modulates the SK autophosphorylation reaction rate (compare the curves for k1 response of Model A to those of Models B and C in ). When the signal modulates the SK autodephosphorylation reaction rate, Model B responds at lower signal intensity and Model C at higher signal intensity than Model A (compare the curves for k2 response of Model A to those of Models B and C in ). However, mostly, the differences in signal intensity at which the systems are turned ON or OFF are small.
Steady state signal-response curves for the various TCS modules.
In addition, the prototypical TCS shown in Model A can show bistable behavior 
, making it possible that a signal can lead to one of two alternative responses, depending on the history of the system. Such a response may have some evolutionary advantages, for example in situations like sporulation where an irreversible developmental decision is made by cells. Bistable regions in the curves of have three values of RRP for a single value of signal intensity. The two extreme values are the alternative stable steady states, while the middle value is a biologically irrelevant unstable steady state that is mathematically required to exist if two stable steady states are present. In the figure one can see that the signaling ranges where bistability exists are different if the environmental signal modulates the autophosphorylation (k1
) or the autodephosphorylation (k2
) of the SK.
Necessary, although not sufficient, conditions for the existence of such bistable behavior in the prototypical TCS are i) the formation of a dead-end complex between the dephosphorylated forms of SK and RR and ii) that a sufficiently high fraction of the flux for the dephosphorylation of RRP is independent of SK. To understand how the presence of a TC affects the possibility of a bistable response in the prototypical TCS, we analyzed Models B and C in search of the existence of multiple steady states, followed by a comparison of the physiological behavior between Models A and B, and between Models A and C.
Given that signals can in principle modulate either the autophosphorylation (k1) or the autodephosphorylation (k2) rate of SK, we performed parallel computational experiments independently modulating their intensity. These experiments were done independently for models with monofunctional and bifunctional SK ().
Our results show that, in an uncontrolled comparison, the range of bistability for the bifunctional prototypical TCS is larger than if a TC binds any of the proteins of the module (compare panel B to panels D and F of ). Bistability for Model B in panel D is only observed for k1 signaling, while no bistability is observed for Model C in panel F. On the other hand, the range of bistability for the monofunctional prototypical TCS is larger than if a TC binds the RR of the module (compare panel A to panel E of ), but smaller than if the TC binds the SK (compare panel A to panel C of ). Differences among the three systems are more pronounced when the signal induces dephosphorylation of the SK (k2), rather than inducing SK autophosphorylation (k1).
An additional definition is needed before further presenting and discussing the results. Hereafter the system is said to be in an ON state if most of its RR is in the phosphorylated RRP form. If most of the RR is in its dephosphorylated form, the system is said to be in its OFF state. With this in mind, and as one might expect, systems with a TCSK are in an ON state for a smaller signaling range (panels C and D) and systems with a TCRR are in an ON state for a larger signaling range (panels E and F), in comparison with the uncontrolled Model A (panels A and B).
When the comparisons are controlled we see that the response of Model A can become similar to that of Model B or C by adjusting the total amount of available SK. If the response of Model B is to be mimicked, the total amount of SK in Model A is decreased (, panels C and D, see methods
for the exact values of the total amount of SK), while mimicking the response of Model C leads to an increase in the concentration of SK (, panels E and F, see methods
for the exact values of the total amount of SK).
The k2-response curves in panels B and C show that the switch from ON to OFF (from high to low levels of RRP) in these models could be irreversible or very difficult to reverse. In other words, modulation of the autodephosphorylation rate of SK by an external signal could generate nearly irreversible biological switches.
Our simulations also show that the necessary conditions for bistability in prototypical TCS remain necessary in the TCS with a TC. If either no independent phosphatase is present in the system (Ph
0) or no dead end complex is formed (k10
0) all TCS modules analyzed here are monostable (see section “Effect of changes in SK-independent RRP dephosphorylation and SKRR affinity on bistability” below).
In summary, a TCRR causes a reduction in the TCS parameter space of bistability and an increase in the signaling range in which the system is in the ON state (responds at lower k1-signal intensity and at higher k2-signal intensity), whether the SK is monofunctional or bifunctional. This can be more effectively compensated by prototypical TCS through a change (an increase) in the concentration of the SK. In contrast, TCSK increases the signaling range in which the TCS can show a bistable response if and only if the SK is monofunctional and the environment modulates k2 (SK dephosphorylation rate). The behavior of TCS with a TCSK can be mimicked by prototypical TCS through a change (a decrease) in the concentration of the SK.
Effect of a third component on TCS response time
In addition to signal amplification, the response time to signals is an important physiological property of TCS. In evolutionary terms, a change in response time may have important consequences to the fitness of the system. Therefore, we analyzed the effect of a TC on the response times of the TCS. To do this we performed four independent sets of experiments for each of the models, and independently considering systems with a monofunctional SK and with a bifunctional SK. In experiments 1 and 2 we instantaneously change the signal k1 and measure how long the system takes to come within 90% of its new steady state. This measures the response time of the system if the physiological signal modulates SK phosphorylation. In experiments 3 and 4, we instantaneously change the signal k2 and measure how long the system takes to come within 90% of its new steady state. This measures the response time of the system if the physiological signal modulates SK dephosphorylation. The details about how the experiments were run are as follows:
- 1 - We set each system to its OFF state, with k1=10−5 s−1. Then, we increased the value of k1 to a value k1 higher and measured how long the system took to get to within 90% of its new steady state value. k1 higher was systematically changed between 10−5 and 10 s−1.
- 2 - We set each system to its ON state, with k1=10 s−1. Then, we decreased the value of k1 to a value k1 lower and measured how long the system took to get to within 90% of its new steady state value. k1 lower was systematically changed between 10−5 and 10 s−1.
- 3 - We set each system to its OFF state, with k2=10 s−1. Then, we decreased the value of k2 to a value k2 lower and measured how long the system took to get to within 90% of its new steady state value. k2 lower was systematically changed between 10−5 and 10 s−1.
- 4 - We set each system to its ON state, with k2=10−5 s−1. Then, we increased the value of k2 to a value k2 higher and measured how long the system took to get to within 90% of its new steady state value. k2 higher was systematically changed between 10−5 and 10 s−1.
Results are shown in . We see that the response times increase by more than two orders of magnitude when the new parameter value k·lower
approaches the threshold value for exiting the bistability region of a system. The peaks of slower response in the curves in are in the region of signal intensity that lies immediately beyond the border of the bistability ranges shown in . Given that the peaks of slower response are located at the exit of the bistable region, there is no peak in the signal-response time curve when the response is monostable or when there is an irreversible turning OFF of the system. Model B and Model A|B (A controlled for B) don't have a peak in their OFF to ON k2
-response times (Panel C of ) because these models irreversibly turn OFF after an increase in k2
(as depicted in panel C). Model C also has no peak in the response time (Panels C and D of ) because this model has a monostable response to changes in k2
(see panel E). In panels G and H of , neither of the three systems shows a peak in their signal-response time curve because of the lack of bistability in their signal-response steady state curve (see panels D and F). When Model A is compared to Model B in an uncontrolled manner, the time response peaks of Model A appear at signal intensities that are always lower than those where the peak appears in the response of Model B. When Model A is compared to Model C in an uncontrolled manner, the time response peaks of Model A appear at signal intensities that are always higher than those where the peak appears in the response of Model C (see Figure S1
Temporal responsiveness curves of Models A, B, and C.
In order to have a proxy of the integral temporal responsiveness of each system, we calculated the area under each of the signal- response time curves shown in . This area is the sum of all the transient response times for each signaling response. The values of these areas are given in and show that overall response times are similar between Models A and B. In contrast, Model A has a faster response than Model C. When the comparison is not controlled, differences between integrated response times of the three models are small, when the signal modulates autophosphorylation of SK. However, if SK dephosphorylation is modulated, Model B has the fastest integrated response, followed by Model A. Model C is, again, the slowest responder (Table S1
Controlled comparison of the overall response times between Models A and B, and between Models A and Ca.
In summary, Model B is a faster overall responder than the prototypical TCS when the system is turned ON by modulating the phosphorylation rate of the SK, and it is a slower responder in any other case. In contrast, Model C is always slower to turn ON or turn OFF than the prototypical TCS, under controlled comparison conditions.
Stochastic effects of a third component
Fluctuations in the amount of proteins that participate in biological reactions can lead to stochastic effects in the system's behavior, when the total number of proteins participating in reactions is small. We performed stochastic simulations to understand the role of stochasticity on the effect of the TC on the physiological response of the TCS networks. These simulations take into account that the number of TCS proteins present in the cell are typically in the 10–1000 molecules range.
The simulation experiments performed were similar to those described in experiments 1–4 of the previous section, although with a smaller number of data points. and show the results of these simulations.
Stochastic time trajectories after an instantaneous change in the signal, for the three systems modeled with a monofunctional SK.
Stochastic time trajectories after an instantaneous change in the signal, for the three systems modeled with a bifunctional SK.
The OFF→ON plots start with the system at the OFF steady-state (low concentration of active RR) corresponding to a low value of k1
) or a high value of k2
), and depict the temporal trajectory of the RRP concentration after an instantaneous increase in k1
or decrease in k2
, for three different values of k1
The ON→OFF plots start with the system at the ON steady-state (high concentration of active response regulator) corresponding to a high value of the signal k1
) or a low value of k2
), and depict the temporal trajectory of the RRP concentration after an instantaneous decrease in k1
or increase in k2
, for three different values of k1
The simulation results for three different signal intensities are plotted in and . Three independent simulations are shown for each signal intensity. The values of k1 and k2 in each trajectory are chosen to be below, next to and above the threshold value at which the system switches from OFF to ON, or from ON to OFF (in the cases in which this threshold exists). Because each system has a different threshold value, the parameter scan is different for each plot.
The results from the analysis of the continuous model are consistent with the stochastic simulations: as discussed in the previous section (), in systems with a signal range of bistability the response times increase when the signal intensity is near the threshold value at which the system exits the bistability region. One can see in and that, in many cases, the curves that correspond to a signal that is just outside of the bistability range do not reach steady state during the simulation time. These curves correspond with the peaks in .
Furthermore, our simulations predict that the systemic response becomes noisier as the signal intensity approaches the threshold value for bistability. Just above and just below this value there is an increase in the stochastic fluctuations of the system. This can be seen because the triplicate curves corresponding to these values in and are much more different among themselves than the triplicate curves for the signals away from this threshold.
The response in the systems A, B and C is noisier when k1 is modulated than when k2 is modulated. The OFF to ON trajectories of Model B after an instantaneous decrease in k2 confirm that the turn OFF of this system due to an increase in k2 is irreversible and the system can't return to the ON state (see panel C). The system C does not have a bistability region in its k2-response curve (see panels E and F). Therefore, we don't find a range of k2 values for which the systemic response becomes slower and noisier.
Robustness of the analysis
The analysis thus far was done using the specific set of parameter values reported in . In order to study the generality of the results we performed sensitivity analyses of the bistability to changes in the different parameter values and concentrations of the systems. The results of the controlled and uncontrolled comparison between Model A and Model B or C with respect to the effect of changing parameter values on a possible bistable response of the TCS are summarized in . The detailed results are shown in Figure S2
, where we show a set of two-dimensional sections of the multidimensional parameter space in which bistability is observed.
Basal values for the parameters and concentrations of the models in .
Percentage of parameter space where bistable responses are possiblea.
Overall, a system with a TCSK appears to have a wider parameter range of bistability if the SK is monofunctional, and a lower parameter range of bistability if the SK is bifunctional, while a system with a TCRR appears to have a lower parameter range of bistability, for systems with either a monofunctional or a bifunctional SK, when either system is compared to a prototypical TCS. However, if the comparison between Model A and Model B or C is controlled, then we see that the robustness of the parameter range of bistability is larger in the prototypical TCS (Model A) with only one exception: in systems with a bifunctional SK, Model C has a more robust parameter range of bistability.
Effect of changes in SK-independent RRP dephosphorylation and SKRR affinity on bistability
SK-independent RRP dephosphorylation and SKRR complex formation are needed for bistable responses to exist in Models A, B, and C. In order to investigate how quantitatively changing these features affects bistability we performed the following computational experiments (). We independently and simultaneously changed the values for k8
(the reaction that regulates dephosphorylation by the SK) and k9
(changing the rate of dissociation between SK and RR) between 10−6
and 10. Then, we calculated the steady state(s) for each system at different values of the signal represented by the parameters k1
were independently and systematically scanned between 10−6
and 10 in logarithmic space at intervals of 0.01 units. The results are shown in and Figure S3
. shows that, overall, bistability is possible in Model C in a smaller interval of parameter values than that for Models A and B. However, the picture changes when we analyze only the parameters that directly influence the necessary conditions for bistability (k8
). For these parameters, Model C is the system where overall bistability is possible in a wider range of parameter values, followed by Model B. Model A is the one where bistability is limited to a smaller region of parameter values. Nevertheless, when Model A is controlled to have signal-response curves that are as similar as possible to those of either Model B or Model C, Model A becomes the system where bistable responses can occur in a larger fraction of the space for k8
, and k10
. For values of k8
below a threshold that depends on the system and is lower in Model B than in Model A, bistability is present in both models. Within the range of k8
values that permit bistability, an increase in k8
causes an increase in the k2
range of bistability (up to approximately six orders of magnitude for k2
at the threshold value for k8
). This is so, despite the enlargement of the fraction of RRP dephosphorylated by SK, because the increase in k8
causes an increase in the concentration of the SKRR dead-end complex (see Figure S4
). As k8
decreases, the range of signal k2
in which the models show bistability decreases steadily for a few orders of magnitude. Then, a lower boundary is reached and bistability is observed for one or less than one order of magnitude of k2
signal, independently of the value for k8
Experiments to analyze the effect of changes in different parameter values and protein concentrations on the range of bistability for the alternative TCS modulesa.
Percentage of parameter space where a bistable response is possible for Models A, B, and Ca.
Given that the formation of a dead-end complex between SK and RR is a necessary condition for bistability, we also want to understand the isolated effect of different fractions of RR and SK being sequestered into this complex on bistability. To understand the effect of changing the amount of SKRR dead-end complex on the signaling range in which the systems can be bistable we performed the following numerical experiment. First, we took each model from . Then, we systematically scanned the values of the parameters k9
, independently and simultaneously, between 10−6
and 10 in logarithmic space at intervals of 0.01 units. These parameters regulate the amount of SKRR that is formed. Finally, for each pair of values for k9
, we independently calculated the steady state(s) for each system at different values of the signal represented by k1
. Each of these parameters was independently and systematically scanned between 10−6
and 10 in logarithmic space at intervals of 0.01 units. The results are shown in and Figure S3
Bistability can be found only for intermediate steady state concentrations of SKRR. If too little or too much SKRR is formed, then no bistable response is possible. Overall, for bifunctional TCS, Model C has the largest range of SKRR steady state concentrations for which bistability is possible, followed by Model B. In its uncontrolled form Model A has the smallest interval of SKRR steady state concentrations where bistability is permitted. This interval of concentrations decreases further when Model A is controlled to be comparable to Model B. However, when Model A is controlled to be comparable to Model C, the range of SKRR steady state concentrations that enable bistability becomes the largest of the three systems. In monofunctional TCS, Model C has a smaller range of SKRR steady state concentrations for which bistability is possible than Model B.
The notion that Model C is the one in which bistable responses are less sensitive to changes in the steady state concentrations of SKRR (in consequence of changing the affinity between SK and RR) is misleading. Bistability is only found in this model if the affinity between the dephosphorylated forms of SK and RR is much larger than that between SKP and RR or SK and RRP. Given that the affinity between all forms of SK and RR was measured as similar, it is not likely that bistability can be found in vivo in systems that are represented by this model.
A similar experiment was made by changing independently and simultaneously the total amount of SK and RR, followed by independent calculation of the steady state(s) for each system at different values of the signal represented by k1
. Again, each of the parameters was independently and systematically scanned between 10−6
and 10 in logarithmic space at intervals of 0.01 units. The results are shown in and Figure S3
. They are consistent with the situation described for changes in k9
Effect of the SK/TCSK and RR/TCRR concentration ratios on bistability
In order to understand how the relationship between the total amounts of SK (RR) and TCSK
) influences the signaling range in which bistable responses are possible, we have performed a number of computational experiments. First, we took Models B and C from . Then, we systematically, simultaneously and independently scanned the total amounts of SK (RR) and TCSK
) in Model B (Model C), as described in . Finally, for each total amount of SK (RR) and TCSK
), we calculated the steady state(s) for each system at different values of the signal represented by k2
. This parameter was also systematically scanned between 10−6
and 10 in logarithmic space at intervals of 0.01 units. The results are shown in Figure S3
. We also performed similar test replacing k2
The range of signal k2 for which Model B can show a bistable response is observed to be dependent on the TC. Bistability is observed only within a narrow band of the SK-TCSK concentration space. Outside of this band, a bistable response cannot be observed. The range of total amount of SK in the system that may lead to a bistable response remains approximately constant for low total amounts of TCSK. However, within the band of total SK and TCSK in which bistability is observed, as total TCSK increases, the range of total SK amount that can generate bistable responses also increases. At concentrations of TCSK between approximately 2 and 7 µM, we find bistability for total SK concentrations between 0.2 and 0.001 µM or lower. At higher total TCSK concentrations, only small amounts of SK are available in free form. This prevents formation of the SKRR dead-end complex that is required for bistability.
As is the case in Model B, bistability in Model C can be achieved in a narrow band of the concentration space. However, within the range of values of this simulation, whatever the concentration of TCRR, the system can always show bistability.