A framework for deriving analytic steady states of biochemical reaction networks

The long-term behaviors of biochemical systems are often described by their steady states. Deriving these states directly for complex networks arising from real-world applications, however, is often challenging. Recent work has consequently focused on network-based approaches. Specifically, biochemical reaction networks are transformed into weakly reversible and deficiency zero generalized networks, which allows the derivation of their analytic steady states. Identifying this transformation, however, can be challenging for large and complex networks. In this paper, we address this difficulty by breaking the complex network into smaller independent subnetworks and then transforming the subnetworks to derive the analytic steady states of each subnetwork. We show that stitching these solutions together leads to the analytic steady states of the original network. To facilitate this process, we develop a user-friendly and publicly available package, COMPILES (COMPutIng anaLytic stEady States). With COMPILES, we can easily test the presence of bistability of a CRISPRi toggle switch model, which was previously investigated via tremendous number of numerical simulations and within a limited range of parameters. Furthermore, COMPILES can be used to identify absolute concentration robustness (ACR), the property of a system that maintains the concentration of particular species at a steady state regardless of any initial concentrations. Specifically, our approach completely identifies all the species with and without ACR in a complex insulin model. Our method provides an effective approach to analyzing and understanding complex biochemical systems.

Introduction (Line 57): It is particularly important to determine properties of the steady states of ODEs because they often describe long-term behaviors of the CRNs.
Methods (Line 378): The long-term behavior of a system is often described by steady states, which could be solved by equating each time derivative to zero, i.e., da dt = db dt = 0, and solving for the concentrations of the species in terms of the rate constants.

Reviewer 2 (Dr. Pavel Kraikivski) Summary
In the manuscript "A framework for deriving analytic long-term behavior of biochemical reaction networks", Hernandez et al. have a computational framework to analyze a molecular interaction network by breaking it into smaller independent subnetworks to derive steady states of the system. Authors have applied their computational tool to analyze toggle switch and insulin signaling models. It is very useful method and the computational tool which can be broadly applied. Thus, this work is well suitable for PLoS Computational Biology readers.
We appreciate the kind comments of the reviewer Dr. Pavel Kraikivski.

Minor comment 1
If there are limitations when the molecular interaction network cannot be analyzed with this approach or limits for the computational package applicability, I suggest to add such information to the discussion section.
We thank the reviewer for pointing this out. We indicated the limitations of the study in the Discussion section (Lines 342-347) as follows: Our work focuses on the derivation of steady states and its usefulness in analyzing biochemical systems. It would be an interesting future work to extend our framework to analyze the other long-term behaviors of networks, such as stability of the steady states [1,4,8], boundary steady states [2,7], and oscillations [3]. In particular, for linear stability analysis of the complex-balanced equilibria of generalized networks, Boros, Müller and Regensburger [5,17] have recently proposed an interesting approach.

Minor comment 2
In line 101: "B + C → C", should it be B + C → B as Figure 1 shows this reaction?
We thank the reviewer for bringing this to our attention. The comment is correct and we have revised this, i.e., "B +C → C" was changed to "B +C → B". This is a typographical error for particular line only and will not affect the rest of the manuscript.
The revision is given as follows: Results (Line 110): the reaction B + C → B is not contained in a cycle.

Reviewer 3 Summary
The authors develop a software package to symbolically parametrize the set of steady states of chemical reaction networks (CRNs) with mass-action kinetics. Theoretically, their approach is based on three ingredients: 1. network decomposition (into independent subnetworks) 2. network translation (into dynamically equivalent generalized chemical reaction networks (GCRNs)) 3. previous results on GCRNs.
Having such a software package is definitely useful for the CRN community and modelers in systems biology.
We thank the reviewer for acknowledging the importance of computational packages in the CRN community and to modelers in systems biology.
However, items 2 and 3 are not (yet) described clearly enough and impossible to understand (for readers without background in GCRNs -and even for me, as a developer of GCRN theory). This should be improved "everywhere": in abstract, intro, results, methods, ...
We appreciate the comments to improve the clarity of the manuscript. We have addressed these in detail (see below).
Another obvious observation concerns the title -it promises too much: Computing steady states is not the same as determining the "analytic long-term behavior". What about stability of steady states? (And periodic orbits? Homo/heteroclinic orbits? Chaos?) The title should be modified, and the stability issue should be discussed; see e.g. the recent paper Müller/Regensburger (2022) arxiv:2212.11039. In my opinion, the paper can be published after a major revision.
We provided the following paragraph in the Discussion section (Lines 337-342) to state the limitations of the study: Our work focuses on the derivation of steady states and its usefulness in analyzing biochemical systems. It would be an interesting future work to extend our framework to analyze the other long term behaviors of networks, such as stability of the steady states [1,4,8], boundary steady states [2,7], and oscillations [3]. In particular, for linear stability analysis of the complex-balanced equilibria of generalized networks, Boros, Müller and Regensburger [5,17] have recently proposed an interesting approach.
Additionally, we also changed the title to: "A framework for deriving analytic steady states of biochemical reaction networks".
In general, however, CRNs arising in applications are rarely WR or DZ. Hence, the method of network translation, that can modify the network structure while maintaining the dynamics, has been developed [11][12][13][14]19]. It involves shifting reactions, i.e., swapping reactions for ones with the same stoichiometric differences while keeping the original reaction rates. For instance, a reaction A+B → B could be shifted to A → 0 while we keep the rate of the original reaction k 1 ab to maintain the original dynamics. This network translation leads to generalized chemical reaction network (GCRN) with two associated structures. The first structure is the network's stoichiometric CRN with nodes that include the new stoichiometric nodes A and 0 from the shifted reaction, and the kinetic-order CRN with nodes that include the original source node A+B from the original reaction. The deficiencies of the former and the latter network structures of the translated network are known as the effective and kinetic deficiencies, respectively. If the stoichiometric network is WR and effective deficiency zero, and the kinetic order network is WR and kinetic deficiency zero, the translated GCRN is WR and DZ. In this case, a parametrization of the steady states can be computed easily [14].

Major comment 2
Results: paragraphs 110-120 and 120-135: the description of item 2 (network translation into dynamically equivalent GCRNs) in Figs 1 and 2 is impossible to follow.
Please improve the presentation! (See my suggestions on the next page!) We did not follow exactly the presentation provided in literature [14] to simplify the description for a wide audience/readers of computational scientists (biologists) of PLOS Computational Biology. Nevertheless, we believe that accurate description is more important as the reviewer pointed out. Thus, we have now followed the terminologies of the GCRN theory as suggested by the reviewer. In particular, we edited both Figs 1 and 2, their captions, and the corresponding text to include the stoichiometric and kinetic-order CRNs arising from the network translation as a generalized network as follows: Results, Lines 119-131: The first subnetwork (N 1 ) (Fig. 1b left) is a WR and DZ network.
On the other hand, the second subnetwork (N 2 ) (Fig. 1b middle) is neither WR nor DZ. Thus, we have to translate N 2 to make it WR and DZ. Specifically, N 2 is not WR because the reaction B + C → B does not belong to a cycle. Thus, we replace it with A + C → A, having the same stoichiometry to obtain the network, via the process of network translation, referred to as the stoichiometric CRN and denoted by N ′S 2 (Fig. 1b middle) (see the Methods section and [11] for details on network translation). We keep the original reaction rate (k 2 bc) for A + C → A so that the dynamics does not change. However, as the network structure changes, and in particular, as all the reactions now belong to a cycle, the stoichiometric CRN is a WR network. Furthermore, the deficiency of the stoichiometric CRN, called the effective deficiency, is zero, i.e., EDZ, because it has three nodes (n = 3), one connected component (ℓ = 1), and a stoichiometric matrix of rank two (s = 2) so that δ = 3 − 1 − 2 = 0.
Results, Lines 120-135 were revised as follows: Results, Lines 132-143: Next, we construct the kinetic-order CRN, which is denoted by N ′K 2 . To do this, we first take all the edges of N ′S 2 . Then, we fill 0 in the tail of the edge associated with k 2 (Fig. 2a (iii)) because the source node of the edge associated with k 2 in the original network (N 2 ) is 0 ( Fig. 2a (i)). Similarly, we fill A in the tail of the edge associated with k 3 (Fig. 2a (iii)). Because the source nodes of the edges with k 1 and k 4 are A + C and B + C in the original network ( Fig. 2a (i)), two source nodes (A + C and B + C) are now placed in a single node (Fig. 2a (iii)). To separate the two source nodes, a phantom edge, which is an edge with a free parameter (σ 1 ), is introduced in such a way that the resulting network is still WR (Fig. 2a (iv)) (see the Methods section and [13] for independent subnetworks (N 1 and N 2 ) in such a way that the rank of the original matrix (s) is equal to the sum of the ranks of the stoichiometric matrices of subnetworks (s 1 + s 2 ). b Since the subnetwork N 2 is not WR and DZ, network translation is performed. The translated subnework (N ′ 2 ) is WR and DZ, i.e., its stoichiometric and kinetic-order CRNs (N ′S 2 and N ′K 2 , respectively) are both WR and DZ, while its dynamics is equivalent to the dynamics of the original subnetwork (N 2 ). Then, the steady states of the original subnetwork (N 1 ) and the translated subnetwork (N ′ 2 ) can be analytically derived since they are WR and DZ (see Fig. 2 for details). c The steady states of subnetworks are merged to identify the steady states of the original networks. In particular, the steady states of the common species of N 1 and N 2 (i.e., species B) are equated, which eliminates the free parameter τ 1 . Then, by combining the steady state of every species, the steady state of the original network can be derived with one remaining free parameter σ 1 . This σ 1 is computed from the conserved quantity in the network, which is the sum of the initial concentrations of species B and D (i.e., b 0 + d 0 ). (ii) have corresponding source nodes A + C and B + C, respectively, in the original network (i). Since there are two source nodes (A + C and B + C) which share a single node in (iii), a phantom edge with a zero stoichiometric vector and a free parameter rate constant (σ 1 ) is introduced (upper right). Notice that this operation maintains the dynamics of the original network. Proceed to the next step since the deficiency of N ′K 2 is zero (i.e., KDZ). b Collect all the spanning trees of N ′K 2 (i.e., connected subgraphs of the N ′K 2 without a cycle) which direct towards each node (B + C, A + C, 0, and A). Then multiply the rate constants associated with the edges of each spanning tree. For instance, the products of the rate constants associated with the edges of the two spanning trees towards node A are σ 1 k 2 k 4 and k 1 k 2 k 4 (bottom). Compute each tree constant (K i ) by adding the product of the rate constants over all the spanning trees towards node (i). Hence, the tree constant K A = σ 1 k 2 k 4 + k 1 k 2 k 4 is obtained. c Choose an arbitrary tree (i.e., a connected graph without a cycle) containing all the nodes in N ′K 2 , which is not necessarily its subgraph. Then, for each edge (i → i ′ ) of the chosen tree, find the ratios of the tree constants . For instance, the edge A + C → 0 has the ratio of the tree Specifically, raise the ratios of the tree constants κ i→i ′ in a component-wise manner to the entries of a row of H, and get their product. For instance, the ratio of the tree constants κ A+C→0 , κ A+C→A , and κ A+C→B+C are raised to the entries in the first row of H (−1, 1 and 0, respectively). This gives κ −1 A+C→0 κ 1 A+C→A κ 0 A+C→B+C . Meanwhile, one free parameter (τ 1 ) is introduced because the number of column of B is one. This τ 1 is raised to the entry of a row in B. Hence, τ 0 1 is obtained. Finally, get the product κ −1 A+C→0 κ 1 A+C→A κ 0 A+C→B+C τ 0 1 , which is the steady state value of species A. The steady state values for species B and C can be computed in a similar manner. details on phantom edges). For the phantom edge, a zero stoichiometric vector is always assigned to maintain the dynamics of the original network. This produces a kinetic-order CRN N ′K 2 .

Major comment 3
Results: remaining paragraphs 135-210: the description of item 3 (main results for GCRNs) is missing in the Methods section. So the "description in words" (of effective deficiency, kinetic deficiency, "free parameters", tree constants, particular solution, generalized inverse, parametrization, ...) in the Results section is impossible to understand.
We added a section about the generalized chemical reaction networks in the Methods (Lines 408-445) section to discuss the Results section in a more understandable manner as suggested by the reviewer. This includes important descriptions, definitions and results that are essential in explaining the required concepts of effective and kinetic deficiencies, free parameters and previous results in parametrization of positive steady states of Johnston et al. [14]. In addition, we also described tree constant and generalized inverse in Theorem 2. Moreover, we indicated further references of tree constants in the Results section.

Generalized chemical reaction networks
In this section, we formally define the important concept of the generalized chemical reaction networks (GCRNs) pioneered by Müller and Regensburger [15,16].
Let G = (V, E) be a directed graph with vertex set V and edge set E ⊆ V × V . Additionally, let V s = {i|i → j ∈ E}, which is the set of all the source vertices.
Definition 1 A generalized chemical reaction network (GCRN) is a directed graph G = (V, E) together with two maps 1. y : V → R m ≥0 that assigns each vertex to a stoichiometric complex, and 2. y : V s → R m ≥0 that assigns each source vertex to a kinetic complex.
Given a CRN with an associated graph G, a dynamically equivalent GCRN (i.e., the associated ODEs of the GCRN agree with the ODEs of the original CRN) is a graph G ′ together with two maps y andỹ, and hence, with two sets of complexes. It gives rise to two associated CRNs, the stoichiometric one (G ′ , y) and a kinetic-order one (G ′ , y).
We call the deficiencies associated to the stoichiometric CRN and kinetic-order CRN, the effective deficiency and kinetic deficiency, respectively. If there is a map between the reactions of the original CRN and the stoichiometric CRN that preserves the reaction vectors and relates the source complexes in the original CRN to the kinetic complexes in the kinetic-order CRN, a GCRN is a network translation of a given CRN (see [11,14]).
We now introduce the following definition of a phantom edge and then proceed with the important result of parametrization of positive steady states introduced by Johnston et al.
Definition 2 For a given GCRN, we call an edge that connects identical stoichiometric complexes a phantom edge. Otherwise, we call it an effective edge.
A phantom edge does not contribute to the associated ODEs because the associated stoichiometric complexes are identical. Hence, the associated dummy reaction rate constant, denoted by σ, can be arbitrary and can be considered a free parameter. We denote the sets of phantom edges and effective edges of the GCRN by E 0 and E * , respectively.
where κ(k * , σ) H ⊤ • τ B ⊤ is the Hadamard product (the number of components of κ is the total number of effective and phantom edges) with the component of κ associated with the edge i → i ′ as κ i→i ′ = K i ′ K i and tree constant K i as the sum (over all the spanning trees of the kinetic-order CRN towards node i) of the products of the rate constants associated with the edges of each spanning tree. In addition, if the effective deficiency is zero, then the set of positive steady states of the original network is preciselyZ.
The theorem stated above covers the case when the kinetic deficiency is zero. For generalized networks with positive kinetic deficiency, additional conditions need to be checked before a steady state parametrization may be constructed (see Theorem 15 of [14]).
We also made revisions in paragraphs in lines 135-210 according to essential terms in the GCRN theory: Results, Lines 144-222: Then, we compute the deficiency of N ′K 2 , i.e., the kinetic deficiency. The differences of the nodes of the edges except for the phantom edge are the same as the reaction vectors of the stoichiometric CRN, and the difference of the nodes of the phantom edge ((B + C) − (A + C) = [−1, 1, 0, 0] ⊤ ) is independent from the reaction vectors of the stoichiometric CRN. Thus, the rank of the matrix of reaction vectors of N ′K 2 is one more than the rank of the stoichiometric matrix of N ′S 2 (s = 3). Furthermore, N ′K 2 has one more node (B + C) compared to N ′S 2 (n = 4). Additionally, the number of connected components are the same (ℓ = 1). Since δ = n − ℓ − s = 4 − 1 − 3 = 0, the kinetic deficiency of N ′K 2 is zero (i.e., KDZ). We proceed to the next step as the kinetic deficiency is zero; otherwise, additional conditions need to be satisfied (see [13,14] To calculate the parametrization of the steady states, we need to get all the spanning trees of N ′K 2 , i.e., connected subgraphs without a cycle directed towards each node as shown in Fig. 2b. For instance, there is only one spanning tree directed towards B + C because there is only one connected subgraph of the kinetic-order CRN without a cycle that points towards B + C (Fig. 2a (iv)). Then, to get the tree constant (K i ) associated with a particular node i with only one spanning tree, we multiply the rate constants corresponding to the edges of this spanning tree (see the description of tree constants in Theorem 2 of the Methods section and [11] for more details). Hence, we obtain K B+C = σ 1 k 2 k 4 as the tree constant directed towards the node B + C. In general, if there is more than one spanning tree directed towards a particular node, we have to repeat the process of getting the corresponding product for each spanning tree and then get the sum of the products over all the spanning trees directed towards the particular node. Specifically, there are two spanning trees directed towards node A because there are two connected subgraphs of N ′K 2 without a cycle that point towards node A. Then, multiply the rate constants associated with the edges of each spanning tree. Hence, we obtain σ 1 k 2 k 4 and k 1 k 2 k 4 (bottom) as the products of the rate constants associated with the two spanning trees directed towards node A. Thus, the tree constant associated with node A is K A = σ 1 k 2 k 4 + k 1 k 2 k 4 . Note that the term "tree constant" was introduced in [11] to refer to the formulas following from the matrix tree theorem presented in [18] and adapted to chemical reaction network theory in [6].
Next, from N ′K 2 ( Fig. 2a (iv)), we form a tree (i.e., a connected graph without a cycle) that contains all the nodes of this CRN. The possible trees that can be formed are given in Fig. 2c (top). Then, for each edge of the tree, find the ratio of the tree constants (center) and the kinetic difference (middle right), the ratio of the tree constants κ A+C→0 = From the matrices H and B together with the ratios of the tree constants, we can derive the analytic steady state of the subnetwork N 2 (Fig. 2d). In particular, to get the analytic steady state of the first species A, we raise each ratio of the tree constants associated with the three edges (i.e., κ A+C→0 , κ A+C→A , and κ A+C→B+C ) to each of the entries in the first row of H (i.e., −1, 1, and 0, respectively), and then multiply them to obtain κ −1 A+C→0 · κ 1 A+C→A · κ 0 A+C→B+C . Additionally, we introduce a number of free parameters, as many as the number of columns of B. We have a single free parameter (τ 1 ) because there is only one column of B. We raise this free parameter to the value of the first row of B (i.e., 0) and we get τ 0 1 . We multiply this value from the previously obtained product, which gives (κ −1 A+C→0 · κ 1 A+C→A · κ 0 A+C→B+C ) · (τ 0 1 ). This is the value of the analytic steady state of the first species A. By following the same procedure, we can get the analytic steady state values of the two remaining species of the translated subnetwork (N ′ 2 ). Because N 1 is a WR and DZ network, we did not need the network translation. In this case, we can get its steady state by applying the procedure (Fig. 2) to N 1 . While we illustrate the parametrization of steady states with the simple example, a rigorous description is provided in the Methods section (see Theorem 2).
Then, we combine the steady states of the two subnetworks (N 1 and N 2 ) to derive the steady state of the original whole network N . Specifically, we equate the steady state values of the common species of the two subnetworks. That is, since species B appears in both subnetworks, we equate the steady state values of species B in N 1 i.e., k 6 k 5 τ 1 and Fig. 1b bottom). As a result, we get τ 1 = σ 1 k 2 k 5 (σ 1 + k 1 ) k 1 k 3 k 4 k 6 . By using this, we can eliminate the free parameter (τ 1 ) in the combined steady state (Fig.  1c upper left) and derive the analytic steady state of the whole network N (Fig. 1c) with one free parameter σ 1 . This free parameter can be solved using the conservation law b + d = b 0 + d 0 where b 0 and d 0 are the initial concentrations of species B and D, respectively. That is, by substituting the steady state values of species B and D to the conservation law, we can solve for the value of σ 1 in terms of the conserved quantity (b 0 + d 0 ). The closed form can be used to easily identify the critical features of a steady state, such as multistability and concentration robustness, which we will illustrate as follows.

Major comment 4
Please extend the Methods section, e.g. as in Johnston et al (2019) Bull Math Biol, state all results of GCRN theory that you are using, and refer to them in the Results section! From my point of view, details of GCRNs and their theory are "hidden" in this manuscript (to make to presentation look simpler?), but in fact, those details would allow to follow the presentation more easily. To this end, the authors should introduce a bit more terminology, as follows: Given a CRN (with graph G), a (dynamically equivalent) GCRN (G ′ , y,ỹ) is a graph G ′ and two maps y andỹ. It gives rise to two associated CRNs, namely a stoichio-metric one (G, y) and a kinetic-order one (G, y). These three objects are constantly confused in the presentation, and should be referred to properly.
Given the network N used in Figs 1  As stated above, please include in the Methods section all parts of theory you are using: network decomposition, network translation into GCRNs, and GCRN theory.
We thank the reviewer for these helpful and constructive comments.
We have extended the methods section to include all the results of the GCRN theory that we are using as we have responded to Major Comment 3.
Furthermore, we rephrased the paragraph suggested by the reviewer to introduce the stoichiometric and kinetic-order CRN structures of a GCRN to easily define the effective and kinetic deficiencies in Methods as follows: Methods, Lines 417-422: Given a CRN with an associated graph G, a dynamically equivalent GCRN (i.e., the associated ODEs of the GCRN agree with the ODEs of the original CRN) is a graph G ′ together with two maps y andỹ, and hence, with two sets of complexes. It gives rise to two associated CRNs, the stoichiometric one (G ′ , y) and a kinetic-order one (G ′ , y). We call the deficiency associated to the stoichiometric CRN and kinetic-order CRN, the effective deficiency and kinetic deficiency, respectively.
Figs 1 and 2 were revised accordingly to indicate the two structures of the given translated network. We indicated the stoichiometric and kinetic-order CRNs as described above.
Moreover, we modified Fig. 3 to include the stoichiometric and kinetic-order CRNs for the first two subnetworks that were translated. We also provided the proper terms such as effective and kinetic deficiencies in the caption. Furthermore, we clarified the GCRN terminologies in the text as follows: Methods, Line 246: To do this, we decompose the model into seven independent subnetworks (Fig. 3a right). Then, for subnetworks that are not WR and DZ (i.e., N 1 and N 2 ), the network translation was performed to obtain WR and DZ generalized networks. The translated subnetworks N ′ 1 and N ′ 2 are WR and DZ because their stochiometric networks (N ′S 1 and N ′S 2 , respectively) are WR and EDZ networks and kinetic order networks (N ′K 1 and N ′K 2 , respectively) are WR and KDZ networks, while their dynamics are equivalent to the original subnetworks N 1 and N 2 , respectively (Fig. 3b left).
The Methods section now includes the network translation into GCRNs, and important concepts that we used in the GCRN theory (see our response to Major Comment 3). Additionally, we specify the decomposition theorem introduced by Martin Feinberg in the Methods section.
Methods, Lines 396-407: It was shown in [9,10] that when the underlying network decomposition of a reaction network is independent, then the set of positive steady states of the whole system is equal to the intersection of the sets of positive steady states of the subsystems as long as the same kinetics is followed by each reaction from the whole network down to its corresponding subnetwork. The theorem is given as follows (see the Supplementary Information for details).
Theorem 1 Let N be a reaction network with kinetics K decomposed into subnetworks N 1 , N 2 , . . . , N α and K i be the restriction of K to reactions in N i . Then where E is the set of positive steady states of the whole network while E i is the set of positive steady states of subnetwork N i . Furthermore, if the network decomposition is independent, then equality holds, i.e., E 1 ∩ E 2 ∩ · · · ∩ E α = E. and network translation. a A CRISPRi toggle switch network (left) that assumes dCas9 with a single guide RNA 1 (CS complex) and dCas9 with a single guide RNA 2 (CT complex) bind to a single site on their target genes H and G, respectively. The CRN has nine species (i.e., CS, CT , C, G, H, P , R, S, and T ) and 14 reactions (i.e., R 1 , . . . , R 14 ), which is decomposed into seven independent subnetworks (N 1 ,. . . , N 7 ) (right). b As the subnetworks N 1 and N 2 are not WR and DZ, network translation is performed. The translated subnetworks (N ′ 1 and N ′ 2 ) are WR and DZ (i.e., both its stoichiometric and kinetic-order CRNs are WR, and it is an EDZ and a KDZ generalized network) while their dynamics are equivalent to the dynamics of the original subnetworks (N 1 and N 2 , respectively). Then, the steady states of the translated subnetworks (N ′ 1 and N ′ 2 ) and the original subnetworks (N 3 , . . . , N 7 ) can be analytically derived (see Fig. 2 for the outline of the steps). c The steady states of the subnetworks are combined by equating the steady states of the common species. For instance, the species CS is common to both subnetworks N 3 and N 5 so the steady state values of CS of both subnetworks are equated to each other (left). Then after solving the steady state of the whole network, two free parameters (τ 1 and τ 2 ) are left, which can be solved in terms of the conserved quantities. This analytic steady state solution could be used to determine the behavior of the steady state concentrations of species with respect to varying rate constants. In particular, the monotonicity of the steady state concentration of species S over varying rate constant k i , when the rest of the rate constants are set to one, is illustrated on the right.