A hidden integral structure endows Absolute Concentration Robust systems with resilience to dynamical concentration disturbances
AA hidden integral structure endows Absolute Concentration Robustsystems with resilience to dynamical concentration disturbances
Daniele Cappelletti ∗ , Ankit Gupta ∗ and Mustafa Khammash ∗ Abstract
Biochemical systems that express certain chemical species of interest at the same level at any positive equilib-rium are called “absolute concentration robust” (ACR). These species behave in a stable, predictable way, in thesense that their expression is robust with respect to sudden changes in the species concentration, regardless thenew positive equilibrium reached by the system. Such a property has been proven to be fundamentally importantin certain gene regulatory networks and signaling systems. In the present paper, we mathematically prove thata well-known class of ACR systems studied by Shinar and Feinberg in 2010 hides an internal integral structure.This structure confers these systems with a higher degree of robustness that what was previously unknown. Inparticular, disturbances much more general than sudden changes in the species concentrations can be rejected, androbust perfect adaptation is achieved. Significantly, we show that these properties are maintained when the systemis interconnected with other chemical reaction networks. This key feature enables design of insulator devices thatare able to buffer the loading effect from downstream systems - a crucial requirement for modular circuit designin synthetic biology.
The network of chemical interactions of a biochemical system of interest can be complex and involve unknown reactionpropensities. One of the main goal of reaction network theory consists in deriving dynamical properties from simplergraphical characteristic of the model, and independently on the specific value of kinetic parameters [22, 45]. Theresults presented in this paper follow this approach.A qualitative property of great interest is the capability of a certain chemical species to be expressed with thesame concentration at any positive steady state, independently on the initial conditions and on how many steadystates are present. Namely, assume that the dynamics of the biochemical system are expressed by the d − dimensionalordinary differential equation ddt x ( t ) = f ( x ( t )) . We say that the i − th species is absolute concentration robust (ACR), if there exists an ACR value q independentof the initial condition x (0) such that, whenever x ( t ) tends to a positive vector x , we have x i = q . In the typicalcases of interest, the positive steady state x that is reached will depend on the initial condition x (0), while the entry x i = q does not. As noted in [5], the property of absolute concentration robustness alone does not imply stability ofthe positive steady states: it only ensures that if a positive steady state x exists, then the value of the ACR speciesat x is the ACR value.Under the assumption of stability, absolute concentration robustness provides a reliable, predictive response toenvironmental changes, since the species of interest reaches the equilibrium level relative to the new environmentalsetting, regardless the previous conditions. The existence and importance of this robustness property for various generegulatory networks and signal transduction cascades is explored in many papers, including [1, 10–12, 35, 37, 39, 42,43].In the Control Theory community, and under the assumption of stability, the absolute concentration robustnessis known as “robustness to disturbances in the initial conditions” [46, 47]. To achieve robustness with respect tosome disturbance, the imbalance caused by the disturbance needs to be measured first. To this aim, a quantity ofinterest is the integrator , which is a function φ of the system variables whose derivative is exactly the imbalance to beeliminated. At steady state, the derivative of φ is zero and so needs to be the imbalance. In the setting of absoluteconcentration robustness, one would like to find an integrator whose derivative is the difference between the ACRspecies and its ACR value. Unfortunately, in general this cannot be done, as shown in Section 4.1 and as discussedin [46].In the present paper, we systematically study for the first time the connection between ACR systems and in-tegrators. Specifically, our first contribution is related to the existence of a linear combination of chemical species ∗ Department of Biosystems Science and Engineering, ETH Zurich, Mattenstrasse 26 4058 Basel, Switzerland. a r X i v : . [ q - b i o . S C ] O c t hose derivative is the difference between the ACR species and its ACR value, multiplied by a monomial. Such linearcombination of species is called constrained integrator (CI), because it behaves similarly to an integrator given thatthe monomial does not vanish [46]. We rigorously prove that such a linear CI always exists for a large class of modelsthat strictly includes the ACR systems introduced in [40]. This result has some important consequences: first of all,under the assumption of stability, it implies that the expression of ACR species is not only robust to changes in theinitial conditions, but also to disturbances that are applied over time.An important application in Synthetic Biology concerns the design of insulators . A number of biochemical systemsare known to express a specific output if given a certain input. The systems can therefore be considered as moduleswith different functions. In cells, different modules are combined so that more complex responses to external stimulibecome possible [30]. In Synthetic Biology, it is desirable to combine different modules to achieve the same level ofcomplexity [38]. However, when connected, the different modules can affect the dynamics of each other and theycan lose the desired dynamical properties they had when considered in isolation [17]. In a simplified framework, an upstream module processes an input, and its output is fed to a second, downstream module to be further processed.Since the information is passed in form of molecules, which are then consumed or temporarily sequestrated by thedownstream module, the equations governing the upstream module dynamics are perturbed and its functionality canbe affected. Such effect is commonly called loading effect [34] or retroactivity [17, 36], and needs to be minimized. Inother words, the upstream module needs to be insulated from the loading effect caused by the downstream module.We propose two ways in which the robustness of the systems studied in this paper can be used to this aim. Thefirst solution is to simply design an upstream module which is robust to loading effects, modeled as a persistentdisturbance over time. The second solution is to design an extra component, called insulator , which transfers thesignal from the upstream module to the downstream module while at the same time shielding the dynamics of theupstream module from retroactivity effects.We will also show how more theoretical results on reaction network models can be obtained as a consequenceof our work. In Reaction Network Theory, the study of steady state invariants consitute an interesting topic ofresearch[18, 22, 33, 40]. In [40], it has been proven that certain graphical properties of the network imply theexistence of an ACR species, regardless the choice of kinetic parameters. Such sufficient conditions are generalized inthe present work while they are maintained simple to check. Moreover, no way to explicitly determine the ACR valuewas given in [40], and we fill the gap by proposing a fast linear method to calculate it. Furthermore, a substantialeffort in the Reaction Network community is devoted to understand under what conditions dynamical properties ofsingle systems can be lifted to larger systems [9, 24, 25, 29, 32]. Our contribution in this sense consist in provingthat, under certain conditions, if an ACR system of the class studied in this paper is part of a larger model, the ACRspecies is still so in the lager system and its ACR value is maintained. Finally, it is worth mentioning that in thepresent work we consider the possibility of time-dependent rates for the occurrence of chemical transformations. Thisis more general than what is usually studied in Reaction Network Theory, with the exception of few works explicitlyallowing for this scenario [2, 13, 14, 16, 28, 31]. Consider two proteins A and B , whose interaction is described by A + B κ BB κ A (1)where the positive constants κ and κ describe the propensity of a reaction to occur. If enough proteins are presentand they are homogeneously spread in space, then a good model for the time evolution of the concentrations ofproteins A and B is given by mass-action kinetics. Specifically, the concentrations of A and B at time t , denoted by x A ( t ) and x B ( t ) respectively, are assumed to solve ddt (cid:18) x A ( t ) x B ( t ) (cid:19) = κ x A ( t ) x B ( t ) (cid:18) − (cid:19) + κ x B ( t ) (cid:18) − (cid:19) (2)It is easy to check that the steady states of (2) are given by states ( x A , x B ) such that either x B = 0 or x A = κ /κ .Hence, A is an ACR species because its expression at any positive steady state is the same. It is common duringbiochemical experiments to be able to control the inflow rate of some species (say B ). Some additional chemicalspecies may also be introduced, with the purpose of degrading some of the present components (in this case, species2 is introduced to faster degrade species B ). After these modifications, (1) becomes A + B κ B u ( t ) B κ A u ( t ) κ CC + B κ ddt x A ( t ) = − κ x B ( t ) (cid:18) x A ( t ) − κ κ (cid:19) , it is still true that the value of x A ( t ) will converge to κ /κ , as long as the functions u , u and the rates κ , κ are such that the species B is not removed from the system too fast. In this paper we will prove a general resultdescribing when such kind of robustness to persistent perturbations is present for ACR systems. Consider the system described in Figure 1. The model is proposed in [40, 41] as osmoregulatory system in EscherichiaColi. It is in accordance with experimental observations discussed in [11, 37, 43]. According to the model, whoseschematics is described in Figure 1, the activation rate of the sensor-transmitter protein EnvZ depends on the mediumosmolarity. Then, an active form of EnvZ transfers its phosphoryl group to the sensory response protein OmpR,which becomes OmpR-P and promotes the production of the outer membrane porins OmpF and OmpC. Hence, itis important that the concentration of OmpR-P responds in a reliable, predictive way to changes in the mediumosmolarity (which the rate constants κ i in the reaction network of Figure 1 depend upon), but not on the initialconcentration of the different chemical species involved. As a matter of fact, it follows from the results developed in[40] that OmpR-P is an ACR species.EnvZ- D κ κ [ADP] EnvZ κ [ATP] κ EnvZ- T κ EnvZ- P EnvZ- P + OmpR κ κ EnvZ-OmpR- P κ EnvZ + OmpR- P EnvZ- D + OmpR- P κ κ EnvZ-OmpR- D - P κ EnvZ- D + OmpRFigure 1: Proposed model for the EnvZ-OmpR signal transduction system in Escherichia Coli, which is able to explainthe experimentally observed robustness in the expression of phosphorylated OmpR. In the first line of reactions, EnvZcan bind to ADP and ATP, but only when bound to ATP it can gain a phosphoryl group, and the resulting speciesis denoted by EnvZ-P. In the second line of reactions, EnvZ-P transfers the phosphoryl group to OmpR, through theformation of an intermediate complex. In the last line of reactions, the phosphoryl group is removed from OmpR-Pthrough the action of EnvZ-D. The concentration of ATP and ADP is assumed to be maintained constant in time. In order to present the theory we develop, we first need to introduce some terminology. The linear combinations ofchemical species appearing on either side of the chemical reactions of interest are called complexes , in accordancewith the Reaction Network Theory literature. Be aware that the word “complex” has usually a different meaningin the Biology literature. We denote by m the number of complexes present in the network, an by d the number ofchemical species. As an example, the complex of (1) are A + B , 2 B , B , and A . Here, d = 2 and m = 4. In (3) thecomplexes are A + B , 2 B , 0, B , A , C , and C + B , hence d = 3 and m = 7. Finally, in the system depicted in Figure 1 d = 8 and m = 10. Since a complex is a linear combination of species, each complex can be regarded as a vector oflength d . For example, for the model (1) we can consider A + B as (1 , B as (0 , B as (0 , A as(1 , stoichiometric subspace as S = span R { y j − y i : there is a reaction from y i to y j } , y n denotes the n th complex, for all 1 ≤ n ≤ m . For example, for (1) we have S = span R (cid:26)(cid:18) − (cid:19) , (cid:18) − (cid:19)(cid:27) = span R (cid:26)(cid:18) − (cid:19)(cid:27) . For (3), we have S = R .In the most general formulation of reaction systems, a (time-dependent) rate function λ ij is associated with thereaction from the i th to the j th complex of the network, and the concentration vector of the different chemical speciesis assumed to solve the differential equation ddt x ( t ) = (cid:88) ≤ i,j ≤ m ( y j − y i ) λ ij ( x ( t ) , t ) , (4)where if a reaction from the i th to the j th complex does not exist, then λ ij is the zero function. Note that (4) simplysums the contributions to the dynamics given by the different chemical reactions. It is also not difficult to show thatevery solution to (4) is necessarily confined within a translation of the stoichiometric subspace. If for all non-zeropropensities λ ij there exists a positive constant κ ij such that λ ij ( x ( t ) , t ) = κ ij d (cid:89) l =1 x l ( t ) y il , then the model is a mass-action system . In this case, (4) can be written as ddt x ( t ) = Y A ( κ )Λ( x ( t )) , where Y is a d × m matrix whose i th column is y i , A ( κ ) is a m × m matrix given by A ( κ ) ij = (cid:40) κ ji if i (cid:54) = j − (cid:80) ml =1 κ il if i = j and Λ( x ( t )) is a vector of length m whose i th entry is (cid:81) dl =1 x l ( t ) y il . Examples of mass-action systems are (1) andthe model in Figure 1.A directed graph can be associated with a reaction network, where the nodes are given by the complexes and thedirected edges are given by the reactions. Such a graph is called reaction graph . As an example, (1) is a reactiongraph, while (3) is not because the complex 0 is repeated. The reaction graph corresponding to (3) is A + B κ BC + B κ u ( t ) B κ Au ( t ) κ C (5)We denote by (cid:96) the number of connected components of the reaction graph associated with the network. For boththe networks (1) and (3) (cid:96) = 2, as well as for the EnvZ-OmpR osmoregulatory system of Figure 1. Then, we definethe deficiency of a network as δ = m − (cid:96) − dim S . The deficiency of a network has important geometric interpretation, and a collection of classical deficiency theoryresults can be found in [21]. The deficiency of (1) is δ = 4 − − δ = 7 − − y is terminal if for all paths in the reaction graph leading from y to anothercomplex y (cid:48) , there is a path leading from y (cid:48) to y . If a complex is not terminal, then it is called non-terminal . As anexample, the only terminal complexes for (1) and (3) are 2 B and A .We recall that a species is said to be absolute concentration robust (ACR) if its concentration at any positivesteady state of (4) is the same. We are ready to state the following result, as presented in [40]. Theorem 3.1.
Consider a mass-action system, and assume the following holds:1. there are two non-terminal complexes y i and y j such that only one entry of y j − y i is non-zero; . the deficiency is 1.3. a positive steady state exists.Then, the species relative to the non-zero entry of y j − y i is ACR. Note that a stronger version of Theorem 3.1 is proven in [40], which detects steady state invariant that are moregeneral than the equilibrium concentration of a single species. The stronger version is stated in the SupplementaryMaterial as Theorem B.3, and an extension of it is proven in the present work.The model (1) has deficiency 1, as already observed, has at least one positive equilibrium and the non-terminalcomplexes A + B and B differ only for the species A . Hence, Theorem 3.1 applies and A is ACR. It is shown in[40] that the EnvZ-OmpR osmoregulatory system in Figure 1 also fulfils the hypothesis of Theorem 3.1, with thenon-terminal complexes EnvZ-D and EnvZ-D+OmpR-P only differing for the species OmpR-P. As a consequence,OmpR-P is ACR. In Section 5 we will develop a method to explicitly calculate the ACR value through symboliclinear algebra. We note that Theorem 3.1 cannot be applied to (3) for two reasons: the model is not a mass-actionsystem and its deficiency is 2.As noted in [5], the positive steady states of a system with an ACR species are not necessarily stable. However, asa consequence of the present work (more precisely, as a consequence of Theorem 6.1 with u being the zero function),we know the following: if a mass-action system as in Theorem 3.1 has an unstable positive steady state, then either thesystem oscillates around it, or some chemical species is completely consumed, or some chemical species is indefinitelyproduced. We give here the formal definition of “oscillation”, as intended in this paper. Definition 3.1.
We say that a function g : R ≥ → R oscillates around a value q ∈ R if for each t ∈ R ≥ there exist t + > t and t − > t such that g ( t + ) > q and g ( t − ) < q. In Control Theory, the focus is usually on systems of differential equations of the form ddt x ( t ) = f ( x ( t ) , u ( t )) , (6)where x : R ≥ → R n x and u : R ≥ → R n u for some n x , n u ∈ Z > , and f is a differentiable function. The function u iscalled the input of the system . Further, a quantity of the form z ( t ) = a ( x ( t )) is of interest, where a is a differentiablefunction with a : R n x → R n z , for some n z ∈ Z > . The function z is called the output of the system. In the usualsetting, one needs to find an appropriate function u such that z is close to a desired level z ∈ R n z , either on averageor for t → ∞ . To this aim, the existence of a function φ : R n x → R such that ddt φ ( x ( t )) = z ( t ) − z is of high importance, an is called an integrator . The name derives from φ ( x ( t )) being the integral of the error thatneeds to be controlled: φ ( x ( t )) = φ ( x (0)) + (cid:90) t ( z ( s ) − z ) ds. If the function is fed back to the system and is used to tune the input, then an integral action or integral feedback is in place [8, 19]. One of the main features of an integrator is that the derivative of φ ( x ( t )) is zero if and only if z ( t ) = z . If a function ˜ φ : R n x → R satisfies ddt ˜ φ ( x ( t )) = r ( x ( t )) (cid:16) z ( t ) − z (cid:17) for some differentiable function r : R n x → R , then ˜ φ is called a constrained integrator (CI) [46]. The name derivesfrom the fact that the derivative of ˜ φ ( x ( t )) is zero if and only if z ( t ) = z , provided that r ( x ( t )) (cid:54) = 0. In biology, it iscommon to find CIs, and the condition r ( x ( t )) (cid:54) = 0 is usually implied by x ( t ) (cid:54) = 0 [46]. Note that in [46] an explicitdistinction between integrators and integral feedbacks is not made.In the setting of systems with ACR species, the output z can be considered to be the concentration of the ACRspecies over time, and z can be their ACR values. In (1), z ( t ) = x A ( t ) and z = κ /κ . A CI (as noted in [46]) isgiven by ˜ φ ( x ( t )) = x B ( t ), since ddt x B ( t ) = κ x B ( t ) (cid:18) x A ( t ) − κ κ (cid:19) . x A ,
0) is a steady state. If an integrator φ existed, then by choosing x (0) = ( x A ,
0) we would have0 = ddt φ ( x ( t )) = x A − κ κ , which cannot hold expect for a specific value of x A . An integrator may still exist in a weaker sense, if we restrict itsdomain. For example, in this case the function ˆ φ ( x ( t )) = κ log x B ( t ) would be an integrator, in the sense that if x B ( t ) > ddt ˆ φ ( x ( t )) = 1 κ x B ( t ) (cid:18) κ x B ( t ) (cid:18) x A ( t ) − κ κ (cid:19)(cid:19) = x A ( t ) − κ κ However, the domain of ˆ φ is not the entire R . Finally, since linear functions could always be extended continuouslyto the boundaries of R > , a linear integrator cannot exist for (1) even if its domain is restricted. We state here our result concerning linear CIs. A stronger version is proved in the Supplementary Material. Theresult is inspired by the analysis carried on in [40], which is here expanded.For any n × l real matrix M and real vector v of length n , we denote by ( M | v ) the n × l + 1 matrix obtained byadding the column v at the right of the matrix A . Let 1 ≤ i, j ≤ m . To present our result, we first need to define thepreimage Γ ij ( κ ) = (cid:8) γ ∈ R d +1 : (cid:0) A ( κ ) (cid:62) Y (cid:62) | e i (cid:1) γ = e j (cid:9) , (7)where e n denotes the n th vector in the canonical basis of R m , whose n th component is 1 and whose other componentsare 0. The role of Γ ij ( κ ) is that of providing vectors γ satisfying ddt (cid:104) ˆ γ, x ( t ) (cid:105) = ˆ γ (cid:62) Y A ( κ )Λ( x ( t )) = Λ j ( x ( t )) − γ i Λ i ( x ( t )) , (8)where ˆ γ is the projection onto the first d components of γ , and (cid:104)· , ·(cid:105) is the standard scalar product. Under certainassumptions, (8) will provide us with a CI. Then, the projection of Γ ij ( κ ) onto the first d coordinates will be ofinterest, and we will denote it by ˆΓ ij ( κ ). The set Γ ij ( κ ) can be calculated with symbolic linear algebra. We willalso prove in the Supplementary Material how ˆΓ ij ( κ ) is connected with S ⊥ , which is a set easily described by linearalgebra and independent on the rate functions. Specifically, we will prove that if ξ ∈ ˆΓ ij ( κ ), then necessarily { ξ + w : w ∈ S ⊥ } ⊆ ˆΓ ij ( κ ) . (9)We will also give sufficient conditions under which the inclusion in (9) is an equality.As an example, consider the model in Figure 1. Using Matlab, we quickly obtain that a vector ξ is in Γ ( κ ), with ξ = κ κ κ ( κ + κ )[ATP] κ ( κ + κ ) κ κ [ADP] , (10)and it is shown in the Supplementary Material thatˆΓ ( κ ) = ˆ ξ + w w w w w w + w w w w : w , w ∈ R , (11)where ˆ ξ is the projection of ξ onto its first d = 8 coordinates.The family of models we study in this paper concerns reaction systems with two non-terminal complexes y i and y j differing in just one entry, for which ˆΓ ij ( κ ) is non-empty. Our first result shows that such a family includes themodels studied in [40]. The proof can be found in the Supplementary Material. Theorem 4.1.
Consider a mass-action system, and assume the following holds:1. there are two non-terminal complexes y i and y j such that only one entry of y j − y i is non-zero; . the deficiency is 1.3. a positive steady state exists.Then, ˆΓ ij ( κ ) is non-empty. As an example, we know already from direct calculation that for the EnvZ-OmpR signaling system (11) holds,which in turn implies that the set ˆΓ ( κ ) is non-empty. However, we could have also derived this information fromTheorem 4.1, without explicitly calculating ˆΓ ( κ ).We note here that the converse of Theorem 4.1 does not hold. We show this with an example of a multisitephosphorylation signaling system in the Supplementary Material, which does not fall in the setting of [40] but forwhich we are able to prove absolute concentration robustness regardless the choice of rate constants, as long as apositive steady state exists. Notably, we are also able to derive information on when this occurs without workingdirectly with the differential equation. As a consequence of this example, the family of models we analyze is provento be strictly larger than that studied in [40]. The following holds. Theorem 4.2.
Consider a mass-action system. Assume that there are two complexes y i and y j only differing in the n th entry, and that ˆΓ ij ( κ ) is non-empty. Let γ ∈ ˆΓ ij ( κ ) , and define q = γ yj − yi ) n d +1 . Then, either no positive steady state exists or the n th species is ACR with ACR value q . Moreover, φ ( x ) = d (cid:88) i =1 β i x i is a linear CI with ddt φ ( x ( t )) = Λ i ( x ( t )) (cid:16) x n ( t ) ( y j − y i ) n − q ( y j − y i ) n (cid:17) for any initial condition x (0) if and only if β ∈ ˆΓ ij ( κ ) . The existence of a linear CI given by Theorem 4.2 is essential to develop the results presented in the next sections.Before unveiling the consequences of Theorem 4.2, however, it is important to stress that a CI does not necessarilyconstitute a feedback, as one may be tempted to think. Consider A + B κ A κ κ B κ B with κ κ = κ κ . It can be shown that the system satisfies the conditions of Theorems 3.1 and 4.1, with the non-terminal complexes A + B and A differing only in species A . Hence, A is ACR and the assumptions of Theorem 4.2hold. A linear CI as in Theorem 4.2 is given by φ ( x ) = − x B /κ , since for this choice ddt φ ( x ( t )) = x B ( t ) (cid:18) x A ( t ) − κ κ (cid:19) . However, the quantity φ ( x ( t )) does not regulate the dynamics of A , since ddt x A ( t ) = κ − κ x A ( t )does not depend on x B ( t ). Since in this case the CI is not acting on the system, it is not surprising that the existenceof positive steady states is lost as soon as κ κ (cid:54) = κ κ .It is also worth mentioning that not all systems with ACR species have a linear CI: consider the mass-actionsystem A + B κ B + C κ κ BB κ κ E κ DC κ AD κ E The model is considered in [3], where it is proven that the species A is ACR. We show in the Supplementary Materialthat there exists no linear function φ whose derivative at time t is of the form r ( x ( t ))( x A ( t ) γ − q ), for some polynomial r and some real numbers γ, q . Note that in this case Theorem 4.1 does not apply because the deficiency of the networkis 2. 7 A method to calculate the ACR value
The first interesting consequence of Theorem 4.2 is that the ACR values of the mass-action systems satisfying theassumption of the theorem can be calculated by finding at least one element of the preimage Γ ij ( κ ), and this canbe done via a simple symbolic linear algebra calculation. As an example, consider the EnvZ-OmpR osmoregulatorysystem in Figure 1. Then, Theorem 4.2 implies that the ACR value of OmpR-P is the value given in (10). Thisvalue is in accordance with the one found in the Supplementary Material of [40], however we found it by calculatinga single element in the preimage of a matrix, as opposed to working with the rather complicated differential equationassociated with the model. An even more involved examples is dealt with in the Supplementary Material. We state an important consequence of Theorem 4.2, a stronger version of which is proven in the SupplementaryMaterial:
Theorem 6.1.
Consider a mass-action system, with associated differential equation ddt x ( t ) = f ( x ( t )) . Assume that there are two complexes y i and y j only differing in the n th entry, and that ˆΓ ij ( κ ) is non-empty. Let q be the ACR value of the n th species. Consider an arbitrary function u with image in R d such that a solution to ddt ˜ x ( t ) = f (˜ x ( t )) + u (˜ x ( t ) , t ) exists. Assume that there exists a ˆ γ ∈ ˆΓ ij ( κ ) which is orthogonal to the vector u ( x, t ) for any x, t . Then, for anyinitial condition ˜ x (0) , at least one of the following holds:(a) the concentration of some species goes to 0 or infinity, along a sequence of times;(b) ˜ x n ( t ) oscillates around q and ˆ γ k (cid:54) = 0 for some k (cid:54) = n ;(c) the integral (cid:90) ∞ t | ˜ x n ( s ) − q | ds tends to , as t goes to infinity. The result implies that if a disturbance orthogonal to a vector ˆ γ ∈ ˆΓ ij ( κ ) is applied over time, then the stabilityof the ACR species is maintained: at most, the ACR species can be forced to oscillate around its original ACRvalue, but it cannot be forced to attain another equilibrium level without causing extinction or overexpression of thechemical species present. We analyze the power of Theorem 6.1 by showing some examples of applications. Example . Consider the mass-action system (1), which fulfills the assumptions of Theorem 6.1 as already observed.Assume the complexes are ordered as A + B , 2 B , B , and A , and the species are ordered alphabetically as A , B .Hence, the two non-terminal complexes differing in the ACR species A are the 1st and the 3rd, and it is shown inthe Supplementary Material that ˆΓ ( κ ) = (cid:26)(cid:18) (cid:19) + (cid:18) ww (cid:19) : w ∈ R (cid:27) . (12)Hence, by choosing w = −
1, we have that (cid:18) − (cid:19) ∈ ˆΓ ( κ ) . (13)This vector is clearly orthogonal to any disturbance acting on the production and degradation rates of the species B . Hence, it follows that the stability and the ACR value of the species A is maintained in (3), provided that nospecies is completely removed or indefinitely expressed. Specifically, since the entry of (13) relative to B is zero, itfollows from Theorem 6.1 that if all the species concentrations are bounded from below and from above by positivequantities, necessarily the concentration of the species A converges to its ACR value as t goes to infinity, despite thedisturbances. 8 xample . Consider the osmoregulatory system in Figure 1, whose featureshave already been discussed in the paper. In particular, we know the species OmpR-P is ACR with ACR value (10).Recall that we ordered the complexes such that the two non-terminal ones differing in OmpR-P are the 1st and the8th. It follows from (11) that for any chemical species, there is a vector in ˆΓ ( κ ) with the associated entry equal to0. It follows that even if the production and degradation of any chemical species in the model is tampered with, thestability and the ACR value of the species OmpR-P are maintained, in the sense described by Theorem 6.1.We can push the disturbances further. By appropriately choosing w and w in (11), we can see that there isvector in ˆΓ ( κ ) whose entries relative to the species OmpR and EnvZ are both 0. Hence, it follows by Theorem 6.1that by tampering with the production and degradation rates of both these species over time, if no extinction andno overexpression occurs, then the concentration of OmpR-P still converges to the value (10), or oscillates around it.As a final remark, we note that (9) can be useful in determining whether a vector in Γ ij ( κ ) exists, with a specificcomponent equal to 0, say the n th one. In fact, the existence of such a vector can be deduced without calculations,if there is a vector w ∈ S ⊥ whose n th component is different from 0. Here we illustrate how the theory we developed can be utilized to solve the Synthetic Biology problem of retroactivity.As explained in the Introduction, the loading effects caused by a downstream biochemical module can disrupt thefunctionality of upstream modules, which prevents the implementation of biochemical circuits by interconnectingbiochemical modules with different functions [17]. A concrete example of loading effect is illustrated in Figure 4.Assume that a mass-action system has two complexes y i and y j , that are only different in the n th component,which corresponds to the species X . Assume further that ˆΓ ij ( κ ) is non-empty and that a positive steady state exists.Hence, the species X is ACR, with some ACR value q . It further follows from Theorem 6.1 that, if there existsˆ γ ∈ ˆΓ ij ( κ ) with ˆ γ n = 0, then the production and degradation rates of the species X can be arbitrarily perturbedover time by an arbitrary function u , without compromising its robustness. Specifically, if the perturbed system isstable and no chemical species is completely consumed, then the concentration of X will still converge to the sameACR value q as in the original mass-action system. The key observation we make here is that the perturbation u can be considered as the loading effect of a downstream module that takes the concentration of species X as input.In this case, the loading effect on the original mass-action is rejected and the concentration of X is maintained at adesired level q at steady state. Further, the concentration of X is maintained approximately constant in the transientdynamics as well, if we assume as done in [17] that a separation of dynamics time scale is in place. Specifically,assume ddt ˜ x ( t ) = 1 ε f (˜ x ( t )) + u (˜ x ( t ) , t ) e n , for some small ε >
0, with f (˜ x ( t )) and u (˜ x ( t ) , t ) being of the same order of magnitude. Under the assumption ofstability, if ε is very small then the perturbed system will quickly approach the slow manifold defined by0 = f (˜ x ( t )) + εu t (˜ x ( t )) e n , where u t is a function from R d to R d defined by u t ( x ) = u ( x, t ). By Theorem 6.1 applied to the disturbance u t , thespecies X assumes its ACR value at any positive point of the slow manifold, which is exactly what we wanted.As an example of application, consider the EnvZ-OmpR osmoregulatory system in Figure 1. It follows from(11) that there exists ˆ γ ∈ Γ ( κ ) with a zero in the entry corresponding to the ACR species OmpR-P. Hence, theproduction and degradation rates of OmpR-P can be arbitrarily changed over time, without altering its robustnessproperty, in the sense described by Theorem 6.1. As observed, the statement still holds true if the perturbation isoriginated by a downstream module that acts on OmpR-P. Hence, the EnvZ-OmpR osmoregulatory system can beused to maintain the expression of OmpR-P at a desired level, which depends on the input rate constants, even ifthe species OmpR-P is used by a downstream module. Moreover, if the downstream module acts on a slower timescale, the concentration of the species OmpR-P is approximately maintained at the target level at any time point. InFigure 2, a diagram describing this situation is proposed.Consider now the case were an upstream module is affected from loading effects. We show how the theorydeveloped in this paper can be used to design an insulator. Assume that the upstream module accepts u ( t ) as input,and modulates the concentration of the species A (cid:63) accordingly. The species A (cid:63) is then used by a downstream module,which returns a function of the concentration of A (cid:63) as output. The action of the downstream module on the species A (cid:63) causes a loading effect on the upstream module. This loading effect should be reduced. To this aim, we proposedto modify the downstream module such that it acts on a species A rather than on the species A (cid:63) , and to include inthe system the following module, where B is a species that is not used by neither the upstream nor the downstreammodule: A + B κ BA (cid:63) + B κ A (cid:63) + A (14)9nvZ-OmpR systemDownstream module κ OmpR-POmpR-P g (OmpR-P)Figure 2: Proposed use of the EnvZ-OmpR signal transduction system of Figure 1 as a controller of a downstreammodule utilizing OmpR-P. The concentration of OmpR-P is regulated by the EnvZ-OmpR signaling system, withequilibrium given by (10). The equilibrium can be adjusted by modifying the parameters κ ij of the EnvZ-OmpRsignaling system (which depend on the medium osmolarity) or the ratio between ADP and ATP present. Theoutput of the downstream module is a function g of the concentration of OmpR-P, which is received as input.Upstream moduleInsulator A + B κ BA (cid:63) + B κ A (cid:63) + A Downstream module u ( t ) A (cid:63) (not changedby the insulator) κ , κ , x B (0) A ≈ κ A (cid:63) κ A (changed by down-stream module) g ( A )Figure 3: The upstream module expresses the chemical species A (cid:63) as output. The insulator transfers a multipleof the signal from the upstream module to the downstream module, which is modified to accept as input theconcentration of A rather than the concentration of A (cid:63) .Assume stability is reached and that the species B is not completely consumed. Then, at steady state the concen-tration level of A (cid:63) is fixed, and the concentration of the ACR species A will converge to its ACR value κ x A (cid:63) /κ regardless any disturbance applied to the production and degradation rate of A . In fact, a linear CI as in the statementof Theorem 4.2 is given by φ ( x ) = x B /κ , and at any time point ddt φ ( x ( t )) = x B ( t ) (cid:18) x A ( t ) − κ κ x A (cid:63) ( t ) (cid:19) . (15)We further note that if the dynamics of (14) occur on a faster time scale than the rest of the system, then a slowmanifold is quickly approached were the concentration of the species A is maintained at the level κ x A (cid:63) ( t ) /κ at anytime point. In this case, the module (14) approximately outputs a multiple of the concentration of A (cid:63) over the wholetime line. The multiplicative constant can be tuned through the parameters κ and κ , as well as the time scale that(14) operates in. The time scale can be further tuned via the concentration of B , as it also follows from (15). Inconclusion, the downstream module receives as input a good approximation of a multiple of the concentration of A (cid:63) ,and its activity does not affect the upstream module, nor (14). Moreover, (14) does not affect the upstream moduleat all, since the species A appears in (14) as a catalyst and is not changed in the catalysed reaction. The proposedinsulating strategy is illustrated in Figure 3, and it is applied to an example discussed in [17] in Figure 4.10solated systems (with x A (cid:63) ( t ) = x A (cid:63) ( t ))0 u ( t )0 . A (cid:63) A (cid:63) + P Cu ( t ) x A (cid:63) ( t ) x A (cid:63) ( t ) x C ( t )Coupled systems0 u ( t )0 . A (cid:63) A (cid:63) + P Cu ( t ) x A (cid:63) ( t ) x C ( t )Insulated system0 u ( t )0 . A (cid:63) A + B . BA (cid:63) + B . A (cid:63) + AA + P Cu ( t ) x A (cid:63) ( t ) x A ( t ) x C ( t )Figure 4: In the example, we let u ( t ) = 0 . . t )) and x A (cid:63) (0) = 5. All the rates are in 1 /min . The ODEsolution is calculated in Matlab with ode23s . In the first panel, the two systems are considered in isolation, with theassumption that the concentration of the species A (cid:63) is maintained at the same level as the concentration of species A (cid:63) , at any time point. In the second panel, the two systems are linked together, and the species A (cid:63) is directly usedby the downstream system. The dynamics are completely disrupted by the loading effects, a plot of the absolutevalue of the difference of the two solutions over time is proposed. In the third panel, the insulator of Figure 3 isutilized, with x A (0) = 0 and x B (0) = 20. The loading effect are practically removed, despite the choice of low rateconstants 0 .
1. Specifically, the difference of the concentration of C between the solution of insulated system and thesolution of the isolated systems spikes quickly to 40, but it decreases to less than 0.5 within 10 minutes, after whichis maintained low as illustrated in the second plot of the third panel. In the previous section, we have seen how the stability properties of an ACR system can be transferred, in a sense,to a larger model including further chemical transformations and external inputs. As a particular case, Here, we11xplicitly state when the property of absolute concentration robustness can be lifted from portions of the biochemicalsystem to the whole system. As a consequence, we further extend the set of sufficient conditions of Theorem 3.1 thatimply the existence of an ACR species. Before stating the relevant result, which is a consequence of Theorem 4.2, weneed a definition. Given a reaction system S , we say that S (cid:48) is a sub-system of S if it can be obtained from S by canceling some reactions, and if the choice of rate functions for the remaining reactions is maintained. Moreover,if S and S (cid:48) have d and d (cid:48) species, respectively, we let π : R d → R d (cid:48) be the projection onto the species of S (cid:48) . Thefollowing holds: Corollary 6.2.
Consider a reaction system S , and let S (cid:48) be a sub-system. Assume that S (cid:48) is a mass-actionsystem with two complexes π ( y i ) and π ( y j ) only differing in the entry relative to the species X , and for which ˆΓ ij ( κ ) is non-empty. Moreover, assume there exists ˆ γ ∈ ˆΓ ij ( κ ) such that π ( y l − y k ) is orthogonal to ˆ γ for all y k → y l thatare reactions of S but not reactions of S (cid:48) . Then, the X is ACR for both S and S (cid:48) , with the same ACR value. The proof of a stronger result is in the Supplementary Material. Here we illustrate how the corollary can beapplied in the case of EnvZ-OmpR signaling system: consider the reaction system S described in Figure 5, whichincludes the EnvZ-OmpR osmoregulation system described in Figure 1 as a sub-system. We assume that a protein canmisfold when the phosphoryl group is transferred from EnvZ to OmpR. Such misfold can be corrected by chaperons,which are independently produced and degraded through a mechanism which we assume unkown, but that does notinvolve EnvZ or OmpR proteins. We also allow for an arbitrary and persistent external control on the expressionlevel of EnvZ sensor-transmitter protein. Finally, we consider the utilization of OmpR-P as transcription regulatoryprotein of the outer membrane porins OmpF and OmpC. For our purposes, we assume the details of the transcriptionmechanism are not known, but that only the protein Ompr-P is involved in the process. As previously done, let thecomplexes of the EnvZ-OmpR osmoregulation system be ordered from left to right and from top to bottom, suchthat EnvZ-D and EnvZ-D+OmpR-P are the 1st and the 8th complex, respectively. Also, let the species be orderedaccording to their appearance from left to right and from top to bottom, as EnvZ-D, EnvZ, EnvZ-T, EnvZ-P, OmpR,EnvZ-OmpR-P, OmpR-P, EnvZ-OmpR-D-P. In particular, Envz is the second species, EnvZ-OmpR-P is the sixthspecies, and OmprR-P is the seventh species. It follows from (11) that, by choosing w = − ξ and w = − ξ = 0, avector ˆ ζ is in ˆΓ ( κ ) with:1. ˆ ζ = ˆ ξ + w = 0;2. ˆ ζ − ˆ ζ = ˆ ξ + w + w − ˆ ξ − w = ˆ ξ − ˆ ξ = 0;3. ˆ ζ = ˆ ξ + w = 0.Denote by E k the vector of R d with the k th entry equal to 1 and the other entries equal to zero. The following holds. Misfolding of OmpR-P.
The projection of the difference between EnvZ + OmpR (cid:63) -P and EnvZ-OmpR-P onto thespecies of the EnvZ-OmpR signaling system is E − E , which is orthogonal to ˆ ζ by 2. The projection of thedifference between OmpR-P + C and OmpR (cid:63) -P + C is E , which is orthogonal to ˆ ζ by 3. Production and degradation of chaperons.
By assumption, any chemical reaction y → y (cid:48) involved in the productionand degradation of chaperones does not consume or produce any chemical species of the EnvZ-OmpR signalingsystem. Hence, π ( y (cid:48) − y ) = 0, which is orthogonal to ˆ ζ . External regulation of EnvZ.
The difference between EnvZ and 0 is ± E , which is orthogonal to ˆ ζ by 1 Transcription of OmpC and OmpF.
We assume that the transcription only involves the species OmpR-P, out of allthe species in the EnvZ-OmpR osmoregulation system. Hence, for all the reactions y → y (cid:48) involved in thetranscription, either π ( y (cid:48) − y ) = 0 or π ( y (cid:48) − y ) is a multiple of E . In either case, π ( y (cid:48) − y ) is orthogonal to ˆ ζ .It follows from Corollary 6.2 that the species OmpR-P is ACR in the reaction system of Figure 5. Moreover, itsACR value is still given by (10), as long as a positive steady state exists. Note that Corollary 6.2 could be appliedeven if not all chemical reactions are known, and even if the model is not mass-action. It is also worth noting thatthe deficiency of the model is not known, due to the lack of information on the precise reactions constituting thenetwork, but is certainly greater than 1. Indeed, the deficiency of the sub-system constituted by the EnvZ-OmpRosmoregulation system and by the misfolding of OmpR-P is 2, and the deficiency of a system is necessarily greaterthan or equal to the deficiency of any sub-system [15, Lemma 5]. We have shown in Theorem 4.2 that a linear CI always exists for a family of ACR systems, and that this familystrictly includes the models studied in [40]. The result, a more general version of which is proven in the Supplementary12nvZ-OmpR osmoregulation system, as described in Figure 1Misfolding of sensory response protein OmpR- P , and recovery through the action of chaperones:EnvZ-OmpR- P κ EnvZ + OmpR (cid:63) - P OmpR (cid:63) - P + C κ OmpR- P + CProduction and degradation of chaperones:Gene regulatory network C 0 u C ( t )External regulation of the sensor-transmitter protein EnvZ:0 u ( t ) u ( t ) EnvZTranscription of OmpC and OmpF:OmpR- P binds to DNA promoterOmpCOmpFFigure 5: Reaction system including the EnvZ-OmpR osmoregulation system depicted in Figure 1 as a sub-system. Parts of the system depicted here are unknown, specifically no model is given for the production ofchaperons or for the transcription of the outer membrane porins OmpF and OmpC.Material, has three main consequences: first, it provides an easy way to calculate the ACR value of ACR species.Secondly, as expressed in Theorem 6.1, the presence of a linear CI implies that the system is robust to arbitrarydisturbances that do not vanish over time, under certain conditions. This fact can be naturally exploited to designperfect insulators, which are able to reject loading effects originated from the downstream modules. Finally, asexpressed in Corollary 6.2, we are able to prove that, under certain conditions, the absolute concentration robustnessof a portion of a system can be lifted to the whole model, and the ACR value of the ACR species remains unchanged.The theory we developed opens the path for future research directions. First, efficient algorithms can be designedin order to check for the existence of portions of the systems that confer absolute concentration robustness to thewhole system. To this aim, the connections of ˆΓ ij ( κ ) with structural properties of the network which we show inthe Supplementary Material can be useful, and theoretical results can be expanded in this direction. Second, furtherdetailed analysis on when stability can be ensured would be welcome. Currently, in the statement of Theorem 6.1 wecannot exclude the possibility that some species is completely consumed or indefinitely produced upon tampering withthe model. Finding structural conditions able to eliminate this possibility would be a nice and useful contribution.As a final remark, we think the study of stochastically modeled systems that satisfy the assumptions of Theorem 4.2would be interesting and fruitful. Stochastic models of reaction systems are tipycally used when few molecules ofcertain chemical species are available [6, 20]. It is proven in [5] that systems satisfying the assumptions of Theorem 3.1,when stochastically modeled, undergo an extinction event almost surely. As a consequence, the desirable robustnessproperties of the ACR systems studied in [40] are completely destroyed in a low molecule copy-number regime. As anexample, the model depicted in (1) undergoes an almost sure extinction of the chemical species B when stochasticallymodeled, regardless the initial conditions. This is caused by the fact that all the molecules of B can be consumed bythe reaction B → A , before the occurrence of a reaction A + B → B . Robustness at finite time intervals of somestochastically modeled ACR systems is recovered, but only in a multiscale limit sense [4]. Moreover, it is shown in[3] that absolute concentration robustness of the deterministic model does not necessarily imply an extinction eventin the corresponding stochastic model, but the connection is still largely unexplored. The results developed in thepresent paper can help in this direction: consider again (1). The extinction of species B cannot occur if production13f B is included in the model as in (3), or as in A + B κ BB κ A κ κ B (16)At the same time, it follows from Theorem 6.1 that the stability properties of the species A are maintained both in(3) and in (16), when deterministically modeled. In particular, the concentration of the species A still converges tothe value κ /κ . It would be interesting to study if in this and in similar cases some form of absolute concentrationrobustness arise in the long-term dynamics of the stochastic models as well. Acknowledgement
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon2020 research and innovation programme grant agreement no. 743269 (CyberGenetics project).
References [1] U Alon, MG Surette, N Barkai, and S Leibler. Robustness in bacterial chemotaxis.
Nature , 397(6715):168, 1999.[2] DF Anderson. A modified next reaction method for simulating chemical systems with time dependent propensitiesand delays.
J Chem Phys , 127(21):214107, 2007.[3] DF Anderson and D Cappelletti. Discrepancies between extinction events and boundary equilibria in reactionnetworks. arXiv preprint arXiv:1809.04613 , 2018.[4] DF Anderson, D Cappelletti, and TG Kurtz. Finite time distributions of stochastically modeled chemical systemswith absolute concentration robustness.
SIAM J Appl Dyn Syst , 16(3):1309–1339, 2017.[5] DF Anderson, GA Enciso, and MD Johnston. Stochastic analysis of biochemical reaction networks with absoluteconcentration robustness.
J R Soc Interface , 11(93):20130943, 2014.[6] DF Anderson and TG Kurtz.
Stochastic analysis of biochemical systems , volume 1. Springer, 2015.[7] D Angeli, JE Ferrell, and ED Sontag. Detection of multistability, bifurcations, and hysteresis in a large class ofbiological positive-feedback systems.
Proc Natl Acad Sci U S A , 101(7):1822–1827, 2004.[8] KJ Astr¨om and RM Murray.
Feedback systems: an introduction for scientists and engineers . Princeton universitypress, 2010.[9] M Banaji and C Pantea. The inheritance of nondegenerate multistationarity in chemical reaction networks.
SIAM J Appl Math , 78(2):1105–1130, 2018.[10] N Barkai and S Leibler. Robustness in simple biochemical networks.
Nature , 387(6636):913, 1997.[11] E Batchelor and M Goulian. Robustness and the cycle of phosphorylation and dephosphorylation in a two-component regulatory system.
Proc Natl Acad Sci U S A , 100(2):691–696, 2003.[12] F Blanchini and E Franco. Structurally robust biological networks.
BMC Syst Biol , 5(1):74, 2011.[13] JD Brunner and G Craciun. Robust persistence and permanence of polynomial and power law dynamical systems.
SIAM J Appl Math , 78(2):801–825, 2018.[14] D Cappelletti, AP Majumder, and C Wiuf. Fixed–time and long–term dynamics of monomolecular reactionnetworks in stochastic environment. in preparation , 2019.[15] D Cappelletti and C Wiuf. Product-form Poisson–like distributions and complex balanced reaction systems.
SIAM J Appl Math , 76(1):411–432, 2016.[16] G Craciun, F Nazarov, and C Pantea. Persistence and permanence of mass-action and power-law dynamicalsystems.
SIAM J Appl Math , 73(1):305–329, 2013. 1417] D Del Vecchio, AJ Ninfa, and ED Sontag. Modular cell biology: retroactivity and insulation.
Mol Syst Biol ,4(1):161, 2008.[18] JP Dexter, T Dasgupta, and J Gunawardena. Invariants reveal multiple forms of robustness in bifunctionalenzyme systems.
Integrative Biology , 7(8):883–894, 2015.[19] JC Doyle, BA Francis, and AR Tannenbaum.
Feedback control theory . Courier Corporation, 2013.[20] P ´Erdi and J T´oth.
Mathematical models of chemical reactions: theory and applications of deterministic andstochastic models . Manchester University Press, 1989.[21] M Feinberg. Chemical reaction network structure and the stability of complex isothermal reactors—i. thedeficiency zero and deficiency one theorems.
Chem Eng Sci , 42(10):2229–2268, 1987.[22] M Feinberg.
Foundations of Chemical Reaction Network Theory . Springer, 2019.[23] M Feinberg and F Horn. Chemical mechanism structure and the coincidence of the stoichiometric and kineticsubspaces.
Arch Ration Mech Anal , 66(1):83–97, 1977.[24] E Feliu, D Cappelletti, and C Wiuf. Node balanced steady states: Unifying and generalizing complex anddetailed balanced steady states.
Math Biosci , 2018.[25] E Feliu and C Wiuf. Enzyme–sharing as a cause of multi–stationarity in signalling systems.
J R Soc Interface ,page 20110664, 2011.[26] S Feng, M S´aez, C Wiuf, E Feliu, and OS Soyer. Core signalling motif displaying multistability through multi-state enzymes.
J R Soc Interface , 13(123):20160524, 2016.[27] E Franco and F Blanchini. Structural properties of the MAPK pathway topologies in PC12 cells.
J Math Biol ,67(6-7):1633–1668, 2013.[28] M Gopalkrishnan, E Miller, and A Shiu. A geometric approach to the global attractor conjecture.
SIAM J ApplDyn Syst , 13(2):758–797, 2014.[29] E Gross, HA Harrington, N Meshkat, and A Shiu. Joining and decomposing reaction networks. arXiv preprintarXiv:1810.05575 , 2018.[30] LH Hartwell, JJ Hopfield, S Leibler, and AW Murray. From molecular to modular cell biology.
Nature ,402(6761supp):C47, 1999.[31] T Jahnke and W Huisinga. Solving the chemical master equation for monomolecular reaction systems analytically.
J Math Biol , 54(1):1–26, 2007.[32] B Joshi and A Shiu. Atoms of multistationarity in chemical reaction networks.
J Math Chem , 51(1):153–178,2013.[33] RL Karp, MP Mill´an, T Dasgupta, A Dickenstein, and J Gunawardena. Complex-linear invariants of biochemicalnetworks.
J Theor Biol , 311:130–138, 2012.[34] D Mishra, PM Rivera, A Lin, D Del Vecchio, and R Weiss. A load driver device for engineering modularity inbiological networks.
Nat Biotechnol , 32(12):1268, 2014.[35] T Miyashiro and M Goulian. High stimulus unmasks positive feedback in an autoregulated bacterial signalingcircuit.
Proc Natl Acad Sci U S A , 105(45):17457–17462, 2008.[36] L Pantoja-Hern´andez and JC Mart´ınez-Garc´ıa. Retroactivity in the context of modularly structured biomolecularsystems.
Front Bioeng Biotechnol , 3:85, 2015.[37] LA Pratt and TJ Silhavy. Porin regulon of escherichia coli. In
Two-Component Signal Transduction , pages105–127. American Society of Microbiology, 1995.[38] PEM Purnick and R Weiss. The second wave of synthetic biology: from modules to systems.
Nat Rev Mol CellBiol , 10(6):410, 2009.[39] G Shinar, U Alon, and M Feinberg. Sensitivity and robustness in chemical reaction networks.
SIAM J ApplMath , 69(4):977–998, 2009. 1540] G Shinar and M Feinberg. Structural sources of robustness in biochemical reaction networks.
Science ,327(5971):1389–1391, 2010.[41] G Shinar, R Milo, MR Mart´ınez, and U Alon. Input–output robustness in simple bacterial signaling systems.
Proc Natl Acad Sci U S A , 104(50):19931–19935, 2007.[42] G Shinar, JD Rabinowitz, and U Alon. Robustness in glyoxylate bypass regulation.
PLoS Comput Biol ,5(3):e1000297, 2009.[43] AM Stock, VL Robinson, and PN Goudreau. Two-component signal transduction.
Annu Rev Biochem , 69(1):183–215, 2000.[44] K Takahashi, S T˘anase-Nicola, and PR Ten Wolde. Spatio–temporal correlations can drastically change theresponse of a mapk pathway.
Proc Natl Acad Sci U S A , 107(6):2473–2478, 2010.[45] J T´oth, AL Nagy, and D Papp.
Reaction kinetics: exercises, programs and theorems . Springer, 2018.[46] F Xiao and JC Doyle. Robust perfect adaptation in biomolecular reaction networks. bioRxiv , page 299057, 2018.[47] TM Yi, Y Huang, MI Simon, and JC Doyle. Robust perfect adaptation in bacterial chemotaxis through integralfeedback control.
Proc Natl Acad Sci U S A , 97(9):4649–4653, 2000.16 upplementary MaterialA Notation
A.1 General notation
We will denote by R , R > , and R ≥ the real, positive real, and non-negative real numbers, respectively. Similarly, wewill denote by Z , Z > , and Z ≥ the integer, positive integer, and non-negative integer numbers, respectively. Givena real number a , we will denote its absolute value by | a | .We denote by e i the vector that has 1 in its i th entry and 0 in all other entries. The dimension of such vector willbe clear from the context, and if not it will be made explicit. Given two real vectors v , w of the same length n , wedenote their scalar product by (cid:104) v, w (cid:105) , and we use the shorthand notation v w = n (cid:89) i =1 v w i i , where 0 is considered to be 1. We will denote the euclidean norm of v by (cid:107) v (cid:107) .Given two subsets V, W ⊆ R n , we denote their direct sum by V ⊕ W = { v + w : v ∈ V, w ∈ W } . Moreover, given a vector v ∈ R n , we define v + W = { v + w : w ∈ W } . A.2 Reaction network terminology
A.2.1 Standard notation
A reaction network G is a triple {X , C , R} , where • X is a finite ordered set of d different symbols, called species ; • C is a finite ordered set of m linear combinations of species with non-negative integer coefficients, referred to as complexes and identified with vectors in Z d ≥ ; • R is a finite ordered set of elements of C × C , referred to as reactions .We will denote by X i the i th species of X , and by y j the j th complex of C , for 1 ≤ i ≤ d and 1 ≤ j ≤ m . A reaction( y i , y j ) will be denoted by y i → y j , and for any y i ∈ C we assume there is no reaction of the form y i → y i .As mentioned, a complex y i can be regarded as a vector of Z d ≥ . Specifically, this is done by considering the j thentry y ij as the coefficient of y i relative to the j th species. In this regards, for any vector v ∈ R d we will denote bysupp v the subset of species such that X i ∈ supp v if and only if v i (cid:54) = 0 . Similarly, for any vector w ∈ R m we will denote by supp w the subset of complexes such that y i ∈ supp v if and only if v i (cid:54) = 0 . We denote by S = span R { y j − y i : y i → y j ∈ R} , S ⊥ = { h ∈ R d : (cid:104) h, y j − y i (cid:105) = 0 for all y i → y j ∈ R} . The subspace S is called stoichiometric subspace , and the elements of S ⊥ are called conservation laws .A choice of kinetics for a reaction network is a set of (time-dependent) rate functions λ ij : R d ≥ × R ≥ → R ≥ forall 1 ≤ i, j ≤ m , such that λ ij is a zero function if and only if y i → y j / ∈ R . A reaction network with a choice ofkinetics S = ( X , C , R , { λ ij } ≤ i,j ≤ m ) is termed reaction system . A reaction system is associated with the system ofdifferential equations ddt x ( t ) = (cid:88) ≤ i,j ≤ m ( y j − y i ) λ ij ( x ( t ) , t ) for all t ∈ R ≥ . We note that rate functions are commonly intended to not depend on time, but we consider this more general settingin this paper. 17 reaction system is called mass-action system , and denoted by S = ( X , C , R , κ ), if there is a matrix κ ∈ R m × m ≥ such that λ ij ( x, t ) = κ ij x y i for all 1 ≤ i, j ≤ m, x ∈ R d ≥ . In this case, the constants κ ij are termed rate constants . Define Λ( x ) ∈ R d ≥ byΛ i ( x ) = x y i for all 1 ≤ i ≤ m, x ∈ R d ≥ , and define the m × m matrix A ( κ ) by A ( κ ) ij = (cid:40) κ ji if i (cid:54) = j − (cid:80) ml =1 κ il if i = j Finally, let Y be the d × m matrix Y with entries Y ij = y ji . Then, for mass-action sytems we have ddt x ( t ) = (cid:88) ≤ i,j ≤ m ( y j − y i ) κ ij x ( t ) y i = Y A ( κ )Λ( x ( t )) for all t ∈ R ≥ . The directed graph {C , R} is called reaction graph . Reaction systems are often presented through the reactiongraph, where indication on the reaction rates are written on top of the arrows that correspond to the related reaction.Specifically, the arrow corresponding to y i → y j is labeled with λ ij ( x, t ) x y , which corresponds to the rate constants for reaction rates of mass-action type. As an example, see (1) and (5) in themain text, or Figure 1. We denote by (cid:96) the number of connected components of the reaction graph. We define the deficiency as the number δ = m − (cid:96) − dim S . A terminal component of the reaction graph is a set of complexes T ⊆ C , such that • if y ∈ T and there is a directed path in the reaction graph from y to another complex y (cid:48) , then y (cid:48) ∈ T ; • for any two complexes y, y (cid:48) ∈ T with y (cid:54) = y (cid:48) , there is a directed path from y to y (cid:48) , and a directed path from y (cid:48) to y .Let τ be the number of terminal components, and denote by T , T , . . . , T τ the different terminal components ofthe network. Since each connected component contains at least one terminal component, we have τ ≥ (cid:96) . We say acomplex is terminal if it is contained in a terminal component, and we say that a complex is non-terminal otherwise.As an example, the terminal components in (1) are { B } and { A } , hence the terminal complexes are 2 B and A . A.3 Absolute concentration robustness
Definition A.1.
Consider a reaction system S = ( X , C , R , { λ ij } ≤ i,j ≤ m ). A species X n ∈ X is said to be absolutelyconcentration robust (ACR) if there exists q ∈ R > such that c n = q for all positive steady state c of S . In this case, q is referred to as the ACR value of X n . Finally, if an ACR species exists, then the reaction system S is said to be ACR . Remark
A.1 . By definition, if a reaction system has no positive steady states, or a unique one, then all the species areACR. We can call this a degenerate case , since, especially if no positive steady states exist, the concept of robustnessis lost.
A.4 Control Theory terminology
Consider a differential equation of the form ddt x ( t ) = f ( x ( t )) + g ( x ( t ) , u ( t )) for all t ∈ R ≥ , where x : R ≥ → R d and u : R ≥ → R n u for some n u ∈ Z > , and f and g are differentiable functions. The function u is called the input of the system . Further, we define output of the system the quantity z ( t ) = a ( x ( t )), for somedifferentiable function a , with a : R d → R n z and n z ∈ Z > . Usually, one needs to find an appropriate function u such18hat z is close to a desired level z ∈ R n z , either on average or for t → ∞ . To this aim, the existence of a function φ : R n x → R such that ddt φ ( x ( t )) = z ( t ) − z is of high importance, an is called an integral feedback (IF) [8, 19]. If a function ˜ φ : R n x → R satisfies ddt ˜ φ ( x ( t )) = r ( x ( t )) (cid:16) z ( t ) − z (cid:17) for some differentiable function r : R n x → R , then ˜ φ is called a constrained integral feedback (CIF) [46].In the setting of ACR species, z ( t ) is usually the concentration of the ACR species at time t , and z is its ACRvalue. B Known results
The following result appears in the Appendix of [23].
Theorem B.1.
All the vectors in ker A ( κ ) have support in the terminal complexes. Specifically, there exists a basis { χ ( κ ) , χ ( κ ) , . . . , χ τ ( κ ) } of ker A ( κ ) such that supp χ i ( κ ) = T i for all ≤ i ≤ m . The following result can be deduced from Section 6 in [23].
Theorem B.2.
We have δ ≥ dim (cid:16) ker Y ∩ Im A ( κ ) (cid:17) = dim ker Y A ( κ ) − dim ker A ( κ ) . Finally, the following result is proven in the supplementary material of [40], where it is stated as Theorem S3.15.
Theorem B.3.
Consider a mass-action system ( X , C , R , κ ) . Assume the following conditions hold:1. y i and y j are non-terminal complexes;2. the deficiency is 1;Then, there exists q ∈ R > such that c y j − y i = q for all positive steady states c . For convenience, we restate here the weaker version given in the main text:
Theorem 3.1.
Consider a mass-action system, and assume the following holds:1. there are two non-terminal complexes y i and y j such that only one entry of y j − y i is non-zero;2. the deficiency is 1.3. a positive steady state exists.Then, the species relative to the non-zero entry of y j − y i is ACR. Note how Theorem 3.1 is a straightforward consequence of Theorem B.3: if the two complexes only differ in the n th entry, then c y j − y i = c ( y j − y i ) n n = q for all positive steady states c , hence X n is ACR. Remark
B.1 . In the original statement of Theorem B.3, the additional hypothesis that a positive steady state exists ismade. In fact, under the condition of Theorem B.3, the existence of positive steady states is not guaranteed. However,if no positive steady state exists then the conclusion of the theorem holds trivially. Hence, our reformulation holdscorrect. See Remark A.1 for the formal relationship between ACR species and the lack of positive steady states.19
Calculation and structural properties of Γ ij ( κ ) Consider a mass-action system S = ( X , C , R , κ ). For any 1 ≤ i, j ≤ m with i (cid:54) = j define the setΓ ij ( κ ) = (cid:8) γ ∈ R d +1 : (cid:0) A ( κ ) (cid:62) Y (cid:62) | e i (cid:1) γ = e j (cid:9) , (C.1)and for any real vector v of length larger than d , let π d ( v ) be its projection onto the first d components. We defineˆΓ ij ( κ ) = (cid:8) ˆ γ ∈ R d : ˆ γ = π d ( γ ) for some γ ∈ Γ ij ( κ ) (cid:9) . (C.2)The set Γ ij ( κ ) can be computed by first calculating a basis forΨ ij ( κ ) = ker (cid:0) A ( κ ) (cid:62) Y (cid:62) | e i | e j (cid:1) . (C.3)Note that this can be easily and quickly done by using a programming language which is able to deal with symboliclinear algebra, such as Matlab. The following holds Proposition C.1.
Consider a mass-action system S = ( X , C , R , κ ) , and ≤ i, j ≤ m with i (cid:54) = j . Let { ψ , ψ , . . . , ψ k } be a basis for Ψ ij ( κ ) , as defined in (C.3) . Then, Γ ij ( κ ) (and consequently ˆΓ ij ( κ ) ) is non-empty if and only if ψ nd +2 (cid:54) = 0 for some ≤ n ≤ k . If this is the case, then Γ ij ( κ ) = (cid:40) − (cid:80) kn =1 a n ψ nd +2 k (cid:88) n =1 a n π d +1 ( ψ n ) : a , a , . . . , a k ∈ R and k (cid:88) n =1 a n ψ nd +2 (cid:54) = 0 (cid:41) , where π d +1 : R d +2 → R d +1 is the projection onto the first d + 1 components, and ˆΓ ij ( κ ) = (cid:40) − (cid:80) kn =1 a n ψ nd +2 k (cid:88) n =1 a n π d ( ψ n ) : a , a , . . . , a k ∈ R and k (cid:88) n =1 a n ψ nd +2 (cid:54) = 0 (cid:41) . Proof.
First, note that γ ∈ Γ ij ( κ ) if and only ( γ, − ∈ Ψ ij ( κ ). This can be easily deduced by the definitions (C.1)and (C.3). The proof simply follows from this equivalence, and from the definition of ˆΓ ij ( κ ) given in (C.2).In Section F we will use Proposition C.1 to calculate ˆΓ ij ( κ ) for the main examples discussed in the main text.We will see that some vectors in the basis of Ψ ij ( κ ) do not depend on the particular choice of rate constants. Asa consequence, some dynamical properties implied by the theory developed in this paper will only depend on thestructure of the model rather than on a fine tuning of the parameters, which is desirable. In the following result, weexplicitly derive structural properties of ˆΓ ij ( κ ). Proposition C.2.
Consider a mass-action system S = ( X , C , R , κ ) , and ≤ i, j ≤ m with i (cid:54) = j . Assume that ˆΓ ij ( κ ) is non-empty, and let ˆ γ ∈ ˆΓ ij ( κ ) . Then, ˆ γ + S ⊥ ⊆ ˆ γ + ker A ( κ ) (cid:62) Y (cid:62) ⊆ ˆΓ ij ( κ ) . (C.4) Moreover, if the mass-action system has a steady state c with c y i > , then γ (cid:48) d +1 = c y j − y i for all γ (cid:48) ∈ Γ ij ( κ ) (C.5) and ˆΓ ij ( κ ) = ˆ γ + ker A ( κ ) (cid:62) Y (cid:62) . (C.6) Finally, if each connected component of the reaction graph contains exactly one terminal component, then kerA ( κ ) (cid:62) Y (cid:62) = S ⊥ . (C.7) Proof.
First, we have that for any ˆ γ ∈ ˆΓ ij ( κ )ˆ γ + ker A ( κ ) (cid:62) Y (cid:62) ⊆ ˆΓ ij ( κ ) . (C.8)Indeed, for any ˆ γ ∈ ˆΓ ij ( κ ), there exists γ ∈ Γ ij ( κ ) with ˆ γ = π d ( γ ). Then, for any for any v ∈ ker A ( κ ) (cid:62) Y (cid:62) we have (cid:0) A ( κ ) (cid:62) Y (cid:62) | e i (cid:1) (cid:18) γ + (cid:18) v (cid:19)(cid:19) = (cid:0) A ( κ ) (cid:62) Y (cid:62) | e i (cid:1) γ + (cid:18) A ( κ ) (cid:62) Y (cid:62) v (cid:19) = e j . γ + (cid:18) v (cid:19) ∈ Γ ij ( κ )which implies that ˆ γ + v ∈ ˆΓ ij ( κ ) and proves (C.8). Since for any h ∈ S ⊥ and any 1 ≤ n ≤ m (cid:16) A ( κ ) (cid:62) Y (cid:62) h (cid:17) n = m (cid:88) l =1 (cid:104) h, y l − y n (cid:105) κ nl = 0 , it follows S ⊥ ⊆ ker A ( κ ) (cid:62) Y (cid:62) . (C.9)(C.4) follows from (C.8) and (C.9).Now assume that the mass-action system has a steady state c , with c y i >
0. Then,0 = ddt (cid:104) ˆ γ, x ( t ) (cid:105)| x ( t )= c = ˆ γ (cid:62) Y A ( κ )Λ( c ) = ( e j − γ d +1 e i )Λ( c ) = c y j − γ d +1 c y i . Since c y i >
0, then necessarily γ d +1 = c y j − y i . For the same argument, for any other γ (cid:48) ∈ Γ ij ( κ ), γ (cid:48) d +1 = γ d +1 = c y j − y i .Hence, (C.5) holds. It follows that any γ (cid:48) ∈ Γ ij ( κ ) is of the form γ (cid:48) = γ + (cid:18) v (cid:19) for some v ∈ R d . Moreover, since γ (cid:48) ∈ Γ ij ( κ ) e j = (cid:0) A ( κ ) (cid:62) Y (cid:62) | e i (cid:1) γ (cid:48) = (cid:0) A ( κ ) (cid:62) Y (cid:62) | e i (cid:1) γ + (cid:18) A ( κ ) (cid:62) Y (cid:62) v (cid:19) = e j + (cid:18) A ( κ ) (cid:62) Y (cid:62) v (cid:19) . Hence, necessarily v ∈ ker A ( κ ) (cid:62) Y (cid:62) , which impliesˆΓ ij ( κ ) ⊆ ˆ γ + ker A ( κ ) (cid:62) Y (cid:62) . The latter, together with (C.8), implies (C.6).To conclude the proof, we need to show that if each connected component of the reaction graph contains exactlyone terminal component (i.e. if (cid:96) = τ ), then (C.7) holds. This follows from (C.9) anddim S ⊥ = d − dim S = d + (cid:96) + δ − m = d + τ + δ − m ≥ d − m + dim ker Y A ( κ )= dim ker A ( κ ) (cid:62) Y (cid:62) , where we utilized Theorems B.1 and B.2 for the forth equality.As a consequence, we have the following. Corollary C.3.
Consider a mass-action system, and assume that y i and y j are two complexes for which ˆΓ ij ( κ ) isnon-empty. Let { v p } Hp =1 be a basis for S ⊥ , where H = d − dim S . Moreover, let X l , X l , . . . , X l n ∈ X such that therank H × n matrix V has rank n , where V pq = v pl q for all ≤ p ≤ H and ≤ q ≤ n . Then, there exists ˆ γ ∈ ˆΓ ij ( κ ) such that (cid:104) ˆ γ, e p (cid:105) = 0 for all ≤ p ≤ n .Proof. Let ˆ γ (cid:63) ∈ ˆΓ ij ( κ ). It follows from Proposition C.2 that for all w ∈ R H ˆ γ (cid:63) + H (cid:88) p =1 w p v p ∈ ˆΓ ij ( κ ) . Hence, the proof is concluded by choosing w such that w (cid:62) V = − (ˆ γ (cid:63)l , ˆ γ (cid:63)l , . . . , ˆ γ (cid:63)l n ) , which is possible because the H × n matrix V has rank n .Further useful structural property of ˆΓ ij ( κ ) are following. Such properties are useful while checking whether theresults of this paper can be applied, and can be used by an algorithm designed to this aim. Before stating thestructural results, it is convenient to prove the following lemma.21 emma C.4. Consider a mass-action system S = ( X , C , R , κ ) , and ≤ i, j ≤ m with i (cid:54) = j . Then, ˆΓ ij ( κ ) isnon-empty if and only if e j − γ d +1 e i ∈ (ker Y A ( κ )) (cid:62) for some γ d +1 ∈ R . Moreover, if that is the case and if themass-action system has a steady state c with c y i > and c y j > , then necessarily γ d +1 > .Proof. By definition of ˆΓ ij ( κ ) given in (C.2), ˆΓ ij ( κ ) is non-empty if and only if Γ ij ( κ ) is non-empty. By (C.1), Γ ij ( κ )is non-empty if and only if there exists γ d +1 ∈ R such that e j − γ d +1 e i is in Im A ( κ ) (cid:62) Y (cid:62) . By the fundamentaltheorem of linear algebra, the latter holds if and only if e j − γ d +1 e i is orthogonal to ker Y A ( κ ).To conclude the proof, we note that if the mass-action system has a steady state c with c y i > ij ( κ ) isnon-empty, then it follows by (C.5) in Proposition C.2 that the quantity γ d +1 discussed in the first part of the proofis necessarily equal to c y j − y i > Proposition C.5.
Consider a mass-action system S = ( X , C , R , κ ) , and ≤ i, j ≤ m . Assume that one of thefollowing holds: • y i is non-terminal and y j is terminal; • y i and y j are in two different terminal components.Then, ˆΓ ij ( κ ) is empty.Proof. By assumption, y j is terminal. By Theorem B.1, there exists a vector χ ∈ ker A ( κ ) ⊆ ker Y A ( κ )such that supp χ is the terminal component containing y j . Hence, if y i is non-terminal or if y i is in a different terminalcomponent than y j , we have that for all γ d +1 ∈ R (cid:104) χ, e j − γ d +1 e i (cid:105) = χ j (cid:54) = 0 . It follows from Lemma C.4 that ˆΓ ij ( κ ) is empty, which concludes the proof. Proposition C.6.
Consider a mass-action system S = ( X , C , R , κ ) , and ≤ i, j ≤ m with i (cid:54) = j . Assume that asteady state c with c y i > and c y j > exists. Then, ˆΓ ij ( κ ) is empty if and only if ˆΓ ji ( κ ) is empty.Proof. By the symmetric role of i and j , it suffices to prove that if ˆΓ ij ( κ ) is non-empty, then necessarily ˆΓ ij ( κ ) isnon-empty.By Lemma C.4, if ˆΓ ij ( κ ) is non-empty then there exists γ d +1 ∈ R with e j − γ d +1 e i ∈ (ker Y A ( κ )) (cid:62) , and γ d +1 > e i − γ d +1 e j ∈ (ker Y A ( κ )) (cid:62) , hence by Lemma C.4 ˆΓ ji ( κ ) is non-empty, which is what we wanted to show. D Proofs of the results stated in the main text
In this section, we state and prove more general versions of the theorems presented in the main text.
D.1 Existence and characterization of a CIF
The following result is stated in the main text.
Theorem 4.1.
Consider a mass-action system, and assume the following holds:1. there are two non-terminal complexes y i and y j such that only one entry of y j − y i is non-zero;2. the deficiency is 1.3. a positive steady state exists.Then, ˆΓ ij ( κ ) is non-empty. The following more general result holds, from which Theorem 4.1 can be immediately deduced.
Theorem D.1.
Consider a mass-action system, and assume the following holds:1. y i and y j are two distinct non-terminal complex; . the deficiency is 1;3. a steady state c with c y i (cid:54) = 0 exists.Then, ˆΓ ij ( κ ) is non-empty.Proof. By condition 3, c is a steady state, hence the vector Λ( c ) is in ker Y A ( κ ). Moreover, Λ i ( c ) = c y i (cid:54) = 0. Bycondition 1, y i is a non-terminal complex. Hence, by combining Λ i ( c ) (cid:54) = 0 with Theorem B.1, it follows that Λ( c ) isnot in ker A ( κ ). Since by assumption δ = 1, it follows from Theorem B.2 thatker Y A ( κ ) = span R { Λ( c ) } ⊕ ker A ( κ ) . (D.1)Define v = e j − c y j − y i e i = e j − Λ j ( c )Λ i ( c ) e i . Clearly, v is orthogonal to Λ( c ). Moreover, by condition 1, both y i and y j are non-terminal complexes. Hence, thevector v is orthogonal to ker A ( κ ) by Theorem B.1. It follows from (D.1) that v is orthogonal to ker Y A ( κ ), whichimplies that ˆΓ ij ( κ ) is non-empty by Lemma C.4.We now proceed to prove the following result, stated in the main text. Theorem 4.2.
Consider a mass-action system. Assume that there are two complexes y i and y j only differing in the n th entry, and that ˆΓ ij ( κ ) is non-empty. Let γ ∈ ˆΓ ij ( κ ) , and define q = γ yj − yi ) n d +1 . Then, either no positive steady state exists or the n th species is ACR with ACR value q . Moreover, φ ( x ) = d (cid:88) i =1 β i x i is a linear CI with ddt φ ( x ( t )) = Λ i ( x ( t )) (cid:16) x n ( t ) ( y j − y i ) n − q ( y j − y i ) n (cid:17) for any initial condition x (0) if and only if β ∈ ˆΓ ij ( κ ) . In order to prove the result, we will show that a more general version holds, which we state here.
Theorem D.2.
Consider a mass-action system, and assume that y i and y j are two complexes for which ˆΓ ij ( κ ) isnon-empty. Let γ ∈ ˆΓ ij ( κ ) . Then, c y j − y i = γ d +1 (D.2) for all steady states c satisfying c y i (cid:54) = 0 . Moreover, φ ( x ) = d (cid:88) i =1 β i x i is a linear CIF with ddt φ ( x ( t )) = Λ i ( x ( t )) (cid:16) x ( t ) y j − y i − γ d +1 (cid:17) (D.3) for any initial condition x (0) if and only if β ∈ ˆΓ ij ( κ ) .Proof. (D.2) follows from (C.5) in Proposition C.2.Let β ∈ ˆΓ ij ( κ ). Then, by Proposition C.2 we have that β − π d ( γ ) ∈ ker A ( κ ) (cid:62) Y (cid:62) . Hence, for any initial condition x (0) ∈ R d ≥ and any t ∈ R ≥ , ddt (cid:104) β, x ( t ) (cid:105) = β (cid:62) Y A ( κ )Λ( x ( t )) = π d ( γ ) (cid:62) Y A ( κ )Λ( x ( t )) = ( e j − γ d +1 e i ) (cid:62) Λ( x ( t )) = x ( t ) y i (cid:16) x ( t ) y j − y i − γ d +1 (cid:17) , which is (D.3).Conversely, assume that (D.3) holds. Then, for all x ∈ R d ≥ we have( β (cid:62) − π d ( γ )) (cid:62) Y A ( κ )Λ( x ) = 0 . Since the entries of Λ are linearly independent monomials on R d , it must be ( β (cid:62) − π d ( γ )) ∈ ker A ( κ ) (cid:62) Y (cid:62) , whichimplies that β ∈ ˆΓ ij ( κ ) by Proposition C.2. 23 .2 Rejection of persistent disturbances We recall here the formal definition of “oscillation”, as intended in this paper.
Definition D.1.
We say that a function g : R ≥ → R oscillates around a value q ∈ R if for each t ∈ R ≥ there exist t + > t and t − > t such that g ( t + ) > q and g ( t − ) < q. The following result holds.
Theorem D.3.
Consider a mass-action system, with associated differential equation ddt x ( t ) = f ( x ( t )) . Assume that y i and y j are two complexes for which ˆΓ ij ( κ ) is non-empty. Then, there exists q ∈ R ≥ such that c y j − y i = q for all steady states c with c y i (cid:54) = 0 . Consider an arbitrary function u : R d ≥ × R ≥ → R d such that asolution to ddt ˜ x ( t ) = f (˜ x ( t )) + u (˜ x ( t ) , t ) exists with ˜ x ( t ) ≥ for all t ≥ . Assume that there exists a ˆ γ ∈ ˆΓ ij ( κ ) which is orthogonal to the vector u ( x, t ) forany x ∈ R d ≥ , t ∈ R ≥ . Then, for any initial condition ˜ x (0) ∈ R d ≥ , at least one of the following statements holds:1. There is a species X k ∈ supp ˆ γ such that lim sup t →∞ ˜ x k ( t ) = ∞ ;2. There is a species X k ∈ supp y i such that lim inf t →∞ ˜ x k ( t ) = 0 ;3. ˜ x ( t ) y j − y i oscillates around q ;4. lim t →∞ (cid:90) ∞ t (cid:12)(cid:12)(cid:12) ˜ x ( s ) y j − y i − q (cid:12)(cid:12)(cid:12) ds = 0 .Proof. It follows from Theorem D.2 that there exists q ∈ R ≥ such that c y j − y i = q for all steady states c with c y i (cid:54) = 0.Note that in this setting c y i (cid:54) = 0 is equivalent to c y i >
0, since the state space is limited to vectors of non-negativeconcentrations.Fix an initial condition ˜ x (0) ∈ R d ≥ . Assume that 1 does not occur. Hence, there exists M ∈ R > such that |(cid:104) ˆ γ, ˜ x ( t ) (cid:105)| < M for all t ∈ R ≥ . By the fact that ˆ γ is orthogonal to u ( x, t ) for all x ∈ R d ≥ and all t ≥
0, and by Theorem D.2, we have ddt (cid:104) ˆ γ, ˜ x ( t ) (cid:105) = (cid:104) ˆ γ, f (˜ x ( t )) (cid:105) = ˜ x ( t ) y i (cid:16) ˜ x ( t ) y j − y i − q (cid:17) . It follows that (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) t ˜ x ( s ) y i (cid:16) ˜ x ( s ) y j − y i − q (cid:17) ds (cid:12)(cid:12)(cid:12)(cid:12) < M − |(cid:104) ˆ γ, ˜ x (0) (cid:105)| for all t ∈ R ≥ . (D.4)If 2 does not hold, then there exists m ∈ R > such that x ( t ) y i > m for all t ∈ R > . Together with (D.4), this wouldimply that there exists M (cid:63) ∈ R > such that (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) t (cid:16) ˜ x ( s ) y j − y i − q (cid:17) ds (cid:12)(cid:12)(cid:12)(cid:12) < M (cid:63) for all t ∈ R ≥ . (D.5)If 3 does not hold, then there exists t (cid:63) ∈ R > such that ˜ x ( t ) y j − y i − q maintains the same sign for all t > t (cid:63) , whichtogether with (D.5) implies 4. The proof is then concluded.Now we prove Theorem 6.1, which is stated in the main text and which we state here again for convenience.Theorem 6.1 follows almost entirely from Theorem D.3. Theorem 6.1.
Consider a mass-action system, with associated differential equation ddt x ( t ) = f ( x ( t )) . Assume that there are two complexes y i and y j only differing in the n th entry, and that ˆΓ ij ( κ ) is non-empty. Let q be the ACR value of the n th species. Consider an arbitrary function u with image in R d such that a solution to ddt ˜ x ( t ) = f (˜ x ( t )) + u (˜ x ( t ) , t ) exists. Assume that there exists a ˆ γ ∈ ˆΓ ij ( κ ) which is orthogonal to the vector u ( x, t ) for any x, t . Then, for anyinitial condition ˜ x (0) , at least one of the following holds: a) the concentration of some species goes to 0 or infinity, along a sequence of times;(b) ˜ x n ( t ) oscillates around q and ˆ γ k (cid:54) = 0 for some k (cid:54) = n ;(c) the integral (cid:90) ∞ t | ˜ x n ( s ) − q | ds tends to , as t goes to infinity.Proof. Assume (a) does not hold. Then, there exists ε ∈ R > such that ε ≤ (cid:107) ˜ x ( t ) (cid:107) ≤ ε for all t ∈ R ≥ . (D.6)Moreover, it follows from Theorem D.3 that at least one of the following holds:1. ˜ x n ( t ) ( y j − y i ) n oscillates around q ( y j − y i ) n , implying that ˜ x n ( t ) oscillates around q ;2. lim t →∞ (cid:90) ∞ t (cid:12)(cid:12)(cid:12) ˜ x n ( s ) ( y j − y i ) n − q ( y j − y i ) n (cid:12)(cid:12)(cid:12) ds = 0.Since 2 clearly implies (c), to complete the proof it suffices to show that 1 implies ˆ γ (cid:54) = ˆ γ n e n .Assume 1 holds. If it were ˆ γ = ˆ γ n e n , then by Theorem D.2 and by the fact that u ( x, t ) is orthogonal to ˆ γ for all x ∈ R d ≥ we would haveˆ γ n ddt ˜ x n ( t ) = ddt (cid:104) ˆ γ, ˜ x ( t ) (cid:105) = (cid:104) ˆ γ, f (˜ x ( t )) (cid:105) = ˜ x ( t ) y i (cid:16) ˜ x n ( t ) − q (cid:17) for all t ∈ R ≥ . By (D.6), there would be m, M ∈ R > such that m (cid:16) ˜ x n ( t ) − q (cid:17) ≤ ˆ γ n ddt (cid:16) ˜ x n ( t ) − q (cid:17) ≤ M (cid:16) ˜ x n ( t ) − q (cid:17) for all t ∈ R ≥ , which would in turn imply that ˜ x n ( t ) − q maintains the same sign for all t ∈ R ≥ . Hence, 1 could not hold and theproof is concluded. D.3 Inclusion in larger systems
Loosely speaking, a subsystem S (cid:48) of a reaction system S is simply the reaction system generated by a subset ofthe reactions of S , with the same choice of rate functions. In order to give a formal definition, some more care isneeded as the number of species involved in the subsystem may be smaller (hence the dimensions of S and S (cid:48) maybe different). Definition D.2.
Let S = ( X , C , R , { λ ij } ≤ i,j ≤ m ) be a reaction system with d species and m complexes. A subsystem of S is a reaction system S (cid:48) = ( X (cid:48) , C (cid:48) , R (cid:48) , { λ (cid:48) ij } ≤ i,j ≤ m (cid:48) ) with d (cid:48) species and m (cid:48) complexes such that, up to reorderingof X and C , • X (cid:48) = { X i ∈ X : 1 ≤ i ≤ d (cid:48) } ; • y in = 0 for all y i ∈ C with 1 ≤ i ≤ m (cid:48) and for all d (cid:48) < n ≤ d ; • C (cid:48) = { y (cid:48) i : y (cid:48) i = π ( y i ) , y i ∈ C , ≤ i ≤ m (cid:48) } , where π : R d → R d (cid:48) is the projection onto the first d (cid:48) components; • R (cid:48) ⊆ { y (cid:48) i → y (cid:48) j : y i → y j ∈ R} ; • for all 1 ≤ i, j ≤ m (cid:48) and all x ∈ R d ≥ λ (cid:48) ij ( π ( x )) = (cid:40) λ ij ( x ) if y (cid:48) i → y (cid:48) j ∈ R (cid:48) S and S (cid:48) , we further saythat y k → y l is a reaction of S but is not a reaction of S (cid:48) if y k → y j ∈ R and either max k, l > m (cid:48) or y (cid:48) k → y (cid:48) l / ∈ R (cid:48) .In this case, we write y k → y l ∈ R \ R (cid:48) .The following result is stated in the main text. 25 orollary 6.2. Consider a reaction system S , and let S (cid:48) be a sub-system. Assume that S (cid:48) is a mass-actionsystem with two complexes π ( y i ) and π ( y j ) only differing in the entry relative to the species X , and for which ˆΓ ij ( κ ) is non-empty. Moreover, assume there exists ˆ γ ∈ ˆΓ ij ( κ ) such that π ( y l − y k ) is orthogonal to ˆ γ for all y k → y l thatare reactions of S but not reactions of S (cid:48) . Then, the X is ACR for both S and S (cid:48) , with the same ACR value. Here, we prove the following, stronger version.
Corollary D.4.
Consider a reaction system S = ( X , C , R , { λ ij } ≤ i,j ≤ m ) , and let S (cid:48) = ( X (cid:48) , C (cid:48) , R (cid:48) , { λ (cid:48) ij } ≤ i,j ≤ m (cid:48) ) be a sub-system. Assume that S (cid:48) is a mass-action system with two complexes π ( y i ) and π ( y j ) for which ˆΓ ij ( κ ) isnon-empty. Moreover, assume there exists ˆ γ ∈ ˆΓ ij ( κ ) such that π ( y l − y k ) is orthogonal to ˆ γ for all y k → y l ∈ R \ R (cid:48) .Then, there exists a value q ∈ R ≥ such that ( c (cid:48) ) π ( y j − y i ) = q for all steady states c (cid:48) of S (cid:48) such that ( c (cid:48) ) π ( y i ) > , and c y j − y i = q for all steady states c of S such that c y i > .Proof. The existence of q ∈ R ≥ such that ( c (cid:48) ) π ( y j − y i ) = q for all steady states c (cid:48) of S (cid:48) with ( c (cid:48) ) π ( y i ) > ddt π ( x ( t )) = f ( π ( x ( t ))) + u ( x ( t ) , t ) , where f ( π ( x ( t ))) = (cid:88) y (cid:48) k → y (cid:48) l ∈R (cid:48) ( y (cid:48) l − y (cid:48) k ) κ ij π ( x ( t )) y i u ( x ( t ) , t ) = (cid:88) y k → y l ∈R\R (cid:48) π ( y l − y k ) λ ij ( x ( t ) , t ) . Note that ˆ γ is orthogonal to u ( x, t ), for all x ∈ R d ≥ and t ∈ R ≥ . Hence, if c is a steady state of S , it follows fromTheorem D.2 that 0 = ddt (cid:104) ˆ γ, π ( c ) (cid:105) = (cid:104) ˆ γ, f ( π ( c )) (cid:105) = π ( c ) π ( y i ) (cid:16) π ( c ) π ( y j − y i ) − q (cid:17) = c y i (cid:16) c y j − y i − q (cid:17) , which concludes the proof. E Example of an ACR system with no linear CIF
Assume that X i is an ACR species, with ACR value q . It is not always possible to find a linear combination ofspecies whose derivative at time t is of the form r ( t )( x i ( t ) α − q ), for some polynomial r ( t ) and some real number α .As an example, consider the following mass-action system:A+B κ B+C κ κ κ κ κ κ AD κ ELet us order the species alphabetically, and the complexes as ( A + B, B + C, B, B, E, D, C, A, D, E ). It is provenin [3] that the species A is ACR, with ACR value q = κ κ κ κ .
26n this case there is no linear combination of species whose derivative is of the form r ( t ) ( x ( t ) α − q α ) , (E.1)for some polynomial function r ( t ) and some real number α . Indeed, for any β ∈ R d , the derivative of (cid:104) β, x ( t ) (cid:105) isgiven by β (cid:62) Y A ( κ )Λ( x ( t )), which in this case is a polynomial of the form w x ( t ) x ( t ) + w x ( t ) + w x ( t ) x ( t ) + w x ( t ) + w x ( t ) + w x ( t ) + w x ( t ) , for some real coefficients w i . The only monomial containing x ( t ) is the first one, so if we want the derivative of (cid:104) β, x ( t ) (cid:105) to be of the form (E.1), necessarily w (cid:54) = 0, α = 1, and r ( t ) = w x ( t ). We can further deduce thatnecessarily w = w = w = w = w = 0, and w = − w q. In matrix form, this is equivalent to β (cid:62) Y A ( κ ) = w (1 , , , − q, , , , , , , but this is not possible because in this case the vector on the right-hand side does not belong to the left image of Y A ( κ ). F Applications to the systems introduced in the main text.
Here we use the theory developed in this work to analyze the reaction systems introduced in the main text.
F.1 Toy example
Consider the mass-action system A + B κ BB κ A Order the species alphabetically, and the complexes in the appearance order from left to right and from top to bottom,as A + B , 2 B , B , and A . In this case, we have Y = (cid:18) (cid:19) and A ( κ ) = − κ κ − κ
00 0 κ (F.1)The two non-terminal complexes differing for the entry relative to A are the first and the third ones. Since the modelhas deficiency 1, it follows from Theorem 4.1 that ˆΓ ( κ ) is non-empty. We can use Matlab to explicitly calculateit. First, in order to define Y and A ( κ ) in Matlab, we define the symbolic variables κ and κ and require they arepositive real numbers via k1 = sym ( ' k1 ' , ' positive ' ) ;k2 = sym ( ' k2 ' , ' positive ' ) ; We then define the matrices Y and Ak as in (F.1). We can calculate Ψ ( κ ) via e = eye (4) ;simplify ( null ([ Ak ' * Y ' e (: ,1) e (: ,3) ]) ) which returns the following matrix, whose columns are a basis for Ψ ( κ ): − κ − κ κ . It then follows from Proposition C.1 thatΓ ( κ ) = κ κ κ + a : a ∈ R , ( κ ) = κ + a : a ∈ R (F.2)and together with Theorem 4.2 that the ACR value of A is κ /κ . Note that in this case S = span R − and S ⊥ = span R . (F.3)Moreover, each connected component contains exactly one terminal component. Hence, Proposition C.2 applies, andas a matter of fact ˆΓ ( κ ) = κ + S ⊥ . Finally, while it is clear from (F.2) that a vector with zero first component and a vector with zero second componentare in ˆΓ ( κ ), this could be derived from Corollary C.3 and from (F.3) without explicitly calculating ˆΓ ( κ ). As aconsequence, due to Theorem 6.1, disturbances could be applied to the production and degradation rates of A (or B )while maintaining the absolute concentration robustness of the species A , its ACR value, and the linear CIF describedin Theorem 4.2. F.2 EnvZ-OmpR osmoregulatory system
Consider the mass-action systemEnvZ-D κ κ [ADP] EnvZ κ [ATP] κ EnvZ-T κ EnvZ-PEnvZ-P + OmpR κ κ EnvZ-OmpR-P κ EnvZ + OmpR-PEnvZ-D + OmpR-P κ κ EnvZ-OmpR-D-P κ EnvZ-D + OmpROrder both the species and the complexes according to their appearance order, from left to right and from top tobottom. Hence, the 8 species are ordered as EnvZ-D, EnvZ, EnvZ-T, EnvZ-P, OmpR, EnvZ-OmpR-P, OmpR-P, andEnvZ-OmpR-D-P. The 10 complexes are ordered as EnvZ-D, EnvZ, EnvZ-T, EnvZ-P, EnvZ-P+OmpR, EnvZ-OmpR-P, EnvZ+OmpR-P, EnvZ-D+OmpR-P, EnvZ-OmpR-D-P, and EnvZ-D+OmpR.It can be checked that dim S = 6 and S ⊥ = span R , . (F.4)The above conservation laws correspond to the fact that the total mass of the chemical species containing some formof EnvZ is conserved, as well as the total mass of the chemical species containing some form of OmpR. We cancalculate the deficiency as δ = m − (cid:96) − dim S = 10 − − . The non-terminal complexes that only differ for the entry relative to OmpR-P are the first and the eighth ones.Therefore, we are interested in the study of ˆΓ ( κ ), which by Theorem 4.1 is not empty. We have Y = A ( κ ) = − κ κ [ADP] 0 0 0 0 0 0 0 0 κ − κ [ADP] − κ [ATP] κ κ [ATP] − κ − κ κ − κ κ κ − κ − κ κ − κ κ
00 0 0 0 0 0 0 κ − κ − κ
00 0 0 0 0 0 0 0 κ . In order to define the corresponding symbolic matrices Y and Ak in Matlab, we can first define the symbolic positivereal variables k = sym ( ' k ' , [1 ,11] , ' positive ' ) ;k (2) = k (2) * sym ( ' ADP ' , ' positive ' ) ;k (3) = k (3) * sym ( ' ATP ' , ' positive ' ) ; Then we calculate Ψ ( κ ), as defined in (C.3), via e = eye (10) ;simplify ( null ([ Ak ' * Y ' e (: ,8) e (: ,1) ]) ) The output of the last command is the following matrix, whose columns are a basis of Ψ ( κ ). − − [ADP] κ κ ( κ + κ )[ATP] κ κ κ ( κ + κ ) − − [ADP] κ κ κ +[ADP] κ κ κ +[ATP] κ κ κ +[ATP] κ κ κ [ATP] κ κ κ ( κ + κ ) − − [ADP] κ κ κ +[ADP] κ κ κ +2 [ADP] κ κ κ +[ATP] κ κ κ +[ATP] κ κ κ [ATP] κ κ κ ( κ + κ ) − − [ADP] κ κ κ +2 [ADP] κ κ κ +[ADP] κ κ κ +2 [ADP] κ κ κ +[ATP] κ κ κ +[ATP] κ κ κ [ATP] κ κ κ ( κ + κ ) [ADP] κ ( κ + κ )[ATP] κ κ κ − [ADP] κ κ κ +[ADP] κ κ κ +[ATP] κ κ κ +[ATP] κ κ κ [ATP] κ κ κ ( κ + κ ) − [ADP] κ κ κ ( κ + κ )[ATP] κ κ κ ( κ + κ ) For convenience, denote the last vector by ζ ( κ ). It follows from Proposition C.1 thatΓ ( κ ) = − π ( ζ ( κ )) + a − − − − + a : a , a ∈ R .
29t follows that ˆΓ ( κ ) = − π ( ζ ( κ )) + a − − − − + a : a , a ∈ R , (F.5)which corresponds to (11) in the main text. Moreover, it follows from Theorem 4.2 that the ACR value of the speciesOmpR-P is − ζ ( κ ) = [ADP] κ κ κ ( κ + κ )[ATP] κ κ κ ( κ + κ )as reported in (10) in the main text.Further things can be noted about this model. First, it follows from (F.4) and (F.5) thatˆΓ ( κ ) = − π ( ζ ( κ )) + S ⊥ , which is in accordance with Proposition C.2 because all connected components contain exactly one terminal com-ponent. Secondly, due to Corollary C.3, we can deduce without explicitly calculating Γ ( κ ) that there is a vectorˆ γ in ˆΓ ( κ ) whose entries relative to species EnvZ and OmpR-P are zero (which are the second and the seventhcomplexes, respectively). Specifically, if we consider the basis of S ⊥ given in (F.4) and we let X l = X =EnvZ and X l = X =OmpR-P, we have V = (cid:18) (cid:19) , which has rank 2. The existence of such vector ˆ γ is used in Section 6.3 and it implies by Theorem 6.1 that the systemis robust to persistent disturbances affecting the production and degradation rates of both EnvZ and OmpR-P. G An ACR signaling system covered by our theory and not by [40]
Consider the double-phosphorylation mass-action system in Figure 6. We will show that the theory developed in [40]stays silent on whether it is ACR. However, our theory covers this case and implies the double-phosphorylated formof the transcriptional regulatory protein is ACR, for any choice of kinetic parameters κ such that a positive steadystate exists. Note that this form of response robustness may seem a bit surprising for a multisite phosphorylationmechanism, since these are often known for their multi-stability properties, notably shown in the case of the MAPKpathway [7, 26, 27, 44].Order the species and the complexes according to their appearance order, from left to right and from top tobottom. The 11 species are then ordered as A , A (cid:63) , A -P, B , A - B -P, B -P, A - B -PP, B -PP, A (cid:63) - B -P, A (cid:63) - B -PP, and B -PP- A . The 13 complexes are ordered as A , A (cid:63) , A -P, A -P + B , A - B -P, A + B -P, A -P + B -P, A - B -PP, A + B -PP, A (cid:63) - B -P, A (cid:63) - B -PP, A -P + B -PP, and B -PP- A . It can be checked that dim S = 9 and S ⊥ = span R , . (G.1)Similarly as for the EnvZ-OmpR osmoregulatory system, the above conservation laws correspond express that thetotal mass of the chemical species containing some form of the protein A is conserved, as well as the total mass ofthe chemical species containing some form of the protein B . The reaction graph is given in Figure 7. In particular, (cid:96) = 2 and the deficiency is δ = m − (cid:96) − dim S = 13 − − . So, this model is not included in the class of models studied in [40], and Theorems 3.1 and B.3 do not apply.30ctivation and phosphorylation of the sensor-transmitter protein:
A κ [I] κ A (cid:63) κ A - P Phosphorylation of the sensory response protein: A - P + B κ κ A - B - P κ A + B - P A - P + B - P κ κ A - B - PP κ A + B - PP Activation and phosphorylation of the protein aggregate: A - B - P κ [I] κ A (cid:63) - B - P A - B - PP κ [I] κ A (cid:63) - B - PP A (cid:63) - B - P κ A (cid:63) - B - PP κ A - P + B - PP Dephosphorylation of the sensory response protein: A + B - PP κ κ B - PP - A κ A - P + B Figure 6: A model of signal transduction. Here, the sensory response protein B is active if two phosphorylgroups are attached. In the first line, the protein A is phosphorylated in response to a stimulus [I]. In the secondand third line, the phosphoryl groups are transferred to the sensory response protein B . In the forth and fifthline, the protein A responds by the stimulus [I] while it is bound to B . The resulting protein aggregate is ableto gain phosphoryl groups until the three phosphoryl sites are occupied (as described in the sixth line). Finally,in the last two lines, the inactive form of the sensor-transmitter protein A acts as a phosphatase on B .The first and the ninth complexes are non-terminal, and they only differ for the entry relative B-PP. Hence, inorder to apply Theorems 4.2 and 6.1, we are interested in studying ˆΓ ( κ ). Since the deficiency of the model is 2,Theorem 4.1 does not apply, so to understand whether ˆΓ ( κ ) is non-empty we need to explicitly calculate it. Tothis aim, we define in Matlab the following positive symbolic variables, which correspond to the rate constants of themodel: k = sym ( ' k ' , [1 ,18] , ' positive ' ) ;k (1) = k (1) * sym ( ' I ' , ' positive ' ) ;k (10) = k (10) * sym ( ' I ' , ' positive ' ) ;k (12) = k (12) * sym ( ' I ' , ' positive ' ) ; κ [I] κ A (cid:63) κ A - P A - P + B κ κ A - B - P κ A + B - P κ [I] κ A (cid:63) - B - P κ A (cid:63) - B - PP κ κ [I] A - B - PP κ A - P + B - PP A + B - PP B - PP - A κ κ A - P + B - P κ κ κ κ Figure 7: Reaction graph for the signaling transduction system proposed in Figure 6.We then define the symbolic matrices Y and Ak corresponding to Y = and to A ( κ ) as described in Figure 8, respectively.In order to calculate a basis for Ψ ( κ ) as defined in (C.3), we type e = eye (13) ;simplify ( null ([ Ak ' * Y ' e (: ,1) e (: ,9) ]) ) The output of the last command is shown in Figure 9. While the output is rather complicated, we did not need toput much effort in its calculation, which was completed in Matlab in a matter of seconds. For convenience we denoteby ζ ( κ ) the last column of the output matrix. 32t follows from Proposition C.1 that Γ ( κ ) is non-empty, and is given byΓ ( κ ) = − π ( ζ ( κ )) + a − − − + a : a , a ∈ R . (G.2)It follows from (G.2) and (G.1) thatˆΓ ( κ ) = − π ( ζ ( κ )) + a − − − + a : a , a ∈ R = − π ( ζ ( κ )) + S ⊥ , which is in accordance with Proposition C.2 because every connected component of the reaction graph in Figure 7contains exactly one terminal component. Moreover, from (G.2) and Theorem 4.2 it follows that the ACR value ofthe ACR species B -PP is − ζ ( κ ) = κ κ [I]( κ + κ )( κ κ + κ κ + κ κ [I])( κ κ + κ κ + κ κ [I]) κ κ ( κ + κ ) g ( κ, [I]) , where g ( κ, [I]) = − κ κ κ κ [I] − κ κ κ κ [I] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ . This means that if a positive steady state c exists, necessarily it entry relative to B -PP (which is the eighth species)satisfies c = − ζ ( κ ). Of course, this cannot occur if − ζ ( κ ) is non-positive, i.e. if g ( κ, [I]) is non-positive, in whichcase no positive steady state exists. Note that we did not need to work directly with the differential equation to derivethis information. We will further show that in fact a positive steady state exists if and only if g ( κ, [I]) is positive. Tothis aim, note that c is a steady state if and only if Λ( c ) ∈ ker Y A ( κ ). Then, we calculate a basis for ker Y A ( κ ) bytyping simplify ( null ( Y * Ak ) ) which returns a matrix of the form b ( κ, [I]) g ( κ, [I])0 0 0 b ( κ, [I]) g ( κ, [I])1 0 0 00 0 0 b ( κ, [I])0 0 0 b ( κ, [I])0 1 0 00 0 0 b ( κ, [I])0 0 0 b ( κ, [I])0 0 0 b ( κ, [I])0 0 0 b ( κ, [I])0 0 0 b ( κ, [I])0 0 1 00 0 0 1 b i ( κ, [I]) which map positive arguments to positive real numbers. Hence, a positive steady stateexists if and only if there exists a positive vector c such that c = a b ( κ, [I]) g ( κ, [I]) c c = a c = a b ( κ, [I]) c = a b ( κ, [I]) g ( κ, [I]) c c = a b ( κ, [I]) c = a b ( κ, [I]) c = a c = a b ( κ, [I]) c c = a c c = a b ( κ, [I]) c c = a b ( κ, [I]) c = a c = a b ( κ, [I])for some a , a , a , a ∈ R > . If g ( κ, [I]) is positive, it is easy to see that such a positive vector c exists. Specifically,the positive steady states are parameterized by c = a b ( κ, [I]) g ( κ, [I]) c = a b ( κ, [I]) c = a b ( κ, [I]) c = a b ( κ, [I]) g ( κ, [I]) c = a a b ( κ, [I]) c = a b ( κ, [I]) c = a c = a b ( κ, [I]) c = a c = a a b ( κ, [I]) c = b ( κ, [I]) b ( κ, [I]) g ( κ, [I])where a and a vary in R > . Note that the entry c is relative to the ACR species B -PP and is fixed.34 ( κ ) = − κ [ I ] κ κ [ I ] − κ − κ κ − κ κ κ κ − κ − κ − κ [ I ] κ
000 0000 κ − κ κ κ − κ − κ − κ [ I ] κ
00 0000000 κ − κ κ κ [ I ] − κ − κ
000 0000000 κ [ I ] κ − κ − κ
00 0000000000 κ
00 00000000 κ − κ − κ F i g u r e : A ( k ) f o r t h e m a ss - a c t i o n s y s t e m i n F i g u r e . − κ − p ( κ , [ I ] ) κ κ ( κ + κ ) ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ ) − p ( κ , [ I ] ) κ κ ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ ) − ( κ + κ )( κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ + [ I ] κ κ κ κ + [ I ] κ κ κ κ + [ I ] κ κ κ κ ) κ κ ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ ) − κ κ κ − κ ( κ + κ )( κ + κ )( κ κ + κ κ + [ I ] κ κ ) κ κ ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ ) p ( κ , [ I ] ) κ κ ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ )
100 01 κ κ κ κ κ − κ κ κ κ κ − κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ κ κ ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ ) p ( κ , [ I ] ) κ κ ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ )
010 00 − [ I ] κ κ ( κ + κ )( κ κ + κ κ + [ I ] κ κ )( κ κ + κ κ + [ I ] κ κ ) κ κ ( κ + κ ) ( − κ κ κ κ [ I ] − κ κ κ κ [ I ] + κ κ κ κ + κ κ κ κ + κ κ κ κ + κ κ κ κ ) w h e r e p ( κ , [ I ] ) = κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ − [ I ] κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ − [ I ] κ κ κ κ κ κ + [ I ] κ κ κ κ κ κ p ( κ , [ I ] ) = κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ p ( κ , [ I ] ) = κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ − [ I ] κ κ κ κ κ p ( κ , [ I ] ) = κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ + [ I ] κ κ κ κ κ F i g u r e : k e r Ψ ( κ ) f o r t h e m a ss - a c t i o n s y s t e m i n F i g u r e ..