Phosphate Sink Containing Two-Component Signaling Systems as Tunable Threshold Devices

Synthetic biology aims to design de novo biological systems and reengineer existing ones. These efforts have mostly focused on transcriptional circuits, with reengineering of signaling circuits hampered by limited understanding of their systems dynamics and experimental challenges. Bacterial two-component signaling systems offer a rich diversity of sensory systems that are built around a core phosphotransfer reaction between histidine kinases and their output response regulator proteins, and thus are a good target for reengineering through synthetic biology. Here, we explore the signal-response relationship arising from a specific motif found in two-component signaling. In this motif, a single histidine kinase (HK) phosphotransfers reversibly to two separate output response regulator (RR) proteins. We show that, under the experimentally observed parameters from bacteria and yeast, this motif not only allows rapid signal termination, whereby one of the RRs acts as a phosphate sink towards the other RR (i.e. the output RR), but also implements a sigmoidal signal-response relationship. We identify two mathematical conditions on system parameters that are necessary for sigmoidal signal-response relationships and define key parameters that control threshold levels and sensitivity of the signal-response curve. We confirm these findings experimentally, by in vitro reconstitution of the one HK-two RR motif found in the Sinorhizobium meliloti chemotaxis pathway and measuring the resulting signal-response curve. We find that the level of sigmoidality in this system can be experimentally controlled by the presence of the sink RR, and also through an auxiliary protein that is shown to bind to the HK (yielding Hill coefficients of above 7). These findings show that the one HK-two RR motif allows bacteria and yeast to implement tunable switch-like signal processing and provides an ideal basis for developing threshold devices for synthetic biology applications.

1 Mathematical analysis of the simple model

Reactions
The model has 6 species: one histidine kinase HK, two response regulators R1, R2 and their corresponding phosphorylated forms: HKp, R1p, R2p. The first model we analyze consists of the following reactions: • Auto-phosphorylation of HK: HK ka − → HKp • Reversible phosphotransfer between HK and R1: • Reversible phosphotransfer between HK and R2: • Auto-dephosphorylation of the response regulators:

System of ODEs
To simplify the notation, we write the concentration of each of the species as: We let x = (x 1 , . . . , x 6 ) be the vector of concentrations. We use mass-action kinetics to model the dynamics of the concentrations of the species in time with a system of ordinary differential equations (ODEs). We write the ODE system as the product of the stoichiometric matrix A and the vector of reaction rates. We order the reactions as shown above and the species analogously to the concentrations order x 1 , . . . , x 6 . The stoichiometric matrix A and the vector of reaction rates v(x) are: We letẋ(t) denote dx(t) dt and drop the dependency on time t in the notation, that is, we writeẋ i and x i foṙ x i (t) and x i (t) respectively. The system of ODEs modeling the dynamics of the concentrations in time is ẋ = Av(x), that is:ẋ 1 = k M x 2 x 5 + k S x 2 x 3 − k rM x 1 x 6 − k rS x 1 x 4 − k a x 1 (S1) At steady state,ẋ i = 0 for all i = 1, . . . , 6.
The system has three conservation laws: x 1 +ẋ 2 = 0,ẋ 3 +ẋ 4 = 0,ẋ 5 +ẋ 6 = 0, which reflect that the sum of species concentration in unphosphorylated and phosphorylated form is constant at each layer. These relations provide three equations at steady state: for some positive total amounts H, R S , R M . Due to these relations, three of the ODEs, one for each conservation law, are redundant for the computation of the steady states. We choose to disregard (S1), (S3) and (S5). Therefore, the steady states of the system are the solutions to the steady state equations corresponding to (S2), (S4) and (S6) and the three conservation laws in (S7). That is, the steady states are the solutions to the following six equations: The equations are easier to analyze if we replace the first equation by the sum of the first three equations. This change does not alter the solutions to the system. The new system to be solved is:

Expressions for the steady states
We recursively write all concentrations at steady state as a function of x 6 . We further impose that all concentrations are positive at steady state. We solve for x 5 at steady state using (S13) and obtain The right-hand side of this equality decreases in x 6 . Further, x 5 , x 6 are positive provided that 0 < x 6 < R M . Let b 1 be this first upper bound for x 6 at a positive steady state, that is We next solve for x 1 , x 2 using (S11) and (S10) and obtain At steady state x 2 is positive provided x 5 , and x 6 are positive. For x 1 to be positive, we need k M Hx 5 − k hM x 6 > 0. Substituting x 5 by the expression in (S14), this is equivalent to require that k M H(R M − x 6 ) − k hM x 6 > 0. Isolating x 6 we obtain that the numerator of x 1 is positive if and only if and hence the upper bound b 1 of x 6 in (S15) is larger than b 2 in (S17). Therefore x 6 < b 1 is satisfied if also x 6 < b 2 and the upper bound for x 6 , b 1 , can be ignored. The expression for x 1 in (S16) decreases in x 6 and increases in x 5 . Since x 5 decreases in x 6 , we conclude that x 1 decreases in x 6 . Since x 2 = H − x 1 (see (S11)), we have that x 2 decreases in x 1 and hence x 2 increases in x 6 .
We solve for x 3 , x 4 using (S9) and (S12) and obtain At steady state both x 3 , x 4 are positive provided x 1 , x 2 are positive. The expression for x 3 decreases in x 2 and increases in x 1 . Since x 1 decreases in x 6 and x 2 increases in x 6 we conclude that x 3 decreases in x 6 . Since from (S12) we have that x 3 = R S − x 4 , we conclude that x 4 increases in x 6 . Summary: By iterative substitution, all concentrations at steady state are expressed as functions of x 6 . Further, all steady state concentrations are positive if and only if (S19)

Signal-response curve and maximal response
The value of x 6 at steady state is determined using the only equation not used so far, (S8), after writing all other concentrations as functions of x 6 . This equation provides the analytical description of the inverse of the signal-response curve, written as k a = f (x 6 ). Solving (S8) for k a we obtain Since x 1 decreases in x 6 and x 4 increases in x 6 , the function f is increasing in x 6 . For positive steady states, this function needs to be evaluated at the values of x 6 satisfying (S19). When x 6 = 0 then k a = 0. As we observed above, when x 6 approaches b 2 we have that x 1 approaches 0, x 2 approaches H, and hence using (S18) we have that x 4 approaches k S HR S k S H+k hS . Therefore it follows from the expression (S20) that k a tends to infinity when x 6 tends to the upper bound b 2 .
We conclude that for x 6 fulfilling (S19), f is an increasing function with range (0, +∞). It follows that given a value of signal, k a > 0, then there exists a unique value of x 6 fulfilling (S19) and k a = f (x 6 ). This is the steady state value of x 6 for the given input k a . The concentrations of the other species at steady state are uniquely determined using their expressions as functions of x 6 . As a consequence, we showed that multistationarity cannot occur.
The explicit expression of the function f (x 6 ) in (S20) is: The signal-response curve is the inverse of the function f . Since the function f (x 6 ) increases, so does the signal-response curve. Further, we conclude that the upper bound b 2 is the level of phosphorylated R2 that the system attains when k a increases infinitely. In other words, is the maximal level of response of the system. The maximal response levels for HK and R1 are H and respectively. Observe that the maximal response (S22) does not depend on any of the rate constants or total amounts involving the sink-RR (R1). Therefore, the presence of a sink does not alter the maximal response, but might alter the shape of the signal-response curve. More generally, the inverse of the signal-response curve without the sink is: This is obtained by showing as above that in this case Alternatively, we can set R S = 0 and all rate constants of the sink system to zero in the expression of f in (S21). Note that the expression of x 1 in (S16) does not depend on the presence of sink. Therefore, for any fixed x 6 , (S24) and (S20) satisfy the inequality after substituting x 1 , x 4 for their expressions in x 6 and for any choice of rate constants of the sink system and R S . We deduce that given a response x 6 , k a in (S20) is larger than k a in (S24). In other words, presence of the sink causes the system to require more signal to achieve the same level of response. If we plot the signal-response curve for a system with a sink, then the graph is below the signal-response curve of the corresponding system without sink, while keeping the common rates the same.

First derivative at zero
The first derivative of the signal-response curve at zero is computed as 1/f (0) and takes the value: The signal-response curve of the corresponding system without sink has first derivative at zero: Observe that (S25) can be rewritten as and it becomes apparent that the first derivative at zero of the signal-response curve when the system has a sink is always smaller than the corresponding first derivative at zero of the signal-response curve of the system without a sign. Therefore, presence of the sink reduces the slope of the curve around zero.

Sigmoidal vs. hyperbolic
A hyperbolic curve has negative second derivative in all its domain, while a sigmoidal curve displays a change of sign of the second derivative. In particular, the second derivative at zero must be positive for a sigmoidal curve and negative for a hyperbolic curve.
The second derivative at zero of the signal-response curve can be computed from its inverse, f . In particular, the sign of the second derivative of the signal-response curve at zero is minus the sign of f (0) and is given by the sign of: where For S to be positive, we necessarily need that k S k hS R S = 0 and k rS < k S as stated in the main text. These conditions are necessary for sigmoidality. In particular, presence of the sink is necessary. A simple system containing only one HK and one response regulator can only display hyperbolic signal-response curves. The above conditions are necessary for sigmoidality. Sufficient conditions can also be given. For example, if k rS k S < k rM k M + k rM then there exist total amounts H, R S large enough and R M small enough such that the signal-response curve is sigmoidal. This follows from a general fact on polynomials. Consider a polynomial p(x) = a n x n + · · · + a 1 x + a 0 . If a n > 0, then for x 0 large enough, p(x 0 ) > 0. Similarly, if a n < 0, then for x 0 large enough, p(x 0 ) < 0. We consider the term S in (S26) as a polynomial in H. This polynomial has degree 2 and the coefficient of the term of degree 2 is If k rS k M − k rM k S + k rM k rS < 0 or, equivalently, k rS k S < k rM k M +k rM , then for R S large and R M small the coefficient of the polynomial is positive and we can use the general fact on polynomials to conclude that for H large enough, S is positive and hence sigmoidality occurs.

Mathematical analysis of the model with intermediates
We consider here the case where the model includes complex formation in the phosphotransfer reactions.

Reactions
The model has now 8 species: one histidine kinase HK, two response regulators R1, R2, their corresponding phosphorylated forms, HKp, R1p, R2p, and two intermediates Y S , Y M . The model consists of the following reactions: • Auto-phosphorylation of HK: HK ka − → HKp • Reversible phosphotransfer between HK and R1 through the formation of a complex: • Reversible phosphotransfer between HK and R2: • Auto-dephosphorylation of the response regulators:

System of ODEs
To simplify the notation, we write the concentrations of each of the species as: We let x = (x 1 , . . . , x 8 ) be the vector of concentrations. We proceed as above to construct a system of ODEs for this model:ẋ At steady state,ẋ i = 0 for all i = 1, . . . , 8. The system has also three conservation laws: x 1 +ẋ 2 +ẋ 7 +ẋ 8 = 0,ẋ 3 +ẋ 4 +ẋ 7 = 0,ẋ 5 +ẋ 6 +ẋ 8 = 0, which provide three equations at steady state: for positive total amounts H, R S and R M . Three of the steady state equations are redundant and are substituted by the conservation law equations. We keep the steady state equations corresponding to (S28), (S30), (S32), (S33) and (S34). Therefore, the steady states constrained to the conservation laws are given as the solutions to the equations:

Sigmoidal vs. hyperbolic
We compute directly the second derivative of the signal-response curve at zero. First of all, we observe that when the signal is zero (k a = 0), the steady state is: We differentiate each equation in (S36) with respect to k a . To simplify the notation, we write and evaluate the concentrations x i at k a = 0. We obtain a new system of equations: The solutions in p i are the derivatives at zero of each of the concentrations at steady state as a function of k a . For instance, the first derivative of x 6 , that is, of the signal-response curve, at zero is We introduce new rate constants that allow us to compare the systems with and without intermediates. That is, we define the inverse of the Michaelis-Menten constants of the intermediates: k yS = k aS k arS + k bS , k yM = k aM k arM + k bM , k yrS = k brS k arS + k bS , k yrM = k brM k arM + k bM and let the new effective reaction rate constants for the phosphotransfer reactions be: With the new constants, the first derivative at zero of the signal-response curve is which is identical to the first derivative of the signal-response curve at zero without modeling of intermediates (given in (S25)). That is, with the identification of constants, the two derivatives are identical.
We compute now the second derivative of each variable at zero (k a = 0) following the same procedure as above. We differentiate the equations (S37), evaluate the concentrations x i at k a = 0 and the first derivatives p i at the solutions obtained in the first step. We obtain a new system of equations on the second derivatives of x i with respect to k a at zero and solve it. In particular, we obtain the second derivative of x 6 with respect to k a at zero and conclude that the sign of the second derivative of the signal-response curve at zero agrees with the sign of the following term: where ω 1 = k rM H + k hM ω 2 = k rS H + k hS and S corresponds to the system without intermediates, and was given in (S26) to be Observe that the terms not involving the Michaelis-Menten constants of the intermediates correspond precisely to the sigmoidality term S without intermediates. For S y to be positive one needs either S in (S26) to be positive or the blue term to be negative. The latter involves the Michaelis-Menten constants of the intermediates and, in particular, the condition k rS < k S is not necessary anymore. What is still necessary is that Let us look closely at the extra term that can be negative. Sigmoidality cannot occur if k rS > k S and The left-hand side of the inequality is written as a polynomial in H as All the coefficients of this polynomial are positive if Hence sigmoidality cannot occur if k rS > k S and (S39) holds. We conclude that necessary conditions for sigmoidality in the model with intermediates are that (i) k S k hS R S = 0, and (ii) k S k rS > min 1, In particular, if k yM < k yS , then we are left with the same necessary conditions as in the case without intermediates. In terms of the original reaction rate constants, we note that Then necessary conditions for sigmoidality are that Finally, we look for conditions that guarantee the existence of sigmoidality for some total amounts. To this end we proceed as above and consider S y as a polynomial in H. The polynomial has degree five and the coefficient of highest degree is: If the coefficient in is positive, then for H large enough, S y is positive and as a consequence the signalresponse curve is sigmoidal. The coefficient is positive if and only if or equivalently (cf. (S40)), if and only if k aM k arM < k aS k arS .

A mathematical model for the two-component system regulating yeast osmoregulation
To model the one HK-two RR system found in yeast osmoregulation, we considered its dynamics in isolation of other cellular components. In this system, the HK (SLN1) is a hybrid protein composed of histidine kinase (HK) and receiver (RC) domains and there exists an additional histidine phosphotransfer (HPT) protein (YPD1). Together, these constitute a phosphorelay, where the YPD1 phosphotranfers to the two RRs, SSK1 and SKN7. We have modeled the hybrid HK as two separate proteins. The reactions in this system are: where HK, RC, HPT, RR1 and RR2 stand for SLN1, its receiver domain, YPD1, SSK1, and SKN7 respectively and the -p suffix represents phosphorylated forms of these proteins/domains. The above reaction scheme can be used to derive a system of ordinary differential equations, which describe the changes in concentrations over time. For the phosphorylated forms, the system is: To analyze the behavior of the phosphate sink motif with increasing signal, we simulated the incoming signals from receptors as an increase in the auto-phosphorylation constant rate of the kinase (k a ). The model was parameterized with data from literature (see Supplementary Table 1). We numerically integrated the model to derive steady state signal-response relationships. The latter analysis gives the steady state level of phosphorylated RR2 at a given signal (k a ), where signal was taken as the constant rate of autophosphorylation of kinase and allows deriving a so-called signal-response curve. This curve is found by numerically integrating the system to steady state at a fixed signal level and then numerically "following" this steady state (i.e. steady state level of phosphorylated RR2), while changing the signal. This analysis is equivalent to allowing the system to reach steady state under different signal values. Both signal-response analyses were performed using the software Oscill8 (http://oscill8.sourceforge.net/).

Analytical comparison of different models
To perform a formal check for the potential of bistability in the different models (discussed in the main text), we have utilized the chemical network theory. This theory provides several analytical tests that can provide a definite answer on the possibility of existence of multiple stationary states in a given reaction network. We have applied these tests to the basic model and model with alternative reaction schemes; we had devised using the Chemical Network Tool v2.2 (http://www.chbmeng.ohio-state.edu/ feinberg/crntwin/). The model files used with this tool and describing the chemical reaction systems, as well as the analytical results from the tool are provided as separate files. Supplementary Text 2 contains the reaction system described in the basic model without complex formation. Supplementary Text 3 contains the reaction system with formation of complexes during phosphotransfer reactions. The four complexes between the phosphorylated/unphosphorylated CheA and the phosphorylated/unphosphorylated CheY1 and CheY2. In these reaction systems A, Y1 and Y2 stand for CheA, CheY1 and CheY2 respectively. P refers to phosphorylated form.