Unlimited multistability and Boolean logic in microbial signalling

The ability to map environmental signals onto distinct internal physiological states or programmes is critical for single-celled microbes. A crucial systems dynamics feature underpinning such ability is multistability. While unlimited multistability is known to arise from multi-site phosphorylation seen in the signalling networks of eukaryotic cells, a similarly universal mechanism has not been identified in microbial signalling systems. These systems are generally known as two-component systems comprising histidine kinase (HK) receptors and response regulator proteins engaging in phosphotransfer reactions. We develop a mathematical framework for analysing microbial systems with multi-domain HK receptors known as hybrid and unorthodox HKs. We show that these systems embed a simple core network that exhibits multistability, thereby unveiling a novel biochemical mechanism for multistability. We further prove that sharing of downstream components allows a system with n multi-domain hybrid HKs to attain 3n steady states. We find that such systems, when sensing distinct signals, can readily implement Boolean logic functions on these signals. Using two experimentally studied examples of two-component systems implementing hybrid HKs, we show that bistability and implementation of logic functions are possible under biologically feasible reaction rates. Furthermore, we show that all sequenced microbial genomes contain significant numbers of hybrid and unorthodox HKs, and some genomes have a larger fraction of these proteins compared with regular HKs. Microbial cells are thus theoretically unbounded in mapping distinct environmental signals onto distinct physiological states and perform complex computations on them. These findings facilitate the understanding of natural two-component systems and allow their engineering through synthetic biology.

the unphosphorylated His-containing phosphotransfer protein, and Hpt P the phosphorylated form.
For notational convenience, we name the rate constants k 1 , . . . , k 6 , which differs from the notation in the main text. We denote the concentration of the species as follows: Under the law of mass-action, we model the dynamics of the concentrations over time by the following system of ordinary differential equations: where we writeẋ * = dx * dt and omit reference to time t, that is, we write x * = x * (t). Observe thatẋ 1 +ẋ 2 +ẋ 3 +ẋ 4 = 0 andẋ 5 +ẋ 6 = 0. It follows that the sums of concentrations x 1 + x 2 + x 3 + x 4 and x 5 + x 6 are constant over time. This leads to the extra equations for some positive total amounts H, T > 0.
The steady state polynomial. All concentrations are expressed as functions of x 5 , and we have not used the equationẋ 5 = 0. We can replace the equationẋ 5 = 0 by any linear combination of the steady state equations that involves this one. By doing so, the solutions to the equations do not change. We replace it by the equationẋ 5 +ẋ 1 −ẋ 4 = 0. This cancels out the quadratic terms in the equationẋ 5 = 0, and we obtain the equation Substituting into (7) the values of x 1 and x 3 in (3), (5), and further letting x 6 = T − x 5 , we obtain that, at steady state, it holds By clearing denominators, the positive solutions to (8) agree with the positive solutions to the polynomial The polynomial p(x 5 ) has degree 3. Any zero of the polynomial between 0 and T corresponds to a positive steady state. From (8), if x 5 ≥ T were a positive zero of the polynomial, then we would have 0 equal to a negative number, which is a contradiction. Therefore, any positive solution to the polynomial equation must fulfil that x 5 < T and hence provide a positive steady state. The polynomial p(x 5 ) has at most 3 positive roots. We show below in subsection 2.3 that there exist choices of rate constants and total amounts such that p(x 5 ) indeed has 3 positive roots. Therefore, there exist choices of rate constants and total amounts such that the system has 3 positive steady states.

Necessary conditions for bistability
Following Descartes' rule of signs, a necessary condition for p(x 5 ) to have 3 positive roots is that the coefficients of the polynomial have alternating signs. Since the leading coefficient is positive and the independent term is negative, a necessary condition is that the coefficient of degree 2 is negative and the coefficient of degree 1 is positive, that is: To derive the necessary condition for bistability stated in the main text, namely k 1 < k 3 (in the main text, k s1 < k s2 ), we consider equation (8) again. We rewrite it as If the right-hand side of the equation, call it ϕ(x 5 ), is an increasing function of x 5 for positive x 5 , then for any value of T there will be a unique value of x 5 such that (10) is fulfilled, and hence a unique positive steady state. Since (10) is derived from (7), the function ϕ(x 5 ) equals (3) and (5). Clearly, k 6 x 5 is increasing in x 5 . The derivative of k 1 x 1 + k 3 x 3 with respect to x 5 is The derivative is not necessarily positive for x 5 > 0. However, if k 1 > k 3 then the derivative of ϕ(x 5 ) is positive, implying that ϕ is an increasing function, and, as a consequence, bistability cannot arise for any value of T . To summarise, k 1 > k 3 implies that there is no bistability, and therefore, a necessary condition for bistability is that k 3 > k 1 .

Necessary and sufficient conditions for bistability
The conditions given above are only necessary for bistability, but their fulfilment does not guarantee bistability. We provide here necessary and sufficient conditions on all the parameters of the system for bistability. The parameters include the reaction rate constants and the total amounts. To obtain them, we apply Sturm's Theorem: Theorem 1 (Sturm). Let p(x) be a real polynomial. Define recursively the Sturm sequence by We are interested in the positive roots of the polynomial p(x) = p(x 5 ) in (9). In this case, a = 0, and we need to take b large enough, which is equivalent to considering instead of the sequence p 0 (b), . . . , p m (b), the leading coefficients of the polynomials p 0 , . . . , p m . This sequence is written as p 0 (+∞), . . . , p m (+∞). Observe that a = 0 is not a root of p(x).
2 The core model for n hybrid HKs competing for the same Hpt

Model description
We study here the core system consisting of n hybrid HKs competing for the same Hpt. We call such a system an nHK-Hpt system for short. In this case, there are n hybrid HKs, which we denote by HK i , for i = 1, . . . , n, and we use subindices 00, P0, 0P, PP to denote the phosphorylation state of each of them.
The set of reactions given in the previous subsection are reproduced for the n hybrid HKs. That is, for i = 1, . . . , n, the reactions for the transfer of phosphate group are as follows: and there is further the dephosphorylation reaction We denote the concentration of the species as follows: x 6 := [Hpt P ], for i = 1, . . . , n. Under the law of mass-action, we model the dynamics of the concentrations over time by the following system of ordinary differential equations: for i = 1, . . . , n. The system has n + 1 conservation laws. Namely, for i = 1, . . . , n we have for some H i > 0, and for T > 0,
Using the equationsẋ i,1 = 0,ẋ i,2 = 0,ẋ i,3 = 0 and (12), namely we obtain the algebraic expressions for x i,1 , x i,2 , x i,3 , x i,4 at steady state, depending on the value of x 5 at steady state, analogous to the expressions (3)-(6): These expressions are positive provided x 5 is positive. From (13) we have that We replace the steady state equationẋ 5 = 0 byẋ 5 + n i=1 (ẋ i,1 −ẋ i,4 ) = 0, and we obtain the equation Substituting into (19) the values of x i,1 and x i,3 in (15), (17), and further letting x 6 = T − x 5 , we obtain that, at steady state, it holds By clearing denominators, that is, by multiplying equation (20) by we obtain a polynomial of degree 2n + 1 in x 5 . Any zero of the polynomial, that lies between 0 and T , corresponds to a positive steady state. We argue again that, if x 5 ≥ T were a positive zero of the polynomial, equation (20) would give a contradiction. Therefore, any positive solution to the polynomial equation must fulfill that x 5 < T and hence provide a positive steady state.

Existence 2n + 1 positive steady states.
We have shown that the positive steady states of the nHK-Hpt system are determined by the positive solutions to (20). Solving for the positive solutions to this equation is equivalent to solving for the positive solutions to a polynomial of degree 2n + 1. By the fundamental theorem of algebra, a polynomial of degree 2n + 1 has 2n + 1 complex roots counted with multiplicity. Therefore, such a polynomial can at most have 2n + 1 distinct positive real roots. We show in this section that there exist choices of rate constants k * and total amounts H i , T such that the polynomial has exactly 2n + 1 distinct positive real roots. As a consequence, this proves that the nHK-Hpt system admits 2n + 1 positive steady states for some choice of rate constants and total amounts.
The proof consists of a series of simplifications and constructions. A key ingredient of the proof is the following theorem: Theorem 2 (Kurtz [5]). Let m ≥ 1 and let p(x) = x 2m+1 −c 1 x 2m +c 2 x 2m−1 +· · ·+c 2m x−c 2m+1 be a polynomial of odd degree 2m + 1 and with c i ≥ 0 for all i. Let c 0 = 1. If for all i = 1, . . . , 2m, then p(x) has 2m + 1 distinct positive real roots.
For clarity, we provide the main arguments of our proof here in the form of lemmas, which are proven in the next section.
First of all observe that the steady states of the system are invariant by multiplication of all rate constants by some scalar λ > 0. Therefore, we assume that k 6 = 1. For simplicity we write x for x 5 .
We let such that we write 5 .
As a consequence of Lemma 1, there exist rate constants and total amounts such that (20) holds if we can find α i,1 , . . . , α i,5 > 0 such that With this notation, we want to determine the positive real roots of the polynomial obtained by clearing denominators in (27): We would like to apply Theorem 2 to such a polynomial. To this end, the coefficients of monomials with even degree should be negative and those of odd degree should be positive. The latter is guaranteed if α i,4 = 0, while the former if α i,1 = 0. Setting these two constants to zero, for i = 1, . . . , n, does not change the degree of the polynomial.
By the continuity of the isolated roots of a polynomial as functions of the coefficients of the polynomial, if we can find α i,2 , α i,3 , α i,5 > 0 such that with α i,4 = 0, α i,1 = 0, the polynomial q(x) has 2n + 1 distinct positive real roots, then for α i,4 , α i,1 small enough, the polynomial q(x) still has 2n + 1 distinct positive real roots. This is what we do next. We set α i,4 = 0, α i,1 = 0 and further α i,3 = 1 for all i = 1, . . . , n, and T = 1. To ease the notation, we write a i = α i,2 and b i = α i,5 , such that the polynomial of interest is We denote by [n] = {1, . . . , n}. In the next lemma we describe the coefficients of p(x). The form of the coefficients depends on the parity of the degree of the coefficient. Therefore the coefficients take two different forms, one for even subindices, that is i = 2k, and one for odd subindices, that is i = 2k + 1.
For example, for n = 1 we have All that is left is to show that we can find b i , a i such that the polynomial p(x) satisfies the inequalities in Theorem 2. This is the content of the following lemmas. We provide in Lemma 3 a choice of constants b i such that the inequalities (21) are fulfilled for even indices i, that is, i = 2k for some k. In Lemma 4 we provide a choice of constants a i such that the inequalities (21) are fulfilled for odd indices i, that is, i = 2k + 1 for some k.  We are ready to prove the main result on the number of positive steady states.
Observe that the proof is constructive. It gives a procedure to find sets of parameters with the maximal number of steady states. The several checks that the proof requires are easily implemented using most available mathematical software to solve equations (e.g. Maple, Mathematica).
In general, we have observed that given any polynomial u(x) with 2n + 1 distinct positive real roots we can find a i , b i such that the coefficients of u(x) agree with c i in Lemma 2, even if u(x) does not fulfill the conditions of Theorem 2. Such a i , b i can be found using mathematical software.

n unstable steady states
In the subsection we show that, considering the 2n + 1 steady states ordered increasingly by their value x = x 5 , then the steady states number 2, 4, . . . , 2n are unstable relative to the stoichiometric compatibility class they belong to, that is, relative to the invariant subspaces described by the conservation laws (12) and (13).
Since the nHK-Hpt system has 4n + 2 variables and n + 1 conservation laws, the Jacobian of f inẋ = f (x) always has n + 1 zero eigenvalues. The remaining 3n + 1 eigenvalues (which could include zero) have corresponding eigenvectors in the stoichiometric subspace and dictate the dynamics around the steady state and within the stoichiometric compatibility class. If the steady state is locally stable relative to the stoichiometric compatibility class, then the product of these 3n + 1 eigenvalues has sign (−1) 3n+1 . Therefore, if the sign of the product of these eigenvalues is (−1) 3n , then the steady state is necessarily locally unstable relative to the stoichiometric compatibility class. We argue in the proof of the next theorem that this is the case for the steady states in even position 2, 4, . . . , 2n. Proof. We order the variables of the system as x 1,1 , x 1,2 , x 1,3 , x 1,4 , . . . , x n,1 , x n,2 , x n,3 , x n,4 , x 6 , x 5 . It follows from [6, Prop. 5.3] that the product of the 3n + 1 eigenvalues of the Jacobian with eigenvectors in the stoichiometric space agrees with the determinant of the Jacobian of the function g : R 4n+2 → R 4n+2 where for i = 1, . . . , n and We now apply the method described in [2] (see the Electronic Supplementary Material), to determine the sign of the determinant of the Jacobian of g from iterative eliminations. One can check that the expressions in (15)-(18) are obtained from iteratively eliminating x i,1 , . . . , x i,4 from the equations g 4(i−1)+1 (x) = · · · = g 4(i−1)+4 (x) = 0, which are equivalent to (14), and which correspond to the conservation law together withẋ i,2 =ẋ i,3 =ẋ i,4 = 0.
Let q(x 5 ) be the polynomial obtained after clearing denominators in (20). Then, by [2], the sign of the determinant of the Jacobian of g at a steady state agrees with the sign of the derivative of q(x 5 ), q (x 5 ), times (−1) 3n . Therefore, if q (x 5 ) is positive, then the corresponding steady state is locally unstable. Since q(0) is positive, the first real root of q(x 5 ) has negative derivative, and then the signs alternate. Therefore, the steady states corresponding to the 2, 4, . . . , 2n-th roots are locally unstable relatively to the stoichiometric compatibility class.

The core model for M hybrid phosphorelays competing for an RR
In the previous section we showed that unlimited multistationarity arises by increasing the number of hybrid HKs that compete for the same Hpt. In this section we show that the same statement holds when the competition occurs at the level of the response regulator RR. For example, we show that the system consisting of two core hybrid phosphorelays, complete with their own Hpt, competing for the same RR can have up to 9 steady states. This system consists of two 1HK-Hpt systems, where the two independent Hpt's donate their phosphate group to the same RR.
We study such a system with full generality. We allow each of the phosphorelays to have multiple hybrid HKs, as studied in the previous sections. For example, consider a system with two His-containing phosphotransfer proteins Hpt 1 and Hpt 2 that transfer the phosphate group to the same RR. Assume, for instance, that Hpt 1 receives the phosphate group from two hybrid HKs, HK 1,1 and HK 1,2 , and that Hpt 2 receives the phosphate group from one hybrid HK, HK 2,1 . The first upper index of HK i,j indicates the Hpt index, and the second index indicates the index of the HK in the nHK-Hpt subsystem. By using the notation introduced in the previous section to denote phosphorylated sites, the reactions of this example system are as follows: (i) Reactions within each HK:

Multiple positive steady states
In general, consider M multiple hybrid phosphorelays competing for the same RR. We let n i be the number of hybrid HKs of the ith system. That is, each phosphorelay consists of a n i HK-Hpt system with a further phosphotransfer to RR.
In the above example, we have M = 2 and n 1 = 2, n 2 = 1. We let HK i,j denote the jth hybrid kinase of the ith multiple HK-Hpt system, where j runs from 1 to n i , and we let Hpt i denote the Hpt of the system. We denote the reaction rates by k i,j, * for those involving HK i,j , the dephosphorylation reaction for Hpt i by k i,6 and the dephosphorylation reaction for RR by k 7 . The system has conserved total amounts of HK i,j , Hpt i and RR.
We denote the concentration of the species as follows: for i = 1, . . . , M , j = 1, . . . , n i . Under the law of mass-action, we model the dynamics of the concentrations over time by the following system of ordinary differential equations: for i = 1, . . . , M , j = 1, . . . , n i . Further, we have the following conservation law equations: For each i, j, solving forẋ i,j,1 = 0,ẋ i,j,2 = 0,ẋ i,j,3 = 0 and the conservation law with H i,j , we get equalities analogous to (15)-(18), where we replace the subindex i by the pair i, j. The steady state concentrations x i,j are expressed in terms of x i,5 and are positive provided x i,5 is positive.
Observe that both x 7 , x 8 are positive provided x i,5 < T i for all i. By substituting the expression for x 7 into (31), we deduce that the steady states of the system are found by finding the positive solutions to (31) in x 1,5 , . . . , x M,5 . The value at steady state of the other concentrations are found using the expressions above. Recall that, as seen in the previous section, a positive solution to (31) must satisfy x i,5 < T i and hence x i,6 is positive.
Theorem 5. Consider the system with M multiple hybrid phosphorelays competing for the same RR, and let n i be the number of hybrid HKs of the ith system. Then there exists a choice of rate constants and total amounts such that the system has M i=1 2n i + 1 positive steady states.
Proof. For each i = 1, . . . , M , fix parameters k i,j,1 , . . . , k i,j,5 , H i,j , T i such that the n i HK-Hpt system has 2n i + 1 steady states, when the dephosphorylation rate constant for Hpt i in the isolated system is set to one. By Theorem 3 such a choice exists. With this choice, let , which only depends on x i,5 . Consider the map Therefore, the Jacobian of ϕ(0, x 1,5 , . . . , x M,5 ) at a solution is nonsingular. This allows us to apply the implicit function theorem to ensure that for S > 0 small enough, the equation ϕ(0, x 1,5 , . . . , x M,5 ) = 0 has precisely M i=1 2n i + 1 positive solutions. Fix any such S > 0 and define k i,6 = S, R = 1 S , and k 7 = 1. Then , which corresponds to the steady state equation (31). With this choice, the system has M i=1 2n i + 1 positive steady states.
The proof gives a constructive way to find sets of parameters with M i=1 2n i + 1 positive steady states. We fix parameters for the individual nHK-Hpt systems that have 2n i + 1 steady states, when the dephosphorylation rate constant for Hpt i in the isolated system is set to one. Then, let k 7 = 1, k i,6 = 1/R and increase R until the system has the desired number of steady states.

From core to full models
In the previous sections we have analysed several systems based on a core model. That is, we have disregarded reverse reactions, hydrolysis reactions, and even complex formation in the phosphotransfer reactions. We provide here arguments that guarantee that the properties on the number of steady states of the different core models that we have considered extend to the full models.

Theoretical results
The results concerning the number of steady states extend to the full models. Two mathematical results, valid for mass-action kinetics, are used for this claim ( [3], [4]): [3] Assume that a network has N (nondegenerate) positive steady states. If complex formation is taken into account, that is, a reaction is split into two by adding an intermediate, then the new extended network also has N (nondegenerate) positive steady states for some choice of rate constants and total amounts. [4] Assume that a network has N (nondegenerate) positive steady states. If reactions are added to the network, in such a way that the conservation laws of the system are preserved, then the new network also has N (nondegenerate) positive steady states for some choice of rate constants and total amounts.
The nondegeneracy condition means that the Jacobian is nonsingular relative to the stoichiometric compatibility class described by the conservation laws (cf. Subsection 2.4). This requirement is fulfilled in our case.

Full hybrid HK
The full model of a hybrid HK with reversible reactions and hydrolysis reactions, consists of the reactions: One might also consider complex formation, that is, substitute each phosphotransfer reaction of the form By Subsection 3.2, the core phosphorelay with a hybrid HK can have 3 positive steady states. Adding reversibility to some reactions does not change the conservation laws, nor does including hydrolysis reactions. Therefore, by [4], the full phosphorelay model with a hybrid HK can have 3 positive steady states for some choice of rate constants and total amounts. This holds true even if hydrolysis reactions are added on the other phosphorylation sites (that is, the first domain of the hybrid HK and Hpt). By [3], adding complex formation also maintains the maximal number of steady states. In both cases, however, a higher number of steady states might be achievable.
The same argument holds for all the models considered in the previous sections.
Using the first equation in (33), we have Finally, using the second equation in (33), we have , which is equivalent to Fix any k 2 > 0. Since the polynomial (37) is a polynomial in k 1 , with positive leading coefficient and negative independent term, it has a unique positive real root. As a consequence, for any k 2 > 0, (33) and (34) hold with k 1 defined such that (37) holds, and k 3 , k 4 , k 5 , H > 0 as in (35), (36).

Proof of Lemma 3
To simplify the presentation, we denote by P k (n) the set of all subsets of [n] with k elements. Similarly, given i ∈ [n], we denote by P k−1 (n, i) the set of all subsets of [n] with k − 1 elements that do not contain i. We want to show that if we define b i = We let Then, to show that c 2 2k > 4c 2k+1 c 2k−1 it is enough to show that B 1 4 ≥ B 2 2 . To this end, observe that given any pair of sets L ∈ P k (n) and J ∈ P k−1 (n), such that J ⊆ L, we can choose i ∈ L such that L := L \ {i} satisfies L = J. Therefore, every summand of B 2 is a summand of B 1 as well. This finishes the proof.

Proof of Lemma 4
Let δ k = c 2 2k+1 − 4c 2k c 2k+2 . For M > 0 and a 1 > 0, we let a i = an M i−1 and b i = a 2 i 4 . Therefore, we have that b i = a 2 1 4M 2(i−1) , for i = 1, . . . , n. With these substitutions, δ k = δ k (a 1 ) becomes a polynomial in a 1 . If there exists a choice of M > 0 such that the coefficient of smallest degree of δ k is positive for all k, then δ k (a 1 ) > 0, for a 1 small enough and all k. Therefore, we compute the coefficient of smallest degree of δ k in a 1 .
We use the notation introduced in the proof of Lemma 3. By definition, we have that The polynomial c 2 2k+1 (a 1 ) consists of one term of degree 4k in a 1 . The polynomial c 2k (a 1 ) is a sum of a term of degree 2k and one of degree 2k − 1. Similarly, the polynomial c 2k+2 (a 1 ) is a sum of a term of degree 2k + 2 and one of degree 2k + 1. Hence, the product c 2k c 2k+2 is a polynomial with lowest degree 4k. If the coefficient of degree 4k of δ k is nonzero, then it is the coefficient if smallest degree. By denoting the coefficient of degree 4k by β k , we have: