Global stability of a class of futile cycles

In this paper, we prove the global asymptotic stability of a class of mass action futile cycle networks which includes a model of processive multisite phosphorylation networks. The proof consists of two parts. In the first part, we prove that there is a unique equilibrium in every positive compatibility class. In the second part, we make use of a piecewise linear in rates Lyapunov function in order to prove the global asymptotic stability of the unique equilibrium corresponding to a given initial concentration vector. The main novelty of the paper is the use of a simple algebraic approach based on the intermediate value property of continuous functions in order to prove the uniqueness of equilibrium in every positive compatibility class.


Introduction
We consider the following class of futile cycles: We give a few words of explanation regarding the above class of futile cycles. There are n reaction chains shown in (1). It is assumed that the kth reaction chain has the first m k reactions reversible and the last reaction irreversible as shown in (1). For the ith reaction chain (i = 1, . . . , n), C i j ( j = 1, . . . , m i ) denote the intermediate complexes and E i denotes the enzyme catalyzing the generation of the product P i from the substrate P i−1 , assuming that P n and P 0 denote the same compound. It is assumed that each of the reactions in (1) is governed by mass action kinetics. Thus each reaction chain in (1) is a reaction mechanism of an enzyme catalyzed single substrate, single product reaction. Since the product of the last reaction chain is the substrate of the first one, (1) is a futile cycle. For i = 1, . . . , n, j = 1, . . . , m i + 1, k i j denotes the forward reaction constant and for i = 1, . . . , n, j = 1, . . . , m i , k −i j denotes the reverse reaction constant of the jth reaction in the ith reaction chain.
One example of a futile cycle model in biochemistry of the type (1) is the model of processive multisite phosphorylation system discussed in Conradi and Shiu (2015). Phosphorylation systems play a vital role in cellular signalling, and their stability properties are of great interest. For a detailed exposition on phosphorylation systems, the reader is referred to Salazar and Høfer (2009). There are two primary mechanisms and related mathematical models of multisite phosphorylation systems. These are distributive and processive mechanisms. The processive mechanism of multisite phosphorylation system discussed in Conradi and Shiu (2015) looks as follows: Note that this mechanism is a specific case of the kind of futile cycles (1) under consideration in this paper.
It has been shown (Wang and Sontag 2008) that a distributive multisite phoshorylation mechanism can exhibit the property of multistationarity if there are at least two sites of phosphorylation, meaning that there can be more than one equilibrium in every positive compatibility class. However, it has been shown (Conradi and Shiu 2015) that a processive multisite phosphorylation system admits a unique equilibrium in every positive compatibility class. We explain in simple words what this means. Any chemical reaction network has a set of conservation laws. For example, it can be shown that the futile cycle (1) under consideration in this paper has the sum of the concentrations of all the substrates P i (i = 0, . . . , n − 1) and the intermediate complexes C i j conserved and positive. Also it can be easily proven that for the ith reaction chain, the sum of the concentrations of the enzyme E i and all the intermediate complexes C i j ( j = 1, . . . , m i ) involved in the ith reaction chain is conserved and positive. If the set of initial concentrations of all the substrates, enzymes and intermediate complexes are given, the reactions will evolve in such a way that the conserved quantities remain conserved. When we say that there is a unique equilibrium in every positive compatibility class, it means that for a fixed set of values of the conserved quantities of a chemical reaction network, there exists only one positive vector of equilibrium concentrations that has the same values of the conserved quantities. Note that equilibrium concentrations are positive concentrations of the species of a reaction network at which the rates of change of the species concentrations are all zeros.
Further, it is shown in Conradi and Shiu (2015) that futile cycles of the type (2) governed by mass action kinetics are globally asymptotically stable. This has been proved using monotone systems theory (Angeli and Sontag 2008). In this paper, we prove the same result for a larger class of futile cycles given by (1). We first prove that (1) has a unique equilibrium in every positive compatibility class. This is proved in a simple way by first of all noting that the equilibrium concentrations of all the species are related, and hence one can obtain a polynomial equation in terms of the equilibrium concentration of one of the species in the network. We then make use of intermediate value property of continuous functions in order to prove the uniqueness of a positive equilibrium of this species and hence of all the other related species. For proving global asymptotic stability of futile cycle (1), we construct a Piecewise Linear in Rates (PWLR) Lyapunov Function. Such Lyapunov functions have been used to prove stability of chemical reaction networks recently in Blanchini and Giordano (2014) and Al-Radhawi and Angeli (2016). A piecewise linear Lyapunov function in the timederivative of the states was also used to prove stability of nonlinear compartmental systems in Maeda et al. (1978) and Maeda and Kodama (1979).

Notation:
The space of n dimensional real vectors is denoted by R n , and the space of m × n real matrices by R m×n . The space of n dimensional real vectors consisting of all strictly positive entries is denoted by R n + . ker(A) and im(A) denote the kernel and image of a real matrix A. A denotes the transpose of a matrix A. y i denotes the ith component of a vector y. 0 p,q denotes a matrix of dimension p × q whose entries are all zeros. span(v) denotes the span of a real vector v.

Stoichiometry, positive compatibility class and conservation laws
In this section, we explain how to model the stoichiometry and how this can be used to derive the conservation laws of a chemical reaction network. Let us assume that a given chemical reaction network has ρ reactions which could be a mix of reversible and irreversible chemical reactions. Let r denote the vector of reaction rates in the forward direction of the ρ reactions. Let us assume that the reaction network involves μ chemical species and let x denote the vector of their concentrations. The dynamical equations for the chemical reaction network can be then obtained from the balance lawsẋ where S ∈ R μ×ρ is called the stoichiometric matrix. Equation (3) holds for any chemical reaction network irrespective of the governing laws of its reactions. However r is a function of x and depends on the governing laws of the reactions in the network. We now show with the help of a simple example how the stochiometric matrix of a network can be constructed. Consider the reaction network The stoichiometric matrix maps the space of reaction rates to the space of species concentrations. Consequently, the entry of S corresponding to the ith row and jth column is the difference between the number of moles of the ith species on the left hand and the right hand sides of the jth reaction. For the reaction network (4), S is given by The network is at equilibrium whenẋ = 0. Let r * denote the vector of reaction rates (in the forward direction) of the network at an equilibrium. Then from Eq. (3), it follows that r * ∈ ker(S). Thus for the reaction network (4), r * 1 = 0 and r * 2 = 0. Let x 0 ∈ R μ + denote the vector of initial concentrations, i.e., x 0 = x| t=0 . Then it is easy to see from Eq. (3) that The space of concentrations has been referred to as the positive reaction simplex in Horn and Jackson (1972) and the positive stoichiometric compatibility class (corresponding to x 0 ) in Feinberg (1995), Anderson (2011) and Siegel and MacLean (2000). In this paper we will use the terminology positive stoichiometric compatibility class to refer to S x 0 .
Consider a vector k ∈ ker(S ). This vector obeys the equation k ẋ = 0, consequently k x is a constant, meaning that k x = k x 0 . The equation is called a conservation law of the network. In general the number of linearly independent conservation laws is equal to the dimension of ker(S ).
Notice that the vector x stays in S x 0 as long as x ∈ R m + and k x = k x 0 for all vectors k ∈ ker(S ).
Thus irrespective of the governing laws of a reaction network, based on the stoichiometry of the network, on the one hand, we can find the relation between reaction rates of the various reactions of the network at equilibrium by computing ker(S), and on the other hand, we can find all the conservation laws by computing ker(S ). Below, we examine the former for the case of the reaction network (1) of interest.

Lemma 1 The futile cycle network (1) is at equilibrium if and only if the reaction rates (in the forward direction) of all the reactions of the network are equal.
Proof To prove the lemma, we first need to construct the stoichiometric matrix for the network. The entries of the stoichiometric matrix depend on the ordering of the species and the reactions of the network. Let us order the species as in the following set.
{P 0 , P 1 , . . . , P n−1 , E 1 , E 2 , . . . , E n , C 11 , C 12 , . . . , C 1m 1 , C 21 , . . . , C 2m 2 , . . . , . . . , C n 1 , . . . , C nm n } Let R i j denote the jth reaction in the ith reaction chain of the network and let r i j denote its rate. For the construction of the stoichiometric matrix, let us order the reactions of the network as in the following set.
Let e j denote the jth element of the standard basis of the vector space R n . For i = 1, . . . , n − 1, define and define P n := −e n 0 n,m n −1 e 1 .
For i = 1, . . . , n, define For i = 1, . . . , n, define C i as a matrix of dimension m i × (m i + 1) given by The stoichiometric matrix S for the network (1) is then given by In order to prove the lemma, it suffices to prove that ker(S) = span(1), where 1 is a vector whose dimension is equal to the number of reactions in the network and each of whose entries is equal to 1. Let r * i j ( j = 1, . . . , m i + 1) denote the value of r i j at an equilibrium. Considering the ordering of reactions in (5), we have r * = [ r * 11 r * 12 · · · r * 1(m 1 +1) r * 21 · · · r * 2(m 2 +1) · · · · · · r * n1 · · · r * n(m n +1) ] ∈ ker(S) It follows that Cr * = 0, which implies that In words, in each chain of reactions in the network, the rates of all the individual reactions at equilibrium are equal. Since Pr * = 0, we get . . .
Combining the sets of Eqs. (7) and (8), we get r * ∈ span(1). It is easy to verify that Ev = 0 for any v ∈ span(1), since every row of E has one entry equal to 1, another entry equal to −1 and the remaining entries equal to 0. This proves that ker(S) = span(1). Hence the proof.
Using a similar method as in the proof above, one can also compute ker(S ) with S given by Eq. (6) and prove that the dimension of the ker(S ) is n + 1. For i = 0, . . . , n−1, letp i denote the concentration of P i and for i = 1, . . . , n, letẽ i denote the concentration of E i . For i = 1, . . . , n, j = 1, . . . , m i , letc i j denote the concentration of the complex C i j . Then the equations with real constants α and β i (i = 1, . . . , n), are n + 1 linearly independent conserved quantities for the network.

Law of mass action kinetics
Recall that the reaction rates depend on the governing laws of the reactions of the network. We now describe this relation for mass action reaction networks, which are reaction networks, each of whose reactions is governed by the law of mass action kinetics. According to this law, the rate of a reaction is proportional to the concentrations of the different species on the substrate side of the reaction. For example, consider the reaction In the reaction above, k f and k r are positive constants known as the forward and the reverse rate constants. Letx i denote the concentration of X i for i = 1, 2, 3. The mass action reaction rate of the forward reaction is k fx1x2 , and the rate of the reverse reaction is k rx3 . Therefore the overall reaction rate in the forward direction of the reversible reaction (11) is r = k fx1x2 − k rx3 . Now consider the reversible reaction The above reaction can also be written as Ifỹ i denotes the concentration of the Y i for i = 1, . . . , 4, then the overall mass action reaction rate in the forward direction of the reversible reaction (12) is given by r = k fỹ 2 1ỹ 3 2 − k rỹ3ỹ 2 4 . Now consider the ith reaction chain in the futile cycle network (1) Assuming that the governing law of each reaction is mass action kinetics, with r i j as defined earlier, we get Using Eqs. (3), (6) and (13), we get the following set of equations describing the rate of change of concentrations of the different species involved in the network (1):

Uniqueness of equilibrium
In this section, we prove that the futile cycle (1) has a unique equilibrium in every positive compatibility class. This proof requires that the nonnegative orthant is invariant with respect to the dynamical equations corresponding to the futile cycle network (1). This property commonly referred to as nonegativity is well known for all mass action reaction networks. The reader is referred to (Vol'pert and Hudjaev 1985, pp. 613-615) for a simple proof of nonnegativity of mass action chemical reaction networks. For the sake of completeness, we provide an alternate proof of nonnegativity specifically for the case of the network of interest.
Lemma 2 If the initial concentration of each of the species of the futile cycle network (1) governed by mass action kinetics is nonnegative, then it remains nonnegative at all future times.
Proof Note that all the forward and the reverse rate constants of the reactions in the network are positive. If we assume that at a particular instant of time, one of the species of the network has zero concentration and the rest have nonnegative concentrations, then from Eq. (14), it follows that the time derivative of the species with zero instantaneous concentration is nonnegative. This proves that the nonnegative orthant is invariant with respect to the differential equations (14), since any point on the boundary of the nonnegative orthant will either be pushed into the positive orthant or will stay in the boundary. In other words, if the initial concentrations are nonnegative, they remain nonnegative at all times.
Hereafter in all the analysis with respect to the mass action chemical reaction network (1) that follows, it will be assumed that the initial concentrations of P 0 , E i , i = 1, . . . , n are positive and the initial concentrations of the remaining species of the network are nonnegative. This ensures that with respect to the conservation laws (9) and (10), we have α > 0 and β i > 0 for i = 1, . . . , n. Further p i , e i and c i j will be used to denote the concentrations of P i , E i and C i j respectively at an equilibrium.
In order to prove the uniqueness of equilibrium in each positive compatibility class, we need to be able to parametrize the equilibrium concentrations of all the species of the network in terms of the equilibrium concentration of one particular species. Below we prove that such a parametrization is possible.
Lemma 3 Choose a q ∈ {1, 2, . . . , n}. For each i ∈ {1, . . . , n}, there exist constants a i ∈ R + , b i ∈ R + and δ i j ∈ R + , j ∈ {1, . . . , m i }, such that where k i j denotes the mass action forward rate constant of the jth reaction in the ith reaction chain of the futile cycle network (1).
Proof Lemma 1 states that at equilibrium, the overall reaction rate in the forward direction of each of the reversible reactions and the reaction rates of each of the irreversible reactions in (1) are equal to each other. This implies by virtue of equations (13) that for i = 1, . . . , n and k 1(m 1 +1) c 1m 1 = k 2(m 2 +1) c 2m 2 = · · · = k n(m n +1) c nm n .
We now prove that the futile cycle (1) cannot have any boundary equilibrium. (1), none of the species concentrations is zero.

Lemma 4 At an equilibrium of the mass action futile cycle network
Proof We use the parametrization of Lemma 3 for the proof. As mentioned earlier, with reference to the conservation laws (9) and (10), α > 0 and β i > 0 for i = 1, . . . , n.
If we assume that c i j = 0 for some admissible numbers i, j, then from Eqs. (15) and (16), it follows that the equilibrium concentrations of all substrates and complexes are zeros. Since the conservation law (9) holds also for equilibrium concentrations, it follows that α = 0, which leads to a contradiction. Assuming that p i = 0 for some i ∈ {0, . . . , n − 1} leads to a similar contradiction. If we assume that e i = 0 for some i ∈ {1, . . . , n}, then from Eq. (17), it follows that From Lemma 2 and conservation law (9), we have 0 ≤ p i−1 < α. Owing to the finiteness of p i−1 , from Eq. (16), it follows that c qm q = 0, which leads to a similar contradiction as before. Hence, the futile cycle (1) cannot have any boundary equilibrium.
Remark 5 Lemma 4 can be also proved using (Vol'pert and Hudjaev 1985, Theorem 1, p. 617) related to positivity of reachable species in a chemical reaction network. This method has been used in Siegel and MacLean (2000) in order to prove the nonexistence of boundary equilibria in enzyme kinetic reaction networks that are similar to but simpler than the network (1). We briefly explain the concept of reachability as in Vol'pert and Hudjaev (1985). Assume that α species S 1 , . . . , S α of a network have nonzero initial concentrations and the rest have zero initial concentrations. Let A 0 denote the set of species with nonzero initial concentrations.
Now consider all the reactions whose substrates consist only of species contained in A 0 . Let A 1 denote the set of all those species that are products of such reactions. Note that in this step, the forward and the reverse parts of each reversible reaction must be considered as two separate reactions. Now update the set A 0 as A 0 = A 0 ∪ A 1 and repeat the process described below Eq. (26) until A 0 cannot be updated with new species anymore. At this point, all the species that are not contained in A 0 are said to be unreachable from the species in the initial set (Vol'pert and Hudjaev 1985, Theorem 1, p. 617) states that all species that are unreachable from the set (27) have zero concentrations at all times. If we assume that for the network (1) under consideration, we begin with equilibrium concentrations of species, then all species must remain at the same equilibrium concentrations. Now assume that c i j = 0 for some admissible i and j. Then by (Vol'pert and Hudjaev 1985, Theorem 1, p. 617), C i j is unreachable from the set of species with nonzero equilibrium concentrations. Thus c ik = 0 for k = 1, . . . , m i and either e i = 0 or p i−1 = 0. e i = 0 violates the conservation law (10), since it is assumed that α > 0, β i > 0, i = 1, . . . , n. Hence p i−1 = 0. This implies that c η j = 0 for j = 1, . . . , m η , where η = i − 1 if i = 1 and η = n otherwise. This in turn implies that p l = 0 and c lk = 0 for all admissible l and k. However this violates conservation law (9), since α > 0. Thus all intermediate complexes C lk have nonzero equilibrium concentrations. In a similar way, it can be proved that p i−1 = 0 and e i = 0 for i = 1, . . . , n.
We now prove the main result of this section pertaining to the uniqueness of equilibrium in every positive compatibility class.

Theorem 6
The mass action chemical reaction network (1) has a unique equilibrium in every positive compatibility class.
The above equation is in terms of only one unknown c qm q . Consider the numbers Since the choice of q in Lemma 3 is arbitrary, let q denote the index for which d q is minimum among the elements of {d 1 , d 2 , . . . , d n }. Now define the positive constants . . . ,n}. Substituting (29) in (28), we get From Lemmas 2 and 4, it follows that the equilibrium concentration of each of the species is positive. Thus e q > 0, and consequently from Eq. (17), it follows that c qm q < b q . We now prove that Eq. (30) in the unknown c qm q has precisely one real root in the admissible open interval (0, b q ). The proof is based on the intermediate value property of continuous real valued functions. Define Assume that v i = b q for exactly k values of i ∈ {1, 2, . . . , n}\{q}. For the remaining values of i ∈ {1, 2, . . . , n}\{q}, v i > b q . Note that 0 ≤ k ≤ n − 1. Without loss of generality assume that v i = b q for i = n − k + 1, n − k + 2, . . . , n and q = n − k. Define and note that ν is a continuous function in the closed interval [0, b q ]. Observe that This implies that there is at least one root of ν in the interval (0, b q ). We prove that there is exactly one root in this interval. Assume by contradiction that there are at least two roots x 1 , x 2 of ν in the interval (0, b q ). Then ν(x 1 ) = 0 and ν(x 2 ) = 0. Since v i ≥ b q for i ∈ {1, 2, . . . , n − k}, this implies that Subtracting Eqs. (32) from (31), we get Notice that (33), (34) and (35) lead to a contradiction. A similar contradiction can be obtained if we assume that x 2 > x 1 . This implies that x 1 = x 2 , in other words, there is exactly one root of ν in the interval (0, b q ). This in turn implies that there is exactly one root of in the interval (0, b q ). It follows that Eq. (30) in the unknown c qm q has exactly one real root in the admissible open interval (0, b q ). This implies that in a particular positive compatibility class, the equilibrium concentration c qm q is unique and positive. From Lemma 3, it follows that the equilibrium concentrations of all the remaining species of the futile cycle are unique and positive in a particular positive compatibility class (with α, β i > 0).

Global stability
We are now in a position to prove the global stability of futile cycle network (1). As mentioned in the introduction, for this purpose, we make use of a PWLR Lyapunov function coupled with LaSalle's invariance principle (LaSalle 1960;Khalil 2014, Section 4.2, Murray et al. 1994. Theorem 7 All solution trajectories to the set of Eq. (14) that belong to a particular positive stoichiometric compatibility class asymptotically converge to the unique equilibrium in that stoichiometric compatibility class.
Proof Note that the set of Eq. (14) describes the rate of change of concentrations of the different species in the futile cycle network (1). With respect to this network, let r denote the vector of all the reaction rates r i j as defined in Sect. 3. We now define the PWLR Lyapunov function It is easy to see that V is a piecewise linear function of the components of r . Observe that V (r ) is continuous in time, since each component of r is continuous in time. We consider three possible cases below and prove that d dt V (r ) ≤ 0 for each of the three cases. Case 1: r ∞ = |r i1 | for some i ∈ {1, . . . , n} in the time-interval t 0 < t ≤ t 1 .
As in (13), we have In case, r i1 > 0, Since r i1 > 0 and r ∞ = r i1 , it follows with the equality occuring if r i1 = r im i +1 = r i−1m i +1 = r i2 . In case r i j < 0, we have In this case r ∞ = −r i1 and hence again In case r i1 = 0, we have r = 0, which implies (from Eq. (13)) that the concentrations of each of the species in the network is zero, which is not possible since it violates the conservation laws (9), (10). Thus Case 2: r ∞ = |r im i +1 | for some i ∈ {1, . . . , n} in the time-interval t 0 < t ≤ t 1 . As in (13), r i(m i +1) = k i(m i +1)cim i ≥ 0. We have r i(m i +1) = 0, since r i(m i +1) = 0 implies that r = 0 leading to all the concentrations being equal to zero as described earlier. Hence r i(m i +1) > 0. Since the rate of the reaction is greater than or equal to that of any other reaction in the network, in particular of the reaction Case 3: r ∞ = |r i j | with 1 < j < m i + 1. As in (13), we have As in the other cases, we cannot have r i j = 0. In both the cases, r i j > 0 and r i j < 0, we can prove that d dt V (r ) = d dt |r i j | ≤ 0 as in Case 1 with the equality occuring if r i( j−1) = r i j = r i( j+1) . Thus Let E denote the set of all vectors r for which d dt V (r ) = 0. In all the three cases considered above, it can be proved that any r ∈ E that is positively invariant satisfies r ∈ span(1) where, as before 1 denotes a vector with the same dimension as r that has every entry equal to 1. In the following, we prove this for Case 1. For the remaining 2 cases, the proof is similar.
From Lemma 1 and Theorem 6, it follows that for a given initial concentration vector x 0 , the only vector r ∈ E that is positively invariant corresponds to the unique equilibrium concentration x * ∈ S x 0 , where S x 0 denotes the positive compatibility class corresponding to x 0 . Since d dt V (r ) ≤ 0 and r is a continuous function of the concentration vector x, by LaSalle's invariance principle, it follows that every x ∈ S x 0 asymptotically converges to x * .

Conclusion
We have proved the global asymptotic stability of a special class of futile cycle networks which includes the model of processive multisite phosphorylation networks described in Conradi and Shiu (2015). We used a simple algebraic method based on intermediate value property of continuous functions in order to prove the existence of a unique equilibrium concentration vector in every positive compatibility class. Furthermore, we constructed a piecewise linear in rates (PWLR) Lyapunov function in order to prove the global asymptotic stability of this equilibrium vector. This construction is inspired by the recent work in Al-Radhawi and Angeli (2016), where a systematic approach towards the construction of PWLR Lyapunov functions in order to prove stability of certain chemical reaction networks is described. Current research is focused on the direction of using the special Lyapunov function described in this paper to prove stability of a larger class of chemical reaction networks.