Structural reduction of chemical reaction networks based on topology
Yuji Hirono, Takashi Okada, Hiroyasu Miyazaki, Yoshimasa Hidaka
KKEK-TH-2296
Structural reduction of chemical reaction networksbased on topology
Yuji Hirono,
1, 2, ∗ Takashi Okada, † Hiroyasu Miyazaki, ‡ and Yoshimasa Hidaka
4, 5, 3, § Asia Pacific Center for Theoretical Physics, Pohang 37673, Korea Department of Physics, POSTECH, Pohang 37673, Korea RIKEN iTHEMS, RIKEN, Wako 351-0198, Japan KEK Theory Center, Tsukuba 305-0801, Japan Graduate University for Advanced Studies (Sokendai), Tsukuba 305-0801, Japan (Dated: February 16, 2021)We develop a model-independent reduction method of open chemical reaction systems basedon the stoichiometry, which determines their network topology. A subnetwork can be eliminatedsystematically to give a reduced system with a fewer number of degrees of freedom. The removalof the subnetwork is accompanied by rewiring of the network, which is prescribed by the Schurcomplement of the stoichiometric matrix. We introduce homology and cohomology groups as acharacterization of the topology of chemical reaction networks, and the change of the topology underthe reduction is traced by the changes of those groups. We prove that, when certain topologicalconditions are met, the steady-state chemical concentrations and reaction rates of the reduced systemis ensured to be the same as those of the original system. The result holds regardless of the modelingof the reactions, namely chemical kinetics, since the conditions involve topological information only.This is advantageous because the details of reaction kinetics and parameter values are difficult toidentify in many practical situations. Our method can reduce a reaction network while preserving itsoriginal steady-state properties, which allows us to study large complex reaction systems efficiently.We demonstrate the reduction method in hypothetical networks and the central carbon metabolismof
Escherichia coli . CONTENTS
I. Introduction 2II. Topology of chemical reaction networks 3A. Chemical reaction networks 4B. Homology, cohomology, and steady states 6C. Subnetworks 9D. Mayer-Vietoris exact sequence 9III. Law of localization 10A. Structural sensitivity analysis 10B. Law of localization 12C. Submodularity of the influence index 13IV. Reduction of chemical reaction networks 15A. Reduction procedure 15B. Simple examples of reduction 17V. Reduction and buffering structures 20A. Decomposition of the influence index 20B. Long exact sequence of a pair of chemical reaction networks 22C. Reduction of buffering structures 24D. Hierarchy of subnetworks 26VI. Example of reduction: metabolic pathway of
E. coli ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] a r X i v : . [ q - b i o . M N ] F e b A. List of reactions 27B. List of buffering structures 29VII. Summary and outlook 30Acknowledgments 31A. Laplace operators and Hodge decomposition 31B. Cycles and conserved charges 321. Interpretation of (cid:101) c ( γ ) and d l ( γ ), and (cid:101) d ( γ ) 322. Embedding of A-matrices 35C. Emergent conserved charges in chemical reaction networks 371. Systems with emergent conserved charges 37a. Pattern A : solutions have arbitrary parameters 37b. Pattern B: no steady-state solution 38c. Pattern C: example with emergent conserved charges and unphysical kinetics 39d. Emergent conserved charges and zero modes of fluctuations 402. Absence of emergent conserved charges in monomolecular reaction networks 40D. Parameter values used in Figure 5 41References 41 I. INTRODUCTION
Chemical reactions in living systems form complex networks and operate in a highly coordinated manner, and theyare responsible for various cellular functions. Experimentally, high-throughput measurements have been conductedto study the response to perturbations for the purpose of elucidating underlying regulatory mechanisms (see, forexample, Refs. [1–3]). One approach to the theoretical studies of biological systems is to build elaborated models,employing particular kinetics, parameter values, and initial/external conditions. Although these models can providedetailed quantitative predictions, a faithful modeling is challenging for most biological systems, because our priorknowledge about kinetics and parameter values is limited and also because many parameter values are difficult tomeasure experimentally. Furthermore, the complexity of models may confound model-independent features withmodel-dependent ones.To address these difficulties, it is desirable to reduce complex reaction systems to simpler ones. A reduction ispractically useful since it can reduce the number of variables and parameters to be dealt with, and can also identifyfeatures essential to focal phenomena or properties of interest (such as biomass production of a metabolic network).It also relates to a conceptual question of the robustness of biochemical processes [4–11]. Chemical reaction networksinside living organisms are highly interconnected, and yet are robust under internal fluctuations and environmentalperturbations. If a system is insensitive to the details of its substructure, it is natural to expect that a reduction ispossible, in the spirit of renormalization. To our best knowledge, the reduction methods [12, 13] of chemical reactionsystems studied so far are, broadly speaking, based on time-scale separation, lumping [14, 15], sensitivity analysis[16–18] or optimizations [19, 20]. Such an approach requires detailed information about the reactions. For example,in order to exploit the time-scale separation, we have to identify fast and slow reactions. The sensitivity analysisalso requires the dependence of the system on various parameters. Accordingly, it is often difficult to distinguish themodel-dependent properties from model-independent ones.In this paper, we develop a systematic method of reducing chemical reaction networks based on their topology (seeFig 1). Our reduction method is motivated by the law of localization [21–23]; if a certain topological index, whichwe call the influence index , is zero for a subnetwork, perturbations inside the subnetwork do not affect the steadystate of the outside (see Sec. III for the precise statement). This observation indicates that certain subnetworksare ‘irrelevant’, as far as their outside is concerned. As we will show, a reduction can be systematically performedthrough the Schur complementation of the stoichiometric matrix with respect to a subset of chemical species andreactions. The well-definedness of the reduction process requires that the subnetwork should satisfy a condition calledthe output-completeness . The behavior of the reduced system depends on the topological nature of the subnetwork.As a central result, we prove that, when the influence index of the subnetwork vanishes, the steady-state chemical e e e e e e e e e e e d x dt = S r ( x ) d x (cid:48) dt = S (cid:48) r (cid:48) ( x (cid:48) ) S (cid:48) := S/S γ v v v v γ Output-completeInfluence index λ ( γ ) = 0 S = v − v − v − − v − − e e e e e e e S γ Reduction v v S (cid:48) = (cid:18) (cid:19) v − v − − e e e e FIG. 1. Schematic of our reduction method. In a chemical reaction network with a stoichiometric matrix S , if a subnetwork γ with S γ satisfies a topological condition called output-completeness, the system can be reduced to a smaller system, whosestoichiometric matrix S (cid:48) is given by the generalized Schur complement S (cid:48) := S/S γ . The reduced system can reproduce thesame steady-state properties of the original system, if the influence index of the subnetwork is zero. Note that S (cid:48) is in generaldifferent from the corresponding submatrix of S (compare the lower-right block of S with S (cid:48) , where the difference betweenthem is indicated by the colored components in S (cid:48) ). This alteration is pictorially represented as ‘rewiring’ of the network (e.g.,the head of e is rewired to v ). concentrations and reaction rates of the reduced system are exactly the same as those of the original system, as faras the remaining degrees of freedom are concerned.We emphasize that those conditions are topological ones and determined only by the network structure, and henceare insensitive to the details of how the reactions are modeled. Thus, the result is broad applicable, because itholds regardless of the kinetics or parameter values. This is of practical merit since the kinetics of reactions or thevalues of parameters are difficult to identify in many situations. For the purpose of characterizing the topology ofreaction networks, we introduce the homology and cohomology groups for chemical reaction networks. The change ofthe topology of chemical reaction networks under the reduction is captured by the change of (co)homology groups,and the tools of algebraic topology are convenient for tracking those changes. We recommend the readers who areinterested in practical aspects of reduction to directly go to Sec. IV where we discuss the reduction procedure withsimple examples.The rest of the paper is organized as follows. In Sec. II, we introduce concepts for characterizing the structureof chemical reaction networks. We introduce the homology and cohomology groups for chemical reaction networks,and the steady-state reaction rates and concentrations are determined by the elements of the cohomology groups. InSec. III, we give a review the structural sensitivity analysis and the law of localization. We also show that the influenceindex is submodular as a function over output-complete subnetworks. In Sec. IV, we introduce the reduction procedureand illustrate the method in simple examples. In Sec. V, we discuss the relation between the structural sensitivityanalysis and the reduction method. We show that the reduction of a buffering structure, that is an output-completesubnetwork with vanishing influence index, has a particularly nice property: The reduced system admits the samesteady states as the original system. In Sec. VI, as an application to realistic networks, we demonstrate the reductionmethod for the metabolic pathway of E. coli . Section VII is devoted to summary and outlook. In Appendix A, wediscuss the Hodge decomposition and Laplace operators for chemical reaction networks. In Appendix B, we giveintuitive interpretations of cycles and conserved charges of various types that appear in the decomposition of theinfluence index. We also illustrate how the decomposition of the index can be seen visually in the structure of theA-matrix, which characterizes the response of the steady state to the perturbations of parameters. In Appendix C,the role of emergent conserved charges in subnetworks is discussed. In Appendix D, we present the parameters usedin the numerical simulations of metabolic pathway of
E. coli.
II. TOPOLOGY OF CHEMICAL REACTION NETWORKS
In this section, we introduce definitions and concepts for characterizing the topology of chemical reaction networks.Those concepts will be used to track the change of reaction networks under reductions.
A. Chemical reaction networks
Definition . A chemical reaction network Γ consists of a quadruple Γ = (
V, E, s, t ),where V is a set of chemical species, E is a set of chemical reactions, and s and t are source and target functions, s : E → N V , t : E → N V , (1)which specify the reactants/products of a reaction. Here, N indicates nonnegative integers, and the elements of N V are maps from V to N .Let us explain the definition in more detail. We will use the indices i, j, k, · · · for chemical species and A, B, C, · · · for chemical reactions. Given a reaction e A ∈ E , we have a map, s ( e A ) : V → N , and s ( e A )( v i ) ∈ N for v i ∈ V indicates how many v i are needed as reactants for the reaction e A . Similarly, t ( e A )( v i ) ∈ N is the number of v i createdin reaction e A . The system can be an open reaction network , since the source or the target function can be zerofor any species (see the example reactions below). When t ( e A )( v i ) = 0 for any v i ∈ V , the product of reaction e A is deposited to the outer world. Similarly, a reaction with s ( e A )( v i ) = 0 for any v i ∈ V is sourced from outside. Areaction is usually represented in the following form, e A : (cid:88) i y iA v i → (cid:88) i ¯ y iA v i , (2)where v i ∈ V , and y iA and ¯ y iA are nonnegative integers. Those integers are given by the source and target functionsas y iA = s ( e A )( v i ) , ¯ y iA = t ( e A )( v i ) . (3)The stoichiometry of the reaction is specified by the stoichiometric matrix S , whose components are given by S iA := ¯ y iA − y iA . (4) Remark . There are several equivalent ways to formulate a chemical reaction network such as a hypergraph [29] ora Petri net [30, 31].
Remark . A reaction that involves at most one chemical species as reactants and products, such as v → v , is called monomolecular . When all the reactions in the system are monomolecular, the corresponding reaction network is ausual directed graph. In this case, the stoichiometric matrix is the incidence matrix of the graph. If we regard Γ as adirected hypergraph, the stoichiometric matrix is the incidence matrix of a directed hypergraph.We consider formal summations of species and reactions with real coefficients, and consider vector spaces whosebases are chemical species/reactions. We denote the resulting vector spaces as C (Γ) := (cid:8) (cid:88) i a i v i | v i ∈ V, a i ∈ R (cid:9) , (5) C (Γ) := (cid:8) (cid:88) A b A e A | e A ∈ E, b A ∈ R (cid:9) . (6)Elements of those spaces are referred to as 0-chains and 1-chains. Higher ( n ≥
2) chains do not exist in the currentsetting. The stoichiometric matrix provides us with natural boundary operators on the space of chains, ∂ n : C n (Γ) → C n − (Γ) . (7)The action of ∂ is defined by its action on the basis e A ∈ C (Γ) and v i ∈ C (Γ), ∂ e A = (cid:88) i ( S T ) Ai v i , ∂ v i = 0 . (8)We often use the notation of linear algebra, where an element (cid:80) i a i v i ∈ C (Γ) is represented by the vector a =( a , a , · · · ) T , and we also write a ∈ C (Γ). For b ∈ C (Γ), the action of the boundary operator is given by themultiplication of the stoichiometric matrix, ∂ : b (cid:55)→ S b ∈ C (Γ) . (9) The compositional aspect of open reaction networks has been studied in the language of category theory [24–26]. Non-equilibriumthermodynamic analysis of open reaction networks with mass-action kinetics and with reversible reactions is performed in Refs. [27, 28].
On the space of chains, let us define inner products by (cid:104) e A , e B (cid:105) = δ AB , (cid:104) v i , v j (cid:105) = δ ij . (10)With these inner products, we can define the adjoint of the boundary operator, ∂ † : C (Γ) → C (Γ) such that (cid:104) ∂ † v i , e A (cid:105) = (cid:104) v i , ∂ e A (cid:105) . The action on the basis v i ∈ C (Γ) is given by ∂ † v i = (cid:88) A S iA e A . (11)In the linear-algebra notation, the action of ∂ † is the multiplication of the transpose of S to a ∈ C (Γ), ∂ † : a (cid:55)→ S T a ∈ C (Γ) . (12) Example . Let us consider a reaction network Γ = ( { v , v , v , v , v } , { e , e , e , e , e , e } ) given by the followingset of chemical reactions, e : (input) → v ,e : (input) → v ,e : v + v → v + v ,e : v → v , (13) e : v → (output) ,e : v → (output) . The stoichiometric matrix of the network is S = − − − − − . (14)It can be drawn as e e e e e v v v v v e (15)We represent a monomolecular reaction by a single arrow, and we use a rectangle to represent a multimolecularreaction. In this network, e is a multimolecular reaction and others are all monomolecular. The action of theboundary operator is, for example, ∂ e = v − v , ∂ e = v + v − v − v , ∂ e = v , (16)and so on. Those are intuitively understood from the figure. The network is open, since we have inputs from theoutside ( e and e ) and outputs to the external world ( e and e ). For example, s ( e )( v i ) = 0 for any v i ∈ V . Theaction of ∂ † is ∂ † v = e − e , ∂ † v = e − e , (17)for example. Namely, the operator ∂ † measures the net inflow of the reactions on a vertex.The chemical concentrations and reaction rates are R -valued linear maps over 0-chains and 1-chains, respectively, C n (Γ) : C n (Γ) → R , (18)for n = 0 ,
1. Given an x ∈ C (Γ), x ( v i ) ∈ R represents the concentration of the chemical species v i . Similarly, for agiven r ∈ C (Γ), r ( e A ) ∈ R represents the rate of the reaction e A . We will also use short-hand notations x i := x ( v i )and r A := r ( e A ). We will also denote an element as a vector as x ∈ C (Γ) and r ∈ C (Γ), where the components of x and r are given by x i and r A , respectively.We define a coboundary operator in a usual way using the boundary operator,(d x )( e A ) := x ( ∂ e A ) = x (cid:32)(cid:88) i ( S T ) Ai v i (cid:33) = (cid:88) i ( S T ) Ai x ( v i ) , (19)where we have used the linearity of the map x . Thus, we can identify the coboundary operator that acts on thechemical concentration x ∈ C (Γ) as the multiplication of the matrix S T .We define the inner product of n -cochains as (cid:104) x, y (cid:105) := (cid:88) i x ( v i ) y ( v i ) , (cid:104) r, s (cid:105) := (cid:88) A r ( e A ) s ( e A ) , (20)where x, y ∈ C (Γ) and r, s ∈ C (Γ). With these inner products, the adjoint of the coboundary operator d n ,d † n : C n +1 (Γ) → C n (Γ), is defined by (cid:104) f n +1 , d n g n (cid:105) n +1 = (cid:104) d † n f n +1 , g n (cid:105) n . (21)Following the definition, we can identify d † as follows, (cid:104) r, d x (cid:105) = (cid:88) A r ( e A ) (d x )( e A )= (cid:88) i,A S iA r ( e A ) x ( v i ):= (cid:88) i (d † r )( v i ) x ( v i )= (cid:104) d † r, x (cid:105) , (22)where r ∈ C (Γ) and x ∈ C (Γ). Thus, the action of d † is given by(d † r )( v i ) := (cid:88) A S iA r ( e A ) . (23)By construction, the adjoint of coboundary operator satisfies (d † r )( v i ) = r ( ∂ † v i ). B. Homology, cohomology, and steady states
With the (co)chains and (co)boundary operators defined above, we can discuss (co)homology groups. We have thefollowing chain complex, 0 (cid:47) (cid:47) C (Γ) ∂ (cid:47) (cid:47) C (Γ) (cid:47) (cid:47) . (24)Noting that the action of ∂ is the multiplication of the stoichiometric matrix S , we can identify the homology groups H (Γ) = C (Γ) /∂ C (Γ) = C (Γ) / im S = coker S, (25) H (Γ) = ker S. (26) Remark . Note that C (Γ) is endowed with a standard inner product, with respect to which we can take theorthogonal linear subspace (im S ) ⊥ . Moreover, the restriction of the quotient map C (Γ) → coker S to (im S ) ⊥ induces an isomorphism (im S ) ⊥ ∼ = −→ coker S . Therefore, we can always regard coker S as a linear subspace of C (Γ).Note also that the orthogonal subspace (im S ) ⊥ is the same as the kernel of the transpose of S , ker S T . Combinedwith the above observation, this implies that we can always identify coker S with ker S T ⊂ C (Γ). More generally, one may define the inner product with a weight function as (cid:104) f, g (cid:105) n := (cid:88) c ∈ C n (Γ) w ( c ) f ( c ) g ( c ) , where w is a R -valued function over C n (Γ). Similarly, with the coboundary operator d , we can define a complex of cochains as0 (cid:47) (cid:47) C (Γ) d (cid:47) (cid:47) C (Γ) (cid:47) (cid:47) . (27)The associated cohomology groups are H (Γ) = { d ∈ C (Γ) | S T d = 0 } = (im S ) ⊥ ∼ = C (Γ) / im S = coker S, (28) H (Γ) = C (Γ) / d C (Γ) = C (Γ) / im S T ∼ = (im S T ) ⊥ = ker S, (29)where ( − ) ⊥ denotes taking orthogonal spaces with respect to the standard inner product on C (Γ) and C (Γ).An Euler number for this complex can be defined as χ (Γ) := | H (Γ) | − | H (Γ) | = | coker S | − | ker S | , (30)where | W | indicates the dimension of the vector space W .Several remarks on the homology and cohomology groups are in order: Remark . Since we consider the R coefficients, the homology and cohomology groups are the same, H n (Γ) ∼ = H n (Γ)for n = 0 , Remark . In the chemistry literature, the elements of H (Γ) are referred to as cycles, and this is consistent with themathematical terminology. Remark . When the network is monomolecular and the corresponding network is a directed graph, the dimension | H (Γ) | is the number of connected components. Remark . Similarly to the homology groups of topological spaces, Laplace operators can be defined and we canperform Hodge decomposition of C (Γ). See Appendix A.The cohomology groups defined above are closely related to the steady states of a reaction network as we seebelow. Let us consider the time evolution of spatially homogeneous chemical concentrations. The change of thechemical concentration is driven by the reactions. The time derivative of the concentration of species v i is given bythe divergence of the reaction rate, ddt x ( v i ) = (d † r )( v i ) , (31)which is more explicitly written as ddt x i ( t ) = (cid:88) A S iA r A . (32)In order to solve the rate equations, we have to specify kinetics of chemical reactions, such as the mass-action kineticsand the Michaelis-Menten kinetics. A reaction’s kinetics gives the reaction rate r A as a function of its substrateconcentrations (i.e., the concentrations of species with y iA >
0) and parameters, r A = r A ( x ; k A ), where k A representsany one of the parameters for the A -th reaction; for example, in the Michaelis-Menten kinetics, k A represents theMichaelis constant or the maximum rate.The elements of H (Γ) and H (Γ) characterize the steady states of chemical reaction networks. The rate equation(32) at the steady state reads (d † r )( v i ) = 0 , or equivalently (cid:88) A S iA r A = 0 , (33)which means that the steady-state reaction rate is an element of the kernel of S , r ∈ ker S ∼ = H (Γ). On the otherhand, the cokernel of S is related to conserved quantities of the system. Given d ∈ coker S ∼ = H (Γ), that satisfies (cid:80) i d i S iA = 0, we have ddt (cid:104) d, x (cid:105) = ddt ( (cid:88) i d i x i ) = (cid:88) i,A d i S iA r A = 0 . (34) Although the natural choice is to consider x and r as the elements of cohomology groups, we can equivalently consider them as elementsof homology groups, since they are isomorphic in the current setting. Thus, the linear combination (cid:80) i d i x i is independent of time and hence is conserved. For this reason, we refer to theelements of coker S as conserved charges . In order to find the steady-state solutions, we have to specify the value ofall the conserved charges. A steady state is specified by an element of H (Γ) and H (Γ), (cid:88) ¯ α (cid:96) ¯ α d ¯ α ∈ H (Γ) , (cid:88) α µ α ( k , (cid:96) ) c α ∈ H (Γ) , (35)where { d ¯ α } and { c α } are basis vectors of H (Γ) and H (Γ), respectively. The coefficients µ α ( k , (cid:96) ) depend on theparameters k and (cid:96) . Example . We consider a network Γ = ( { v , v , v , v } , { e , e , e , e , e } ) with the following reactions, e : (input) → v ,e : v → v ,e : v → (output) , (36) e : v + v → v + v ,e : v + v → v + v . The network structure can be drawn as e e e v v e e v v (37)We here take the mass-action kinetics, and the equations of motion are written as ddt x x x x = − − − − −
10 0 0 1 − r r r r r , r r r r r = k k x k x k x x k x x , (38)where x i = x ( v i ) and r A = r ( e A ) are the concentration and reaction rate for the species v i and reaction e A , respectively.The kernel and cokernel of the stoichiometric matrix are given byker S = span { (cid:0) (cid:1) T , (cid:0) (cid:1) T } , (39)coker S = span { (cid:0) − (cid:1) T } . (40)The cokernel is one-dimensional and the system has one conserved charge. To find the steady states, we need tospecify the value of the charge as (cid:96) = x − x . (41)The steady-state reaction rates and concentrations are¯ r = k (cid:0) (cid:1) T + k k k k (cid:0) (cid:1) T , (42)¯ x = (cid:16) k k k k (cid:0) (cid:96) + √ (cid:96) + 4 K (cid:1) (cid:0) − (cid:96) + √ (cid:96) + 4 K (cid:1)(cid:17) T , (43)where we set K := k k /k k . The vector ¯ r is spanned by the basis vectors of ker S and their coefficients are µ α . “Conserved moiety” may be more chemistry-oriented terminology. C. Subnetworks
Let us consider a subset of chemicals and reactions, γ ⊂ Γ, that we specify by γ = ( V γ , E γ ) with V γ ⊂ V and E γ ⊂ E . Correspondingly, we have a submatrix S γ of the stoichiometric matrix S , whose components are given by( S γ ) iA = S iA , (44)where the indices are restricted to those of the subnetwork, i ∈ V γ , A ∈ E γ . We denote the space of relative chainsby C n ( γ ) := C n (Γ) /C n (Γ \ γ ) , (45)where Γ \ γ := ( V \ V γ , E \ E γ ) is the complement of the subnetwork γ . The homology and cohomology groups forthe subnetwork can be defined similarly. The chain complex for a subnetwork γ is0 (cid:47) (cid:47) C ( γ ) ∂ (cid:47) (cid:47) C ( γ ) (cid:47) (cid:47) , (46)where the action of the boundary operator ∂ on the basis of C ( γ ) is defined with the partial stoichiometric matrix S γ , ∂ e A = (cid:88) i ( S Tγ ) Ai v i . (47)The associated homologies with the complex (46) are H ( γ ) = C ( γ ) /∂ C ( γ ) = C ( γ ) / im S γ = coker S γ , (48) H ( γ ) = ker S γ . (49)The Euler number for a subnetwork is given by χ ( γ ) = | H ( γ ) | − | H ( γ ) | . (50)Note that χ ( γ ) = | H ( γ ) | − | H ( γ ) | = | C ( γ ) | − | C ( γ ) | = | V γ | − | E γ | . (51)The value of the concentrations and reaction fluxes inside a subset γ are given by R -valued functions over the spaceof chemicals and reactions, C n ( γ ) : C n ( γ ) → R . (52)The cohomology for subnetworks can be defined similarly to the homology. D. Mayer-Vietoris exact sequence
In this subsection, we give a long exact sequence of homology groups that connects local and global information.Suppose that there are two subnetworks γ , γ ⊂ Γ, which consist of γ = ( V γ , E γ ) and γ = ( V γ , E γ ). We canconsider the intersection and union of the subnetworks, γ ∩ γ := ( V γ ∩ V γ , E γ ∩ E γ ) , γ ∪ γ := ( V γ ∪ V γ , E γ ∪ E γ ) . (53)The exact sequence (56) below explains the relationship among cohomology groups of γ ∪ γ , γ , γ and γ ∩ γ .Regarding the family { γ , γ } as a “covering” of γ ∪ γ , we can think of Eq. (56) as an analogue of the Mayer-Vietorissequence associated with an open covering of a topological space. Following the usual technique in topology, we willderive the long exact sequence from a short exact sequence of chain complexes. We have the following short exact0sequence of chain complexes, 0 (cid:15) (cid:15) (cid:15) (cid:15) (cid:15) (cid:15) (cid:47) (cid:47) C ( γ ∩ γ ) f (cid:47) (cid:47) ∂ (cid:15) (cid:15) C ( γ ) ⊕ C ( γ ) g (cid:47) (cid:47) ∂ (cid:15) (cid:15) C ( γ ∪ γ ) (cid:47) (cid:47) ∂ (cid:15) (cid:15) (cid:47) (cid:47) C ( γ ∩ γ ) f (cid:47) (cid:47) (cid:15) (cid:15) C ( γ ) ⊕ C ( γ ) g (cid:47) (cid:47) (cid:15) (cid:15) C ( γ ∪ γ ) (cid:47) (cid:47) (cid:15) (cid:15)
00 0 0 (54)where the horizontal maps are given by f n : c (cid:55)→ ( c, − c ) , g n : ( c , c ) (cid:55)→ c + c . (55)By applying the snake lemma to Eq. (54), we obtain0 (cid:47) (cid:47) H ( γ ∩ γ ) (cid:47) (cid:47) H ( γ ) ⊕ H ( γ ) (cid:47) (cid:47) H ( γ ∪ γ ) (cid:47) (cid:47) H ( γ ∩ γ ) (cid:47) (cid:47) H ( γ ) ⊕ H ( γ ) (cid:47) (cid:47) H ( γ ∪ γ ) (cid:47) (cid:47) . (56)In general, if there is an exact sequence of finite-dimensional vector spaces, the alternating sum of the dimensions ofthem is equal to zero. Therefore, the exact sequence (56) implies χ ( γ ∪ γ ) = χ ( γ ) + χ ( γ ) − χ ( γ ∩ γ ) . (57) III. LAW OF LOCALIZATION
A sensitivity analysis studies the response of the system to the perturbations of reaction parameters or initialconditions (conserved charges). In the context of metabolic networks, a theoretical framework called the metaboliccontrol analysis has been developed [32–36]. Under the mass-action framework, biologically insightful results havebeen obtained regarding the sensitivity to conserved charges [8, 11] as well as stability properties of stable states[37–40], although the mass-action law is not necessarily appropriate for some biological systems. Among studies onsensitivity analysis, the structural sensitivity analysis [41–43] aims at constraining the response of reaction systemsfrom the network structure alone.In this section, we first review the structural sensitivity analysis and the law of localization [21–23]. For a givensubnetwork, we assign a nonnegative integer, which we call the influence index. The influence index is determinedfrom the topology of the subnetwork, and plays a decisive role in structural sensitivity. When the influence index iszero, the perturbation of the parameters and conserved charges inside the subnetwork does not affect the rest of thenetwork. Such a structure is called a buffering structure. In Sec. III C, we prove that the influence index is submodularas a function over subnetworks. As a corollary of this property, we show that buffering structures are closed underintersection and union.
A. Structural sensitivity analysis
At the steady state, the reaction rates and the chemical concentrations satisfy (cid:88) A S iA r A ( x ( k , (cid:96) ) , k A ) = 0 , (58) (cid:88) i d ¯ αi x i ( k , (cid:96) ) = (cid:96) ¯ α , (59)1where the second equation specifies the values of conserved charges. Considerable effort has been devoted to the studyof the existence or uniqueness of steady states under the mass-action kinetics [44]. In the current analysis, we assumethe existence of a steady state, and we focus on how it is perturbed under the change of parameters. The steady-statevalues of the concentrations and reaction rates are determined by the values of rate parameters and conserved charges, { k A , (cid:96) ¯ α } . The reaction rates r A ( x ( k , (cid:96) ) , k A ) have explicit dependence on k A , and also dependence on k and (cid:96) through x i ( k , (cid:96) ). Equation (58) means that the reaction rates are in the kernel of S and can be expanded using a basis { c α } of ker S as r A ( x ( k , (cid:96) ) , k A ) = − (cid:88) α µ α ( k , (cid:96) ) c αA . (60)We are interested in the sensitivity of the reaction rates and concentrations under the perturbation of the parameters, k A (cid:55)→ k A + δk A , (cid:96) ¯ α (cid:55)→ (cid:96) ¯ α + δ(cid:96) ¯ α . (61)By taking the derivative of Eqs. (59) and (60) with respect to k B and (cid:96) ¯ β , we obtain the following equations, (cid:88) i ∂r A ∂x i ∂x i ∂k B + ∂r A ∂k B = − (cid:88) α ∂µ α ∂k B c αA , (62) (cid:88) i ∂r∂x i ∂x i ∂(cid:96) ¯ α = − (cid:88) α ∂µ α ∂(cid:96) ¯ α c αA , (63) (cid:88) i d ¯ αi ∂x i ∂k A = 0 , (64) (cid:88) i d ¯ αi ∂x i ∂(cid:96) ¯ β = δ ¯ α ¯ β . (65)Note that r A ( x ( k , (cid:96) ) , k A ) depends explicitly on k A and also depends implicitly on k and (cid:96) through x . The equationscan be compactly written in the matrix form, A (cid:18) ∂ B x i ∂ B µ α (cid:19) = − (cid:18) ∂ B r A (cid:19) , A (cid:18) ∂ ¯ β x i ∂ ¯ β µ α (cid:19) = (cid:18) δ ¯ α ¯ β (cid:19) , (66)where ∂ B := ∂/∂k B , ∂ ¯ β := ∂/∂(cid:96) ¯ β , and we have introduced a partitioned square matrix, A := (cid:18) ∂ i r A c αA d ¯ αi (cid:19) , (67)where the upper-left block is an | A | × | i | matrix whose ( A, i )-th element is given by ∂r A ∂x i evaluated at the steady state,the upper-right one an | A | × | α | matrix consisting of the basis { c α } of ker S , the lower-left one d ¯ αi an | ¯ α | × | i | matrixconsisting of the basis { d ¯ α } of coker S , and the lower-right one the | ¯ α | × | α | zero matrix. Here, we use the notationthat index i for chemicals runs from 1 to | i | . The matrix A is square due to the identity, | i | + | α | = | A | + | ¯ α | . (68)One can see from Eq. (66) that the response to the change of the parameter is determined by the inverse of thematrix A , (cid:18) ∂ B x i ∂ B µ α (cid:19) = − A − (cid:18) ∂ B r A (cid:19) , (cid:18) ∂ ¯ β x i ∂ ¯ β µ α (cid:19) = A − (cid:18) δ ¯ α ¯ β (cid:19) . (69)We refer to A as the A-matrix (“A” indicates that it is an augmented matrix). Its inverse, − A − determines thesensitivity of the system and is called the sensitivity matrix . If we partition A − as A − = (cid:18) ( A − ) iA ( A − ) i ¯ α ( A − ) αA ( A − ) α ¯ α (cid:19) , (70)and noting that ∂ B r A is a diagonal matrix, ∂ B r A ∝ δ BA , the responses of steady-state concentrations and reactionrates (or equivalently, the coefficients µ α in Eq. (60)) to the perturbations of k A and (cid:96) ¯ α are given by ∂ A x i ∝ ( A − ) iA , ∂ ¯ α x i ∝ ( A − ) i ¯ α , ∂ ¯ α µ α ∝ ( A − ) α ¯ α , ∂ A µ α ∝ ( A − ) αA . (71)In this paper, we consider the following class of chemical reaction systems:2 Definition . A chemical reaction network with kinetics iscalled regular , if it admits a stable steady state and the associated A -matrix is invertible.Note that whether a reaction system is regular or not depends on the choice of kinetics. Throughout the paper, weassume the regularity unless otherwise stated so that A is invertible and the response of the system is well-defined. Theregularity implies the asymptotic stability of the steady state, through the relation between det A and the determinantof the Jacobian [23]. B. Law of localization
Definition . When a subnetwork γ = ( V γ , E γ ) satisfies the condition that E γ includes all thechemical species affected by V γ , γ is called output-complete . Definition . For an output-complete subnetwork γ , the influence index is defined by λ ( γ ) := −| V γ | + | E γ | − | (ker S ) supp γ | + | P γ (coker S ) | . (72)The definitions of the spaces that appear in the influence index are given as follows:(ker S ) supp γ := (cid:8) c (cid:12)(cid:12) c ∈ ker S, P γ c = c (cid:9) = ker S ∩ C ( γ ) , (73) P γ (coker S ) := (cid:8) P γ d (cid:12)(cid:12) d ∈ coker S (cid:9) , (74)where S is the stoichiometric matrix, P γ and P γ are the projection matrices to γ in the space of chemical speciesand reactions, respectively. Namely, (ker S ) supp γ is the space of vectors of ker S supported inside γ , and P γ (coker S )is the projection of coker S to γ . Here, recall from Remark 3 that we regard coker S as a subspace of C (Γ) via theidentification coker S ∼ = (im S ) ⊥ . We will use similar identifications throughout this paper. Remark . The influence index is nonnegative, λ ( γ ) ≥
0, for a regular chemical reaction network. It will be shown inthe proof of Theorem 1.
Theorem . Let γ be an output-complete subnetwork of a regular chemical reaction network Γ.When γ is a buffering structure, λ ( γ ) = 0, chemical concentrations and reaction rates outside γ do not change underthe perturbation of rate parameters or conserved charges inside γ . Definition . For a given chemical reaction network Γ, an output-complete subnetwork γ withthe vanishing influence index, λ ( γ ) = 0, is called a buffering structure . Example . The influence index of the empty subnetwork is zero. The index of the whole network Γ is also zero, λ (Γ) = −| C (Γ) | + | C (Γ) | − | ker S | + | coker S | = | H (Γ) | − | H (Γ) | − ( | C (Γ) | − | C (Γ) | )= 0 . (75)This is natural in the sense that there is no “outside” of the whole network. Example . Let us take the same network as Example 2. The A-matrix of this system is A = r , r , r , r , r , r , − , (76)where r A,i := ∂r A /∂x i and it is evaluated at the steady state. With this matrix A , the responses of the concentrationand reaction rates to the change of parameters and the value of conserved charges can be obtained by Eq. (69). Thesubnetwork γ = ( { v , v } , { e } ) is output-complete and is a buffering structure, since λ ( γ ) = − − γ = ( { v , v } , { e , e } ) is also a buffering structure, λ ( γ ) = − − γ . This explains the fact that x and x do not depend on the value of conserved charge (cid:96) = x − x . The subnetwork γ = ( { v , v , v } , { e , e , e } ) is also a buffering structure, λ ( γ ) = − − x does not depend on k , k , k , and (cid:96) .3 d | E γ | c | V γ | cycles of supported in Γ γ Projected conserved charges A = chemical species in γ reactions in γ FIG. 2. Structure of the A-matrix. The numbers c and d are given by c = | (ker S ) supp γ | and d = | P γ (coker S ) | . Proof.
The law of localization follows from the structure of the matrix A . Given an output-complete subnetwork γ , we can bring the rows and columns associated with γ in the way shown in Fig. 2. All the component of thelower-left part is zero, because the reaction rate r A outside γ does not depend on the chemical species in γ (since γ is output-complete), and those cycles are supported in γ . The index λ ( γ ) measures how far the black rectangle onthe upper-left corner is from a square matrix. The numbers c and d in Fig. 2 are given by c = | (ker S ) supp γ | and d = | P γ (coker S ) | , which appear in Eq. (72). Because of the assumption of regularity, we have det A (cid:54) = 0, and theblack rectangle on the upper-left corner should be vertically long (if it is horizontally long, the determinant vanishes),which is equivalent to the condition λ ( γ ) ≥ λ ( γ ) = 0, the black box in the upper-left corner is a square matrix. Then, A − inherits the same structure, A − = (cid:18) ∗ ∗ ∗ (cid:19) . (77)Namely, if we denote the generic index of ( A − ) as µ, ν, · · · and write the index inside and outside γ as µ (cid:63) and µ (cid:48) ,respectively, we have ( A − ) µ (cid:48) ν (cid:63) = 0. Because of this structure, ∂ A (cid:63) x i (cid:48) ∝ ( A − ) i (cid:48) A (cid:63) = 0 , (78)which means that the concentrations out of γ do not depend on the parameter k B inside γ . Consequently, we have ∂ A (cid:63) r A (cid:48) ∝ (cid:88) i (cid:48) ∂ i (cid:48) r A (cid:48) ∂ A (cid:63) x i (cid:48) = 0 , (79)where we used the fact that r A (cid:48) only depends on the concentrations outside γ because of the output-completeness.The same is true for the perturbation of the conserved charge, ∂ ¯ α (cid:63) x i (cid:48) ∝ ( A − ) i (cid:48) ¯ α (cid:63) = 0 , ∂ ¯ α (cid:63) r A (cid:48) ∝ (cid:88) i (cid:48) ∂ i (cid:48) r A (cid:48) ∂ ¯ α (cid:63) x i (cid:48) = 0 . (80) C. Submodularity of the influence index
The influence index λ ( γ ) can be regarded as a function over subnetworks. We here show that the influence indexsatisfies an inequality. As a corollary, we show that the buffering structures are closed under union and intersection.This fact is useful in enumerating buffering structures in large reaction networks.We first note that:4 • Given output-complete subnetworks γ , γ ⊂ Γ, the union and intersection, γ ∪ γ and γ ∩ γ , are also output-complete. This follows from the definition of output-completeness. • A function f ( γ ) over a set is called submodular , when it satisfies f ( γ ∪ γ ) ≤ f ( γ ) + f ( γ ) − f ( γ ∩ γ ) . (81)When ≤ is replaced with ≥ , the function satisfying the replaced equation is called supermodular . Theorem . Let γ , γ ⊂ Γ be output-complete subnetworks. The influence index satisfies λ ( γ ∪ γ ) ≤ λ ( γ ) + λ ( γ ) − λ ( γ ∩ γ ) . (82)Namely, λ ( γ ) is a submodular function over output-complete subnetworks. Proof.
We show that λ ( γ ) = −| V γ | + | E γ | − | (ker S ) supp γ | + | P γ (coker S ) | (83)is submodular. Recall that χ ( γ ) = | V γ | − | E γ | = | H ( γ ) | − | H ( γ ) | is the Euler number for subnetwork γ = ( V γ , E γ ).We note that χ ( γ ) is a modular function, meaning that it satisfies χ ( γ ∪ γ ) = χ ( γ ) + χ ( γ ) − χ ( γ ∩ γ ) , (84)which is derived from the Mayer-Vietoris exact sequence (56). Thus, it suffices to show that the last two terms onthe RHS of Eq. (83) are submodular. In fact, we show that each of them is submodular.Let us first look at | P γ (coker S ) | . If denote W := coker S , the submodularity of | P γ (coker S ) | reads | P γ ∪ γ W | ≤ | P γ W | + | P γ W | − | P γ ∩ γ W | . (85)We prove this equation just after this proof. Thus, we have shown the submodularity of | P γ (coker S ) | .Next, we show that | (ker S ) supp γ | is supermodular. Consider the following vector space, Z := (ker S ) supp γ + (ker S ) supp γ . (86)Its dimension is given by | Z | = | (ker S ) supp γ | + | (ker S ) supp γ | − | (ker S ) supp γ ∩ (ker S ) supp γ | = | (ker S ) supp γ | + | (ker S ) supp γ | − | (ker S ) supp γ ∩ γ | . (87)Since any element of Z is supported in γ ∪ γ , we have (ker S ) supp γ ∪ γ ⊃ Z , which implies | (ker S ) supp γ ∪ γ | ≥ | Z | .Thus, we have shown that | (ker S ) supp γ | is a supermodular function, and −| (ker S ) supp γ | is submodular.Therefore, | P γ (coker S ) | and −| (ker S ) supp γ | are both submodular function, and we obtain the claim. Proof of Eq. (85) . Let us pick a basis of the space W as { w , · · · , w n } and we denote the basis as a matrix, w ij .The set of possible row and column indices are denoted as V and N , respectively. For a subset of indices v ⊂ V , wedenote the corresponding submatrix of w ij as w [ v, N ]. With this notation, the dimension of a projected subspace of W is written as | P γ W | = rank w [ v γ , N ] for a subnetwork γ = ( v γ , e γ ).Let us pick two subnetworks γ and γ , and denote the sets of chemical species by v and v , respectively. Weconsider the submatrix w [ v ∩ v , N ]. We can pick a row basis as { a Ti | i ∈ α ∩ } , where α ∩ ⊂ v ∩ v . Here, | α ∩ | = rank w [ v ∩ v , N ]. We can form a row basis of w [ v , N ] by adding row vectors from w [ v \ v , N ] to α ∩ .Let α be the picked indices, then | α ∩ | + | α | = rank w [ v , N ] . (88)We can further pick row vectors from w [ v \ v , N ] and form a basis of w [ v ∪ v , N ]. Let us denote the added indicesas α , then | α ∩ | + | α | + | α | = rank w [ v ∪ v , N ] . (89)Since the vectors specified by the indices α ∩ ∪ α are linearly independent and α ∩ ∪ α ⊂ v , we have | α ∩ | + | α | ≤ rank w [ v , N ]. This can be written asrank w [ v ∪ v , N ] ≤ rank w [ v , N ] + rank w [ v , N ] − rank w [ v ∩ v , N ] . (90)This is equivalent to Eq. (85).5 (cid:3) Corollary . Let Γ be a regular chemical reaction network. The union and the intersection of two buffering structuresinside Γ are also buffering structures.
Proof.
Suppose γ and γ are buffering structures inside Γ. Then λ ( γ ) = λ ( γ ) = 0. From the submodularity of theinfluence index, we have λ ( γ ∪ γ ) + λ ( γ ∩ γ ) ≤ . (91)Since influence indices are nonnegative for a regular chemical reaction network, we have λ ( γ ∪ γ ) = λ ( γ ∩ γ ) = 0.Thus, we obtain the claim. IV. REDUCTION OF CHEMICAL REACTION NETWORKS
Generically, a reduction is a process to reduce the number of degrees of freedom while keeping some features of theoriginal system. A reduction of chemical reaction networks involves two ingredients:(1) Identify a subnetwork to be reduced.(2) Perform the reduction for given a subnetwork.As a result, we obtain a new reaction network with fewer chemical species and reactions,Γ −→ Γ (cid:48) . (92)The reduced network is characterized by a new stoichiometric matrix, S −→ S (cid:48) . (93)Crucial points are how to identify a subnetwork elimination of which does not affect the chemical functionalities ofthe original network, and how to obtain the new stoichiometric matrix S (cid:48) . In this section, we introduce a procedureof the reduction of chemical reaction networks based on the network topology, and mainly discuss step (2). We willdiscuss more on the choice of a subnetwork in later sections. A. Reduction procedure
Let us here illustrate a method of reduction based on the network topology. We denote the whole reaction networkby Γ = (
V, E ), where V and E are the sets of chemical species and reactions, respectively. . We choose a subnetwork γ = ( V γ , E γ ), where V γ ⊂ V and E γ ⊂ E , and eliminate the degrees of freedom inside γ . We refer to the chemicalspecies and reactions inside γ as internal , and those in Γ \ γ as boundary . For the given subnetwork γ , we separatethe chemical concentrations and reaction rates as x = (cid:18) x x (cid:19) , r = (cid:18) r r (cid:19) , (94)where 1 and 2 correspond to internal and boundary degrees of freedom, respectively. The stoichiometric matrix S can also be partitioned into four blocks as S = (cid:18) S S S S (cid:19) . (95)Note that the submatrix S is the same matrix as S γ that appeared in Sec. II C. Hereafter we use S for notationalconvenience. The rate equation of the whole reaction system reads ddt (cid:18) x x (cid:19) = (cid:18) S S S S (cid:19) (cid:18) r r (cid:19) = (cid:18) S r + S r S r + S r (cid:19) . (96)6While the internal reaction rates r = r ( x , x ) in general depend on both of the internal and boundary chemicalconcentrations, when γ is chosen to be output-complete, the boundary reaction rates r = r ( x ) do not depend onthe internal chemical concentrations x . We can solve the first equation of Eq. (96) for r as r = S +11 ddt x − S +11 S r + c , (97)where S +11 is the Moore-Penrose inverse of S , and c ∈ ker S . Substituting this to the second equation of Eq. (96),we get ddt (cid:0) x − S S +11 x (cid:1) = ( S − S S +11 S ) r + S c . (98)When the following condition is satisfied, ker S ⊂ ker S , (99) S c = and the second term of the RHS of Eq. (98) vanishes. The rate equation is now written as ddt (cid:0) x − S S +11 x (cid:1) = S (cid:48) r . (100)Here, S (cid:48) is the generalized Schur complement, S (cid:48) = S/S := S − S S +11 S . (101)As long as steady states are concerned, the subnetwork ( x , r ) satisfies the rate equation whose stoichiometric matrixis S (cid:48) . This motivates us to consider the subnetwork ( x , r ) whose rate equation is given by ddt x = S (cid:48) r ( x ) . (102)Based on the considerations above, we define the reduction as follow: Definition . Let Γ = (
V, E ) be a chemical reaction network with stoichiometric matrix S and γ =( V γ , E γ ) be an output-complete subnetwork whose stoichiometric matrix is denoted by S . We define a reducednetwork Γ (cid:48) = ( V \ V γ , E \ E γ ) obtained by eliminating γ from Γ, by a stoichiometric matrix S (cid:48) given by the generalizedSchur complement, S (cid:48) = S/S := S − S S +11 S . (103)The second term represents the rewiring of the network associated with the elimination of γ . We denote the resultantreaction network by Γ (cid:48) = Γ /γ . Accordingly, the chemical concentrations and reaction rates of the reduced system( x (cid:48) , r (cid:48) ) are obtained from the original ones ( x , r ) as x = (cid:18) x x (cid:19) −→ x (cid:48) = x , (104) r ( x ) = (cid:18) r ( x , x ) r ( x ) (cid:19) −→ r (cid:48) ( x (cid:48) ) = r ( x ) , (105)and the rate equation of the reduced system is given by ddt x (cid:48) = S (cid:48) r (cid:48) ( x (cid:48) ) . (106) Remark . The reduced system can be always defined if γ is output-complete; otherwise, the reduction is ill-definedsince the reduced system would depend on x through r . We emphasize that the output-completeness is a topologicalcondition determined by the stoichiometry and the details of the reactions, namely the kinetics, are irrelevant. Thus,the reduction is applicable to any kind of kinetics. How the reduced system is related to the original system dependson further nature of γ . In the following sections, we will discuss more on the features of the subnetworks that behavenicely under reductions. In Sec. V C, we prove that, when γ has a vanishing influence index (see Sec. III), which isdetermined by the network topology, the steady state of the reduced system is assured to be the same as the steadystate of the original system. As we discuss later, this condition is the same as the absence of emergent cycles in γ . See the text around Eq. (146). Remark . We note that the elements of the matrix S (cid:48) are rational, since the Moore-Penrose inverse of an integralmatrix is rational [45]. The matrix S (cid:48) can be always transformed into an integral matrix by column-wise rescaling of S (cid:48) together with the rescaling of the reaction rates. Remark . The current method is different from the ones discussed for reaction networks with the mass-action kineticsin detailed balanced [46] and complex balanced [47] situations, where the Schur complementation is performed forthe weighted Laplacian similarly to the Kron reduction of electrical circuits [48–50]. In the current formulation, theSchur complementation is performed for the stoichiometric matrix.
B. Simple examples of reduction
To illustrate the reduction procedure, here we discuss simple examples of reductions. In Sec. VI, we discuss thereduction of the metabolic pathway of
E. coli as a more realistic example.
Example . We consider a monomolecular reaction network that consists of (
V, E ) = ( { v , v } , { e , e , e } ). We takea subnetwork γ = ( { v } , { e } ) to be reduced. Under the reduction, the stoichiometric matrix changes as S = (cid:18) (cid:19) v − v − e e e S S (cid:48) = ( ) v − e e (107)where we have brought the reduced part to the upper-left part. The reduction looks like e e e e e v v γ v (108)The original rate equation is ddt (cid:18) x x (cid:19) = (cid:18) − − (cid:19) r ( x ) r r ( x ) , (109)where x = x ( v ), r = r ( e ) and so on. If we eliminate r ( x ), ddt ( x + x ) = (cid:0) − (cid:1) (cid:18) r r ( x ) (cid:19) . (110)The reduced equation of motion is obtained by replacing x + x with x on the LHS.In order to compute the steady-state solutions, let us for example employ the mass-action kinetics, r ( x ) r r ( x ) = k x k k x . (111)The steady-state reaction rates and concentrations are given by ¯ r ¯ r ¯ r = k , (cid:18) ¯ x ¯ x (cid:19) = (cid:18) k /k k /k (cid:19) . (112)The steady-state reaction solutions of the reduced system are (cid:18) ¯ r ¯ r (cid:19) = k (cid:18) (cid:19) , ¯ x = k k . (113)8Note that this solution of the reduced system is exactly the same as the solution (112) of the original system for theboundary concentrations and rates. Indeed, this is a special property of buffering structures. In this example, thesubnetwork has a vanishing influence index, λ ( γ ) = − − Example . ( V, E ) = ( { v , v , v } , { e , e , e } ). The stoichiometric matrices of the original and reduced system aregiven by S = (cid:32) (cid:33) v − v − v − e e e S (cid:48) = (cid:18) (cid:19) v − v − e e S where we reduced the subnetwork γ = ( { v } , { e } ). The reduction is visually expressed as e e e e e v v v γ v v (114)Suppose that we take the mass-action kinetics, r ( x ) r ( x ) r ( x ) = k x k x k x . (115)The system has one conserved charge and we specify the value as (cid:96) = x + x + x . The steady-state reaction ratesfor the original system are ¯ r ¯ r ¯ r = (cid:96)K , (116)where K is defined by K := k + k + k . In the reduced system, (cid:96) (cid:48) = x + x is a conserved charge. The steady-staterates in the reduced system are (cid:18) ¯ r ¯ r (cid:19) = (cid:96) (cid:48) K (cid:48) (cid:18) (cid:19) , (117)where K (cid:48) := k + k . Note that, if we want to have the same steady-state in the reduced system as the one in theoriginal system, we have to choose the parameters so that (cid:96)K = (cid:96) (cid:48) K (cid:48) . This is in contrast to Example 5, whereno fine-tuning of the parameters is needed. The difference is attributed to the fact that the subnetwork γ is not abuffering structure and the index is nonzero, λ ( γ ) = − − Example . ( V, E ) = ( { v , v , v , v } , { e , e , e , e , e , e } ). The stoichiometric matrix changes under reduction as S = v − v − v − − v − e e e e e e S S (cid:48) = (cid:18) (cid:19) v − − v − e e e e (118)9where we chose the subnetwork γ = ( { v , v } , { e , e } ) to be reduced. The reduction is visually expressed as e e e e e e e e e e v v v v γ v v (119)The subnetwork is a buffering structure, λ ( γ ) = − − Example . ( V, E ) = ( { v , . . . , v } , { e , . . . , e } ) with the following stoichiometric matrix, S = v − − v − v − v − − v − v − v − v − − − v − e e e e e e e e e e e e e e S .We choose the subnetwork γ = ( { v , v , v , v , v } , { e , e , e , e , e , e , e } ) to be reduced. The reduced subnetworkis given by S (cid:48) = v − v − v − − v − e e e e e e e .The subnetwork γ is a buffering structure: λ ( γ ) = − − e e e e e e e e e e e e e e e e e e e e v v v v v v v v e v γ v v v v (120)0We note that, under the reduction, the stoichiometries for reactions e , e , e are changed from the original ones.In particular, e , which is originally monomolecular, becomes non-monomolecular, 2 v → v . To reproduce steadystates of the original system, the rate r is required to be the same as before the reduction; for example, in themass-action kinetics, r is given by r ( x ) = k x rather than by r ( x ) = k x , even after the reduction. V. REDUCTION AND BUFFERING STRUCTURES
There is a close connection between the structural sensitivity analysis and the reduction method we introducedin the previous section. The structural sensitivity analysis works as a guide to identify ’unimportant’ subnetworks.In this section, we present a key result of this paper: we show that, when a subnetwork is a buffering structure,the reduced network has exactly the same steady-state solution as the original reaction network. The proof will becompleted in subsection C.The structure of this section is as follows: In subsection A, we show that the influence index allows for a decompo-sition in terms of the numbers of cycles and conserved charges. In subsection B, we construct a short exact sequenceof the chain complexes for a subnetwork γ ⊂ Γ, under some conditions. This short exact sequence automaticallyderives a long exact sequence of homology groups. Using this exact sequence, we can describe the relationship amongcycles and conserved charges of γ , Γ and Γ (cid:48) . In subsection C, we show the main result, that is, that the steady stateof the reduced network is the same as the one of the original network, under some conditions. In the proof, the longexact sequence prepared in subsection B plays an important role. In subsection D, we study the situation where wehave nested subnetworks γ (cid:48) ⊂ γ ⊂ Γ. In this case, we have a subnetwork γ/γ (cid:48) ⊂ Γ /γ (cid:48) . We will show that the reducednetwork Γ /γ is the same as (Γ /γ (cid:48) ) / ( γ/γ (cid:48) ). This ensures that the eventual network does not depend on the orderingof the reductions. A. Decomposition of the influence index
As we detailed in Sec. II, steady-state properties are captured by the cycles and conserved charges, which are theelements of homology groups. In this subsection, we study their meaning in more detail, and discuss the relationbetween the influence index λ ( γ ) and cycles/conserved charges in γ , Γ, and Γ (cid:48) . We introduce a decomposition of theinfluence index in terms of the spaces of cycles/conserved charges of certain classes.We first note that the index can be written as λ ( γ ) = −| V γ | + | E γ | − | (ker S ) supp γ | + | P γ (coker S ) | = | ker S | − | (ker S ) supp γ | + | P γ (coker S ) | − | coker S | , (121)where we used Eq. (51). With the first two terms, we define (cid:101) c ( γ ) := | ker S | − | (ker S ) supp γ | = | ker S | − | ker S ∩ C ( γ ) | , (122)where we regard C ( γ ) as a subspace of C (Γ) through the inclusion map. The number (cid:101) c ( γ ) is a nonnegative integer,because there is an injective map from (ker S ) supp γ to ker S . Indeed, an element of (ker S ) supp γ is written as c = (cid:18) c (cid:19) satisfying the condition (cid:18) S S S S (cid:19) (cid:18) c (cid:19) = (cid:18) S c S c (cid:19) = . (123)Consider an injective map c (cid:55)→ c . Equation (123) indicates that the image of this map is always included in ker S ,(ker S ) supp γ (cid:51) c (cid:55)→ c ∈ ker S . (124)Thus, we have (cid:101) c ( γ ) ≥ | P γ (coker S ) | = | coker S | − | im ¯ P γ ∩ coker S | , (126) For a vector space V and a projection matrix P , | P V | = |{ P v | v ∈ V }| = | V/ (im ¯ P ∩ V ) | = | V | − | im ¯ P ∩ V | , (125)where ¯ P := 1 − P . P γ := 1 − P γ . The second term of the RHS of Eq. (126) is the number of conserved charges of Γ with supporton Γ \ γ , ¯ d (cid:48) ( γ ) := | im ¯ P γ ∩ coker S | = | ¯ D (cid:48) ( γ ) | , (127)where the space ¯ D (cid:48) ( γ ) is given by ¯ D (cid:48) ( γ ) := (cid:26)(cid:18) d d (cid:19) ∈ coker S (cid:12)(cid:12)(cid:12)(cid:12) d = (cid:27) . (128)We divide the space coker S according to the following distinctions: • Projection to γ is also a conserved charge in γ . • Projection to γ is not a conserved charge in γ .Correspondingly to the two distinctions above, we introduce the following spaces, D ( γ ) := X ( γ ) / ¯ D (cid:48) ( γ ) , (129) D (cid:48) ( γ ) := coker S/D ( γ ) ∼ = (coker S ) /X ( γ ) ⊕ ¯ D (cid:48) ( γ ) , (130)where we defined X ( γ ) := (cid:26)(cid:18) d d (cid:19) ∈ coker S (cid:12)(cid:12)(cid:12)(cid:12) d ∈ coker S (cid:27) . (131)The dimension of coker S is written as | coker S | = d ( γ ) + d (cid:48) ( γ ) , (132)where d ( γ ) := | D ( γ ) | , and d (cid:48) ( γ ) := | D (cid:48) ( γ ) | . We now have the expression | P γ (coker S ) | = d ( γ ) + d (cid:48) ( γ ) − ¯ d (cid:48) ( γ ) . (133)In order to rewrite | coker S | , we introduce the following spaces, D ( γ ) := (cid:26) d ∈ coker S (cid:12)(cid:12)(cid:12)(cid:12) ∃ d such that (cid:18) d d (cid:19) ∈ coker S (cid:27) , (134) (cid:101) D ( γ ) := coker S /D ( γ ) . (135)The element of D ( γ ) are conserved charges in γ that can be extended to a global conserved charge, while those in (cid:101) D ( γ ) are emergent conserved charges that are only conserved in the subnetwork γ .Observe that D ( γ ) ∼ = D ( γ ). Indeed, there is a surjection X ( γ ) (cid:51) (cid:18) d d (cid:19) → d ∈ D ( γ ) . (136)The kernel of this map is ¯ D (cid:48) ( γ ), and the induced map D ( γ ) = X ( γ ) / ¯ D (cid:48) ( γ ) → D ( γ ) is an isomorphism. Thus, | D ( γ ) | = | D ( γ ) | = d ( γ ) and we have the decomposition, | coker S | = d ( γ ) + (cid:101) d ( γ ) , (137)where (cid:101) d ( γ ) := | (cid:101) D ( γ ) | is the number of charges that cannot be obtained as the projections of conserved charges in Γ.Combining Eqs. (122), (133), and (137), we find that the influence index is written as λ ( γ ) = (cid:101) c ( γ ) + d l ( γ ) − (cid:101) d ( γ ) , (138)where we defined d l ( γ ) := d (cid:48) ( γ ) − ¯ d (cid:48) ( γ ) = | (coker S ) /X ( γ ) | . (139)The decomposition (138) is the central result of this subsection. Each term of Eq. (138) allows for the followingintuitive interpretations: Note that we can regard the element of D ( γ ) as a vector in coker S by the isomorphism X ( γ ) / ¯ D (cid:48) ( γ ) ∼ = X ( γ ) ∩ [ ¯ D (cid:48) ( γ )] ⊥ . The isomorphismin Eq. (130) can be derived as follows:(coker S ) /D ( γ ) = (coker S ) / (cid:16) X ( γ ) ∩ [ ¯ D (cid:48) ( γ )] ⊥ (cid:17) ∼ = (coker S ) ∩ (cid:16) X ( γ ) ∩ [ ¯ D (cid:48) ( γ )] ⊥ (cid:17) ⊥ = (coker S ) ∩ (cid:16) [ X ( γ )] ⊥ + ¯ D (cid:48) ( γ ) (cid:17) = (coker S ) /X ( γ ) ⊕ ¯ D (cid:48) ( γ ) , where we used the relations V/W ∼ = V ∩ W ⊥ , ( V ∩ W ) ⊥ = V ⊥ + W ⊥ for vector spaces W ⊂ V , and [ X ( γ )] ⊥ + ¯ D (cid:48) ( γ ) = [ X ( γ )] ⊥ ⊕ ¯ D (cid:48) ( γ )since [ X ( γ )] ⊥ ∩ ¯ D (cid:48) ( γ ) = . • The first term (cid:101) c ( γ ) = | ker S / (ker S ) supp γ | represents the number of emergent cycles in γ . Namely, (cid:101) c ( γ ) is thenumber of cycles in γ , which are not cycles in Γ. • The second term d l ( γ ) = | (coker S ) /X ( γ ) | is the dimension of the space of lost conserved charges by focusingon γ , namely those that are conserved in Γ but their projection to γ are not. • The third term (cid:101) d ( γ ) = | (cid:101) D ( γ ) | is the number of emergent conserved charges in γ . It is the number of conservedcharge in γ that cannot be extended to conserved charges in Γ. The meaning becomes evident if we note that (cid:101) D ( γ ) is isomorphic to the space that consists of d ∈ coker S that are orthogonal to the vectors that can beextended to conserved charges in coker S .For more detailed explanations with examples, see Appendix B 1. In Appendix B 2, we show that the decomposition(138) can be visually understood from the structure of A-matrices. B. Long exact sequence of a pair of chemical reaction networks
The reduction of a reaction network naturally induces the reduction of (co)homology groups, which are the steady-state characteristics of reaction networks. Suppose that we have a reaction network Γ, and choose a subnetwork γ ⊂ Γ, and reduce it to obtain Γ (cid:48) = Γ /γ . The inter-relations of homologies of γ , Γ, and Γ (cid:48) , can be systematicallydiscussed using a long exact sequence for a pair of chemical reaction networks, which we define momentarily. Weconsider the following short exact sequence of chain complexes,0 (cid:15) (cid:15) (cid:15) (cid:15) (cid:15) (cid:15) (cid:47) (cid:47) C ( γ ) ψ (cid:47) (cid:47) ∂ γ (cid:15) (cid:15) C (Γ) ϕ (cid:47) (cid:47) ∂ (cid:15) (cid:15) C (Γ (cid:48) ) (cid:47) (cid:47) ∂ (cid:48) (cid:15) (cid:15) (cid:47) (cid:47) C ( γ ) ψ (cid:47) (cid:47) (cid:15) (cid:15) C (Γ) ϕ (cid:47) (cid:47) (cid:15) (cid:15) C (Γ (cid:48) ) (cid:47) (cid:47) (cid:15) (cid:15)
00 0 0 (140)where the space of chains in Γ (cid:48) is given by C n (Γ (cid:48) ) := C n (Γ) /C n ( γ ). In the linear-algebra notations, the boundarymaps are given by the following multiplications of matrices on vectors, ∂ γ : c (cid:55)→ S c , ∂ : c = (cid:18) c c (cid:19) (cid:55)→ S c , ∂ (cid:48) : c (cid:55)→ S (cid:48) c . (141)We define the horizontal maps by ψ : c (cid:55)→ (cid:18) c (cid:19) , ϕ : (cid:18) c c (cid:19) (cid:55)→ c , (142) ψ : d (cid:55)→ (cid:18) d S S +11 d (cid:19) , ϕ : (cid:18) d d (cid:19) (cid:55)→ d − S S +11 d . (143)The exactness of the rows of Eq. (140) can be checked easily. One can see that the diagram (140) commutes whenthe following condition is satisfied: S (1 − S +11 S ) c = , (144)where c ∈ C ( γ ). The matrix (1 − S +11 S ) is the projection matrix to ker S , and Eq. (144) is equivalent toker S ⊂ ker S . (145)This condition is the same as the condition that an arbitrary term in Eq. (98) vanishes.3The condition (145) has a natural interpretation in terms of cycles: Eq. (145) is equivalent to (cid:101) c ( γ ) = 0, namely theabsence of emergent cycles, which can be checked as follows. When (cid:101) c ( γ ) = | ker S / (ker S ) γ | = 0, any c ∈ ker S isa cycle in Γ by an inclusion to C (Γ). Thus, c satisfies (cid:18) S S S S (cid:19) (cid:18) c (cid:19) = . (146)This implies S c = and we have ker S ⊂ ker S . Conversely, when ker S ⊂ ker S is true, the map ker S (cid:51) c (cid:55)→ (cid:18) c (cid:19) ∈ (ker S ) supp γ is a bijection. This implies (cid:101) c ( γ ) = 0. Thus, we have shown that the diagram (140)commutes if and only if γ has no emergent cycle.Applying the snake lemma to Eq. (140), we obtain a long exact sequence,0 (cid:47) (cid:47) H ( γ ) ψ (cid:47) (cid:47) H (Γ) ϕ (cid:47) (cid:47) H (Γ (cid:48) ) δ (cid:47) (cid:47) H ( γ ) ¯ ψ (cid:47) (cid:47) H (Γ) ¯ ϕ (cid:47) (cid:47) H (Γ (cid:48) ) (cid:47) (cid:47) , (147)where ¯ ψ and ¯ ϕ are induced maps of ψ and ϕ . The map δ : H (Γ (cid:48) ) → H ( γ ) is called the connecting map. For agiven c ∈ H (Γ (cid:48) ), the connecting map is given by δ : c (cid:55)→ [ S c ] ∈ H ( γ ) = coker S , (148)where [ ... ] means to identify the differences in im S .Let us look at the consequences of the long exact sequence (147). Suppose that we choose γ so that its homologygroups are trivial, H ( γ ) ∼ = ker S ∼ = , H ( γ ) ∼ = coker S ∼ = . (149)Then, we have the isomorphisms, H (Γ) ∼ = H (Γ (cid:48) ) , H (Γ) ∼ = H (Γ (cid:48) ) , (150)equivalently, ker S ∼ = ker S (cid:48) coker S ∼ = coker S (cid:48) . (151)Thus, the spaces of cycles and conserved charges before and after the reduction are isomorphic when γ has trivialhomologies. Example 6 in Sec. IV corresponds to this situation, where the partial stoichiometric matrix is given by S = (cid:0) − (cid:1) , whose kernel and cokernel are trivial.The exact sequence applies as long as the commutativity condition, ker S ⊂ ker S , is satisfied, and we canconsider more general cases with ker S (cid:54) = . If the connecting map δ : H (Γ (cid:48) ) → H ( γ ) is a zero map, the longexact sequence (147) results in the following two exact sequences,0 (cid:47) (cid:47) H ( γ ) (cid:47) (cid:47) H (Γ) (cid:47) (cid:47) H (Γ (cid:48) ) (cid:47) (cid:47) , (152)0 (cid:47) (cid:47) H ( γ ) (cid:47) (cid:47) H (Γ) (cid:47) (cid:47) H (Γ (cid:48) ) (cid:47) (cid:47) . (153)This implies the isomorphisms, ker S/ ker S ∼ = ker S (cid:48) , coker S/ coker S ∼ = coker S (cid:48) . (154)Note that ker S consists of only locally supported global cycles, due to the assumption ker S ⊂ ker S . Theisomorphisms (154) represent equivalence of chemical reaction networks up to locally supported global cycles andlocally supported global conserved charges (emergent conserved charges are also absent when δ is a zero map, as wesee below). The connecting map is identified as follows. An element c ∈ H (Γ (cid:48) ), can be included in C (Γ (cid:48) ). ϕ is surjective and there exists c = (cid:18) c c (cid:19) such that ϕ ( c ) = c . From the commutativity of the diagram (140), we have ϕ ( S c ) = S (cid:48) c = . From the exactness ofthe row of Eq. (140), there exists d ∈ C ( γ ) such that ψ ( d ) = S c . We obtain [ d ] ∈ H ( γ ) by identifying the differences in im S .More explicitly, [ d ] = [ S c + S c ] = [ S c ]. The mapping c (cid:55)→ [ S c ] is the connecting map. The well-definedness of the map(indifference to the choice of c ) is obvious in this expression. δ is a zero map if S c ∈ im S = (coker S ) ⊥ , (155)for any c ∈ H (Γ (cid:48) ). Below we show that, if every conserved charge in γ is obtained by the projection of a globalconserved charge in Γ (namely, when there is no emergent conserved charge), the connecting map δ is a zero map.For a given d ∈ coker S , there exists an element of coker S , d T = ( d T , d T ). The condition d T S = reads d T S = , (156) d T S + d T S = , (157)where we used d T S = . Let us pick c ∈ H (Γ (cid:48) ) = ker S (cid:48) . The quantity d T S c can be shown to vanish as follows: d T S c = (157) − d T S c = S (cid:48) c = − d T S S +11 S c = (156) . (158)Therefore, we have shown d T S c = 0 for any d ∈ coker S and c ∈ ker S (cid:48) . This is equivalent to S c ∈ (coker S ) ⊥ .The relation between the long exact sequence and the numbers of cycles and conserved charges of various types issummarized in Fig. 3. The vertical lines represent the spaces, and the kernels are shown in black. Since it is an exactsequence, the kernel and image coincide at each space, such as im ψ = ker ϕ and so on. The exactness is the key tothe connections between cycles and conserved charges of particular types. Let us see an example. The image of δ is the space of emergent conserved charges, im δ = (cid:101) D ( γ ). They are emergent, because the image of δ is the kernelof ¯ ψ , and there is no counterpart in Γ. The connecting map δ provides us with a one-to-one mapping between anemergent cycle in Γ (cid:48) and an emergent conserved charge in γ (elements of ker δ are not emergent, since they can bewritten as an image of ϕ due to the exactness). The numbers d ( γ ), d (cid:48) ( γ ) in Fig. 3 are the same as the dimensionsof the spaces (129) and (130) that we defined previously.Compare Fig. 3 also with Fig. 6 in the Appendix B 2, where we discuss the relation between the numbers of cyclesand conserved charges and the structure of the A-matrix. The long exact sequence is valid when (cid:101) c ( γ ) = 0 (i.e., whenthe diagram (140) commutes). This implies (cid:101) d (cid:48) ( γ ) = 0 and there is not emergent conserved charge in Γ (cid:48) , since ¯ ϕ issurjective. H ( γ ) → H (Γ) → H (Γ′ ) → H ( γ ) → H (Γ) → H (Γ′ ) → 0 : emergent CCs of ˜ d γ : global cycles supported on c γ : cycles of that go through c ′ ΓΓ\ γ : emergent cycles of ˜ c ′ Γ′ : local CCs obtained by projections of global CCs d : CCs of that have no counterpart in d ′ Γ γ : conserved charges of d ′ Γ′ δ = i m δ k e r ¯ ψ ¯ ψ ¯ φ = i m ¯ ψ k e r ¯ φ = i m φ k e r δ φ ψ = i m ψ k e r φ CC = conserved charges
FIG. 3. Long exact sequence and conserved charges/cycles of various types.
C. Reduction of buffering structures
Here we present the main result, regarding the reduction of buffering structures. The following theorem representsa particularly nice property of buffering structures under reductions. We show that the steady-state concentrations5and rates of the network obtained by reducing a buffering structure are exactly the same as those of the networkbefore reduction, without any modification of parameters. Thus, the reduction of a buffering structure preserves thesteady-state properties of the boundary degrees of freedom. The theorem only relies on topological information of thenetwork and is true regardless of the kinetics.Theorem . Let Γ be a regular chemical reaction network with kinetics r ( x ) and let γ be an output-complete sub-network of Γ. We assume that the subnetwork γ does not have an emergent conserved charge. We consider a reducednetwork Γ (cid:48) = Γ /γ . If γ is a buffering structure, we have the isomorphisms,ker S/ ker S ∼ = ker S (cid:48) , coker S/ coker S ∼ = coker S (cid:48) , (159)Furthermore, when ( r , x ) is steady-state reaction rates concentrations of Γ, whose components we separate into thosein γ and Γ \ γ as r = (cid:18) r r (cid:19) , x = (cid:18) x x (cid:19) , (160)then, ( r , x ) is a steady-state solution of Γ (cid:48) . Remark . Let us comment on the assumption of the absence of emergent conserved charges. Under the assumptionof the regularity, the appearance of emergent conserved charges in an output-complete subnetwork γ is quite unlikely.In fact, in the case of monomolecular reaction networks, we can prove (cid:101) d ( γ ) = 0 for a connected and output-completesubnetwork γ , assuming that Γ is regular (see Appendix C 2), and this condition is redundant. So far, the examplesof buffering structures with nonzero emergent conserved charges are pathological in some sense. Presently, we havenot been able to prove the absence of emergent conserved charges for a generic (sound) reaction network, and thus itis assumed. We have more discussions on this point in Appendix C. Proof.
The regularity of Γ requires λ ( γ ) ≥ λ ( γ ) = (cid:101) c ( γ ) + d l ( γ ) = 0 . (161)Since (cid:101) c ( γ ) and d l ( γ ) are nonnegative integers, we have (cid:101) c ( γ ) = 0 and d l ( γ ) = 0. Since (cid:101) c ( γ ) = 0, we can use the longexact sequence (147). Because there is no emergent conserved charge, (cid:101) d ( γ ) = 0, by assumption, the connecting map δ is a zero map in the long exact sequence. This proves Eq. (159).Let us proceed to the latter part of the claim. The steady-state condition of Γ is written as S r ( x ) = 0 , (162) d ¯ α · x = (cid:96) ¯ α . (163)As usual, we divide the degrees of freedom to those in γ and Γ \ γ . Then Eq. (162) is written as (cid:18) S S S S (cid:19) (cid:18) r ( x , x ) r ( x ) (cid:19) = . (164)The reactions r ( x ) only depend on x , because γ is chosen to be output-complete. The first equation can be solvedfor r as r = − S +11 S r + c , with c ∈ ker S , and we have S (cid:48) r ( x ) = − S c = , (165)where the last equality is due to (cid:101) c ( γ ) = 0, that is equivalent to ker S ⊂ ker S .Let us turn to the conserved charges. Recall that d l ( γ ) is written as d l ( γ ) = | (coker S ) /X ( γ ) | . When d l ( γ ) = 0, wehave coker S = X ( γ ) = (cid:26)(cid:18) d d (cid:19) ∈ coker S (cid:12)(cid:12)(cid:12)(cid:12) d ∈ coker S (cid:27) . (166)Correspondingly, we can divide the basis vectors of coker S into two classes, { d ¯ α } = { d ¯ α γ , d ¯ α (cid:48) } , according to whether d is zero or not, d ¯ α γ = (cid:18) d ¯ α γ d ¯ α γ (cid:19) with d ¯ α γ (cid:54) = , d ¯ α (cid:48) = (cid:18) d ¯ α (cid:48) (cid:19) . (167)6With this basis, Eq. (163) is written as d ¯ α γ · x + d ¯ α γ · x = (cid:96) ¯ α γ , (168) d ¯ α (cid:48) · x = (cid:96) ¯ α (cid:48) . (169)In fact, d ¯ α (cid:48) is a conserved charge in Γ (cid:48) , d ¯ α (cid:48) ∈ coker S (cid:48) , as we see in the following. Since d ¯ α (cid:48) ∈ coker S , it satisfies (cid:0) ( d ¯ α (cid:48) ) T (cid:1) (cid:18) S S S S (cid:19) = (cid:0) ( d ¯ α (cid:48) ) T S ( d ¯ α (cid:48) ) T S (cid:1) = . (170)This implies that d ¯ α (cid:48) satisfies ( d ¯ α (cid:48) ) T S (cid:48) = ( d ¯ α (cid:48) ) T ( S − S S +11 S ) = , (171)hence d ¯ α (cid:48) ∈ coker S (cid:48) . Thus we have obtained a mapping,coker S (cid:51) d ¯ α (cid:48) = (cid:18) d ¯ α (cid:48) (cid:19) (cid:55)→ d ¯ α (cid:48) ∈ coker S (cid:48) . (172)This map is nothing but the induced map ¯ ϕ . It is important to note that this map is a surjection, that is evidentfrom the long exact sequence (147) . The equations satisfied by the boundary part (denoted by 2) of the concentra-tions/rates of Γ are Eqs. (165) and (169). Since all the conserved charges in Γ (cid:48) is given as a image ¯ ϕ , we find thatthe set of Eqs. (165) and (169) are exactly the same as the steady-state condition for the reduced network Γ (cid:48) , S (cid:48) r (cid:48) ( x (cid:48) ) = , (173) d ¯ α (cid:48) · x (cid:48) = (cid:96) ¯ α (cid:48) . (174)where x (cid:48) = x and r (cid:48) ( x (cid:48) ) = r ( x ). Thus, the steady-state solution of Γ (cid:48) should also be the steady-state solution ofΓ for the boundary degrees of freedom. This concludes the proof. D. Hierarchy of subnetworks
Let us consider nested subnetworks γ (cid:48) ⊂ γ ⊂ Γ. Given the stoichiometric matrix S of the whole network, we denotethe stoichiometric matrices of the subnetworks γ and γ (cid:48) by S γ and S γ (cid:48) , respectively. The submatrices are included inthe following form, S = (cid:18) S γ ∗∗ ∗ (cid:19) , S γ = (cid:18) S γ (cid:48) ∗∗ ∗ (cid:19) , (175)where ∗ indicates arbitrary matrices. Let us consider the situation where γ and γ (cid:48) has no emergent cycle and emergentconserved charge in Γ, namely (cid:101) c ( γ ) = (cid:101) d ( γ ) = 0, and (cid:101) c ( γ (cid:48) ) = (cid:101) d ( γ (cid:48) ) = 0. Under those assumptions, the quotient formulaof the generalized Schur complement [51] holds, S/S γ = ( S/S γ (cid:48) ) / ( S γ /S γ (cid:48) ) , (176)This indicates the isomorphisms of homology groups, H n (Γ /γ ) ∼ = H n ((Γ /γ (cid:48) ) / ( γ/γ (cid:48) )) , (177)for n = 0 ,
1. Thus, when we perform the reductions of nested subnetworks that have no emergent cycles and emergentconserved charges, the order of the reduction of them does not matter. Another way to see this is by Eq. (B13). We have (cid:101) c ( γ ) = (cid:101) d ( γ ) = 0 from the assumption, and (cid:101) c (cid:48) ( γ ) = (cid:101) d ( γ ) holds by the connecting map δ . Thus, we have (cid:101) d (cid:48) ( γ ) = 0, which means that there is no emergent charge in Γ (cid:48) . VI. EXAMPLE OF REDUCTION: METABOLIC PATHWAY OF
E. COLI
As an application of the reduction method, let us examine the central metabolism of
E. coli . We use the stoichiom-etry matrix presented in Ref. [42], which is constructed based on Ref. [1] with minor modifications. The networkstructure is shown in Fig. 4a, which consists of the glycolysis, the pentose phosphate pathway (PPP), and the tricar-boxylic acid cycle (TCAC). The list of the reactions for this system is given in Sec. VI A. Buffering structures in thisnetwork have been identified in Ref. [21] and there are in total 17 buffering structures, which we list in Sec. VI B. Aswe showed in Sec. III C, the intersections or unions of buffering structures are also buffering structures. They form ahierarchy, and such an architecture can be regarded as a source of robustness against perturbations, since bufferingstructures work as a kind of firewalls.Let us now perform reductions of buffering structures, under which the steady state is ensured to be the same asthe original network as we showed in Sec. V C. We denote the whole network by Γ. We can pick a buffering structure γ , which is a part of the pentose phosphate pathway (the yellow subnetwork in Fig. 4a) and given by γ = ( { X5P, S7P, E4P } , { , , , , } ) , (178)and perform reduction to obtain Γ := Γ /γ . The stoichiometric matrix of the reduced reaction network can becomputed by Eq. (103). The resulting network is shown in Fig. 4b. The reduction procedure induces rewiring ofreactions, which are colored in magenta in Fig. 4. The reaction 15 and 22 are rewired, and 22 is now an emptyreaction. The fraction 1 / necessary if we want the steady state to be the same as those of the original network. Otherwise,the steady state is changed in general. We can proceed further and reduce the subnetwork ( { G3P,R5P } , { , , } )(colored in red and orange in Fig. 4a). This reduction is the same as reducing γ ∪ γ from Γ. The result of thereduction is shown in Fig. 4c. Again, rewiring occurs and the reactions 5,6,15,36 are modified from the original system.Finally, let us focus on the part colored in green in Fig. 4a, which consists of the following subsets of chemical speciesand reactions, ( { G6P , F6P , F16P , , Ru5P , DHAP } , { , , , , , , , , , , } ) . (179)The complement of the subset (179) is given by γ ∪ γ ∪ γ , which hence is a buffering structure, and a reductioncan be performed. The structure of the reduced network Γ = Γ / ( γ ∪ γ ∪ γ ) is given in Fig. 4d. Compared to theoriginal network, we notice that the reactions 15 and 43 are rewired.To demonstrate our theoretical prediction, we numerically solve the rate equations for the four systems in Fig. 4 (theoriginal network Γ and the reduced ones, Γ , Γ , Γ ), using the same initial condition and reaction rate constants inall of the four cases (see Appendix D for details of parameter values). The time-series of concentrations are presentedin Fig. 5. After the initial transient dynamics, the original system approaches a (stable) steady state [Fig. 5(a)]. Wecan see that the reduced systems can reproduce the steady-state concentrations that the original system eventuallyreaches, although they have distinct short-time dynamics [Fig. 5(b-d)].In this way, buffering structures work as a guide as to how to perform the reduction and simplify a complex reactionnetwork. As long as the reduced part is a buffering structure and we use the generalized Schur complement (103) asa stoichiometric matrix of a reduced network, the steady-state concentrations and rates of the remaining part staythe same as the original ones regardless of the details of the kinetics, as a consequence of Theorem 3. A. List of reactions
1: Glucose + PEP → G6P + PYR.2: G6P → F6P.3: F6P → G6P.4: F6P → F16P.5: F16P → G3P + DHAP.6: DHAP → G3P.7: G3P → → PEP.9: PEP → In fact, γ = γ ∩ γ , and if we allow taking intersections of buffering structures, γ is redundant. This is consistent with Corollary 1. G6P 6PGRu5PF6PF16PG3PGlucose3PGPEPDHAPCO2 LactateAcCoA AcetateEthanolOAA CIT ICT2-KGSUCFUMMAL
Glyoxylate
X5PS7PE4P R5PPEP PYR
43 13 1416152,3456 3621,2219,20 17,18 4078,910,1142 41 33 31,3238 45 444624 25
28 35 23 3739341 12
CO2CO2PYR (a) Central metabolic pathway Γ of
E. coli . G6P 6PGRu5PF6PF16PG3PGlucose3PGPEPDHAP LactateAcCoA AcetateEthanolOAA CIT ICTSUCFUMMAL
Glyoxylate
R5PPEP PYR
43 13 16152,3456 3622 4078,910,1141 33 38 45 444624
28 35 23 3739341 12 CO22-KG CO2CO2
42 31,32
PYR (b) Reduced network Γ = Γ /γ , where γ is colored inyellow in (a). The fraction 1 / G6P 6PGRu5PF6PF16PGlucose3PGPEPDHAP LactateAcCoA AcetateEthanolOAA CIT ICTSUCFUMMAL
Glyoxylate
PEP PYR
43 13 16152,3456 368,910,1141 33 38 45 444624
35 23 3739 CO22-KG
CO2CO2
42 31,32
PYR (c) Γ = Γ / ( γ ∪ γ ), where γ ∪ γ is colored in yellow,red, and orange in (a). G6P 6PGRu5PF6PF16PDHAP
43 13 1416152,3456 36 (d) Γ = Γ / ( γ ∪ γ ∪ γ ), where γ ∪ γ ∪ γ is colored inyellow, red, orange, and blue in (a). FIG. 4. Central metabolic pathway (a) of
E. coli and reduced networks (b, c, d). In the reduced networks, rewired reactionsunder the reductions are colored in magenta.
10: PEP → PYR.11: PYR → PEP.12: PYR → AcCoA + CO2.13: G6P → → Ru5P + CO2.9
FIG. 5. Time-series of concentrations of DHAP, F16P, F6P, G6P, PG6, Ru5P computed by solving the whole system Γ (a) andthe reduced systems Γ , Γ , Γ (b, c, d). The same initial condition and reaction rate constants are used in the four cases.
15: Ru5P → X5P.16: Ru5P → R5P.17: X5P + R5P → G3P + S7P.18: G3P + S7P → X5P + R5P.19: G3P + S7P → F6P + E4P.20: F6P + E4P → G3P + S7P.21: X5P + E4P → F6P + G3P.22: F6P + G3P → X5P + E4P.23: AcCoA + OAA → CIT.24: CIT → ICT.25: ICT → → SUC + CO2.27: SUC → FUM.28: FUM → MAL.29: MAL → OAA.30: OAA → MAL.31: PEP + CO2 → OAA.32: OAA → PEP + CO2.33: MAL → PYR + CO2.34: ICT → SUC + Glyoxylate.35: Glyoxylate + AcCoA → MAL.36: 6PG → G3P + PYR.37: AcCoA → Acetate.38: PYR → Lactate.39: AcCoA → Ethanol.40: R5P → (output).41: OAA → (output).42: CO2 → (output).43: (input) → Glucose.44: Acetate → (output).45: Lactate → (output).46: Ethanol → (output). B. List of buffering structures γ = ( { Glucose } , { } ), γ = ( { Glucose , PEP , G6P , F6P , F16P , DHAP , G3P , , PYR , , Ru5P , X5P , R5P , S7P , E4P , AcCoA , OAA , CIT , ICT , , SUC , FUM , MAL , CO2 , Glyoxylate , Acetate , Lactate , Ethanol } , { , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , } ), γ = ( { F16P } , { } ), γ = ( { DHAP } , { } ),0 γ = ( { G3P , X5P , S7P , E4P } , { , , , , , , } ), γ = ( { } , { } ), γ = ( { Glucose , PEP , , PYR , AcCoA , OAA , CIT , ICT , , SUC , FUM , MAL , CO2 , Glyoxylate , Acetate , Lactate , Ethanol } , { , , , , , , , , , , , , , , , , , , , , , , , , , , } ), γ = ( { X5P , S7P , E4P } , { , , , , } ) , γ = ( { CIT } , { } ), γ = ( { } , { } ), γ = ( { SUC } , { } ) , γ = ( { FUM } , { } ) γ = ( { Glyoxylate } , { } ), γ = ( { X5P , R5P , S7P , E4P } , { , , , , , } ) , γ = ( { Acetate } , { } ), γ = ( { Lactate } , { } ), γ = ( { Ethanol } , { } ). VII. SUMMARY AND OUTLOOK
The main theme of the present paper has been the relation between the structure and functions of the chemical reac-tion network. As a characterization of the structure, homology and cohomology groups for chemical reaction networksare introduced, in which the actions of boundary and coboundary operators are determined by the stoichiometry. Theelements of homology groups correspond to cycles and conserved charges of chemical reaction networks. Steady statesare shown to be determined by the elements of the cohomology groups. In a similar way to the homology and coho-mology groups of topological spaces, the Mayer-Vietoris sequence and the long exact sequence of a pair of chemicalreaction networks are introduced. The latter one is particularly useful in studying the reduction of reaction networks.We have proposed a method of reduction of chemical reaction networks. The reduced network is characterized by thestoichiometric matrix obtained by eliminating the chemical species and reactions of an output-complete subnetworkvia the Schur complementation. The reduction relies only on the stoichiometry, which determines the topology ofthe reaction networks, and hence is applicable to any kind of kinetics . This is advantageous because it is difficult todetermine the kinetics and parameters of the reactions experimentally in many biological systems. For tracking thechange of cycles and conserved charges under the reductions, the tools of algebraic topology, such as the long exactsequence, have been useful. We have studied how the law of localization can be understood from this perspective. Wehave shown that the influence index is expressed in terms of the numbers of cycles/conserved charges of particulartypes, as in Eq. (138). We have also shown that the influence index is a submodular function over output-completesubnetworks. A corollary of this is that buffering structures are closed under intersection and union, which is usefulwhen we enumerate the buffering structures of a large reaction network. As a central result of the paper, we haveshown that, buffering structures, which are subnetworks with vanishing influence index, behave nicely under thereduction. Namely, under the reduction of a buffering structure, the steady state of the rest of the network stays thesame as the original network. This is the content of Theorem 3. The theorem justifies the intuition that bufferingstructures are regarded as an ‘irrelevant’ substructures: they can be eliminated safely through the reduction methodproposed here without changing the long-time behavior of the system. The reduction procedure introduces rewiringof reactions, and it is necessary so that the steady state is not modified under the reduction. As an application ofthe reduction method, we discussed the reduction of the central metabolic pathway of
E. coli and illustrated thatreactions are rewired non-trivially under the reduction. We also demonstrated the invariance of the steady state underthe reduction of buffering structures by numerically solving the rate equations before and after the reduction.Our results highlight that special care should be taken when we simplify a reaction network. A naive eliminationof a subnetwork not of interest would alter steady-state properties of the original system. As long as the subnetworkhas the vanishing influence index and reactions are rewired appropriately using the generalized Schur complement, itcan be eliminated while keeping the steady state intact.We believe that the mathematical formulation that we discussed to characterize the topology of chemical reactionnetworks will be useful for understanding the static and dynamical properties of reaction systems. The eigenvalues ofthe Laplacian operators entail the information of the topology of the network connectivity. Steady states correspond tothe eigenvectors with zero eigenvalues and they incorporate crudest topological information of the reaction network.The eigenvectors with higher eigenvalues are going to be needed if we want to extend the reduction method toapproximate the dynamics as well as the steady states. In Ref. [52], morphisms of chemical reaction networks are considered and a condition is given as to when a reaction network candynamically emulate another one. ACKNOWLEDGMENTS
This work is in part supported by RIKEN iTHEMS Program. Y. Hirono is supported by the Korean Ministry ofEducation, Science and Technology, Gyeongsangbuk-do and Pohang City at the Asia Pacific Center for TheoreticalPhysics (APCTP) and by the National Research Foundation (NRF) funded by the Ministry of Science of Korea(Grant No. 2020R1F1A1076267). H. M. is supported by JSPS KAKENHI Grant (19K23413). Y. Hidaka is supportedby JSPS KAKENHI Grant Numbers 17H06462. Y. Hirono is grateful to Benjamin for helpful discussions and theconstant encouragement.
Appendix A: Laplace operators and Hodge decomposition
In this section, we discuss the Hodge decomposition and Laplace operators, which are closely related to the coho-mology groups introduced in the main text.We can define Laplace operators, ∆ n : C n (Γ) → C n (Γ), as∆ := d † d , ∆ := d d † . (A1)Recall that the coboundary operator (19) and its adjoint (23) are given by (d a )( e A ) = (cid:80) i ( S T ) Ai a ( v i ) for a ∈ C (Γ)and (d † a )( v i ) = (cid:80) A S iA a ( e A ) for a ∈ C (Γ). The action of the Laplacians are written in the matrix form as(∆ a )( v i ) = (cid:88) j ( SS T ) ij a ( v j ) , (∆ a )( e A ) = (cid:88) B ( S T S ) AB a ( e B ) , (A2)for a ∈ C (Γ) and a ∈ C (Γ). Those are generalizations of the graph Laplacian to hypergraphs. The properties ofhypergraph Laplacians were discussed recently in Refs. [53–55]. When all the reaction is monomolecular, the Laplacianreduces to the graph Laplacian of the directed graph.The space C (Γ) admits the following orthogonal decomposition, C (Γ) = im d ⊕ ker ∆ . (A3)This is a natural generalization of the Hodge decomposition of flows on networks [56] to the case of a hypergraph.Thus, given a 1-cochain f ∈ C (Γ), we can decompose it in a unique way as f = d a + c, (A4)where c ∈ ker d † ∩ ker d is a harmonic cochain and a ∈ C (Γ). This is the Hodge decomposition associated with thecomplex (27). By acting d † on Eq. (A4), we haved † f = d † d a = ∆ a. (A5)We can solve this for the potential a as a = ∆ +0 d † f + a . (A6)Here, ∆ +0 : C (Γ) → C (Γ) is the operator defined by (∆ +0 b )( v i ) := (cid:80) j ( SS T ) + ij b ( v j ) for b ∈ C (Γ), where M + indicates the Moore-Penrose inverse of a matrix M , and a ∈ ker ∆ . The harmonic component c can be obtained by c = f − d a = (1 − d ∆ +0 d † ) f. (A7)Using the properties of the Moore-Penrose inverse, the action of the operator that appears on the RHS of Eq. (A7)is written as [(1 − d ∆ +0 d † ) b ]( e A ) = (cid:88) B (1 − S + S ) AB b ( e B ) , (A8)for an arbitrary b ∈ C (Γ). The matrix 1 − S + S is the projection matrix to ker S . Thus, the harmonic componentcan be identified by the projection to ker S , c ( e A ) = (cid:88) B (1 − S + S ) AB f ( e B ) . (A9)2This is consistent with the fact that c ∈ H (Γ) = ker S . The potential a can be obtained by the multiplication of theMoore-Penrose inverse of S to f , a ( v i ) = (cid:88) A ( S T ) + iA f ( e A ) + a ( v i ) , (A10)where a ∈ ker ∆ . Appendix B: Cycles and conserved charges1. Interpretation of (cid:101) c ( γ ) and d l ( γ ) , and (cid:101) d ( γ ) Let us here discuss intuitive interpretations of the integers appearing in the decomposition (138). We will refer tothe elements of ker S as “global cycles” and those of coker S as “global conserved charges.” Similarly, the elementsker S and coker S are referred to as “local cycles” and “local conserved charges.” With this terminology, theelements of (ker S ) supp γ are called as “locally supported global cycles.”We first look at (cid:101) c ( γ ). We denote the space by (cid:101) C( γ ) := ker S / (ker S ) supp γ . (B1)Then, (cid:101) c ( γ ) = | (cid:101) C( γ ) | . To clarify its meaning, we represent the space (B1) as follows, (cid:101) C( γ ) = { c ∈ C (Γ) | S c = ¯ P γ v , P γ c = c } / { c ∈ C (Γ) | S c = 0 , P γ c = c , } = { c ∈ C (Γ) | S c = v (cid:54) = , v ∈ C (Γ \ γ ) , P γ c = c } . (B2)An element of (cid:101) C( γ ) is a local cycle that is not a global cycle, by which we mean that c ∈ (cid:101) C( γ ) has its boundary inΓ \ γ . For example, let us take the subnetwork γ = ( { v } , { e , e } ) of a monomolecular reaction network, e e e v v v v γ (B3)Although e + e ∈ C (Γ) has its support in γ , its boundary, ∂ ( e + e ) = − v + v , (B4)is outside of γ . The element e + e is a local cycle, since ∂ ( e + e ) is zero as a relative chain in C ( γ ) = C (Γ) /C (Γ \ γ ).Note that the network (B3) as a whole does not have a cycle and ker S = . Thus, we can identify c ∈ (cid:101) C( γ ) to bea local cycle whose boundary is out of γ . When c is viewed in Γ, it may be extended to a global cycle, but it doesnot have to be. Considering its meaning, we will refer to the elements of (cid:101) C( γ ) as emergent cycles , which only appearwhen we focus on a subnetwork.Let us illustrate the space (cid:101) C( γ ) pictorially. The matrix S works as a boundary operator on the space of chemicalreactions. Thus, the kernel of S are linear combinations of reactions without boundaries. Cycles and non-cycles canbe drawn pictorially as Cycles
Non-cycles , Recall that the boundary of each reaction is specified by the stoichiometric matrix as ∂ e A = (cid:80) i ( S T ) Ai v i . γ . The space ker S is spanned by local cycles in γ , for example, ker S = γ ,where the inner box represents a subnetwork γ and the red lines constitute the basis of the space. Here, the symbolmeans that the cut ends are reactions and not chemical species. See the following two choices for example: e e e v v v v γ e e e v v v v γ For the left one, cut ends are both reactions. For the right one, the cut ends are a species and a reaction. The bothends have to be reactions so that the cut cycle can be a local cycle. An element of ker S may be extended to aglobal cycle, or it may be a part of a global noncycle.On the other hand, the space (ker S ) supp γ for the same configuration takes into account only the global cyclessupported on γ , (ker S ) supp γ = γ .Therefore, the coset space is generated by the following elements, ker S / (ker S ) supp γ = γ .As we see in the figure, the space ker S / (ker S ) supp γ consists of local cycles that are not global cycles.We can similarly interpret conserved charges. The transpose of the stoichiometric matrix, S T , can be regarded asa boundary operator acting on C (Γ), which is the space of chemical species. In this sense, an element of coker S hasno boundary, with respect to this boundary operator. We here visualize this in a similar way to the cycles, Conserved charges
Not conserved .Note that the boundary of the box is identified. The filled circles represent a source or a drain of chemical species,because of which the charge is not conserved.The space coker S represents the global conserved charges, and P γ (coker S ) is the projection of coker S to γ ,4 P γ (coker S ) = γ Here, how the conserved charges are cut does not matter. On the other hand, the space coker S is generated by thered and green elements in the following figure, coker S = γ Emergent conserved charge ,where filled and open rectangles mean boundaries with chemical species and reactions, respectively. The parts wedenoted by green lines, , are emergent conserved charges , which are conserved when it is seen in a subnetworkbut not conserved in Γ. In fact, the appearance of emergent conserved charges typically leads to “unphysical” systems,in the sense that either a steady state does not exist or the matrix A is not invertible and the response of the systemto the perturbation of parameters is not well-defined. For example, one can consider the following network andsubnetwork, e e e v v v γ The whole network does not have a conserved charge, but the subnetwork γ has one, v + v + v . However, sucha reaction network cannot reach a steady state, since the concentration of v continues to increase. When we takethe difference | P γ (coker S ) | − | coker S | , we can count the number of lost conserved charges minus the number ofemergent conserved charges (if any), | P γ (coker S ) | − | coker S | = γ − γ ,where the part colored in red in the first term indicates lost conserved charges, that are conserved in Γ but theirprojections to γ are not. This equation is equal to the latter two terms of the decomposition (138), d l ( γ ) − (cid:101) d ( γ ). Example . Consider a monomolecular network Γ = (
V, E ) = ( { v , v , v } , { e , e } ) with the following structure, e e v v v γ We discuss more on this point in Appendix C. γ = ( { v } , { e , e } ) that is indicated by a box. The whole network does not have a cycle, andthe subnetwork γ has one emergent cycle given by c = e + e . Also, Γ has one conserved charge, d = v + v + v .Its projection to γ is given by v and it is not a conserved charge in γ . So we have one lost conserved charge. Eachinteger appearing in the decomposition of λ ( γ ) is (cid:101) c ( γ ) = 1 , d l ( γ ) = 1 , (cid:101) d ( γ ) = 0 , (B5)and λ ( γ ) = 2. For the same Γ, let us consider a different choice of a subnetwork, e e v v v γ The subnetwork γ does not have a cycle, and there is one lost conserved charge, so we have (cid:101) c ( γ ) = 0 , d l ( γ ) = 1 , (cid:101) d ( γ ) = 0 , λ ( γ ) = 1 . (B6) Example . Consider a network (
V, E ) = ( { v , v , v } , { e , e , e } ) with the following structure, e e e v v v γ If we choose a subnetwork γ = ( { v } , { e , e } ), (cid:101) c ( γ ) = 1 , d l ( γ ) = 1 , (cid:101) d ( γ ) = 0 , λ ( γ ) = 2 . (B7)The subnetwork γ has one emergent cycle and one lost conserved charge and the influence index is 2.
2. Embedding of A-matrices
It is useful to look at the A-matrix to visualize the relations among cycles/conserved charges of various types insubnetworks and reduced networks.Let us first summarize the notations. In this section, we suppress the dependence on γ for notational simplicity.General rules are as follows. Quantities with a tilde are emergent ones, and we use character c for cycles and d forconserved charges. The numbers with a prime are associated with Γ (cid:48) . The relevant numbers are listed as follows: • v, v (cid:48) : number of chemical species in γ, Γ (cid:48) • e, e (cid:48) : number of chemical reactions in γ, Γ (cid:48) • (cid:101) c, (cid:101) c (cid:48) : number of emergent cycles of γ, Γ (cid:48) • (cid:101) d, (cid:101) d (cid:48) : number of emergent conserved charges of γ, Γ (cid:48) • c, c (cid:48) : number of cycles of Γ, whose projections to γ, Γ (cid:48) are also cycles of γ, Γ (cid:48) • d, d (cid:48) : number of conserved charges of Γ, whose projections to γ, Γ (cid:48) are also conserved in γ, Γ (cid:48) • ¯ c (cid:48) , ¯ d (cid:48) : number of cycles/conserved charges of Γ that are locally supported in Γ (cid:48) • c, d : number of cycles/conserved charges of Γ that have nonzero support in γ d e c v d ′ cycles of supported in Γ γ v ′ e ′ ˜ d ′ c ′ = v ′ e ′ ˜ c ′ c ′ d ′ emergent cycles of Γ′ c ˜ c d ˜ d ve emergent cycle in γ emergent conserved charge in γ → → ¯ d ′ ¯ d ′ ¯ c ′ ¯ c ′ local conserved charge in obtained by the projection of global conserved charges γ conserved charges of that have no counterpart in Γ γ conserved charges of obtained from those of Γ′ Γ conserved charges of supported inside ΓΓ′ cycles of supported in
Γ Γ′ A A γ A Γ′ d conserved charges of that have nonzero support in Γ γ measures non-squareness of this matrix λ cycles of obtained from the cycles of Γ′ Γ emergent conserved charge in Γ′ c cycles of that have nonzero support in Γ γ FIG. 6. Embedding of the A-matrices for a generic output-complete subnetwork.
In Fig. 6, we illustrate a more detailed structure of the matrix A than Fig. 2. In the center is the matrix A ofthe total system Γ. We choose an output-complete subnetwork γ , and bring the rows/columns related to γ to theupper-left part. Then the matrix A looks like one in the center of Fig. 6. We consider an output-complete subnetwork γ , and the A-matrix of γ , which we denote by A γ , is shown in the upper-left part of Fig. 6. The part surrounded by apink rectangle is the common part of A γ and A . The subnetwork γ can in general contain additional (i.e., emergent)cycles and conserved charges, whose numbers are denoted by (cid:101) c and (cid:101) d . Because the matrix A γ is square, we have therelation, e + d + (cid:101) d = v + c + (cid:101) c. (B8)This equation is in fact the same as Eq. (51). Similarly, we can consider the matrix A for the network Γ (cid:48) = Γ /γ obtained by reducing γ from Γ. The numbers of the emergent cycles and emergent conserved charges in Γ (cid:48) aredenoted by (cid:101) c (cid:48) and (cid:101) d (cid:48) . The matrix A Γ (cid:48) is also square and we have e (cid:48) + d (cid:48) + (cid:101) d (cid:48) = v (cid:48) + c (cid:48) + (cid:101) c (cid:48) . (B9)The influence index is given by λ := e + d + d (cid:48) − ¯ d (cid:48) − v − c, (B10) When (cid:101) c ( γ ) >
0, the equation of motion contains some terms that cannot be determined, as in Eq. (98). Here, we formally considera reduced network Γ (cid:48) which is defined with the generalized Schur complement S (cid:48) . In this sense, the property of the reduced networkdefined this way cannot be fully constrained from the properties of Γ and γ . d = d + d (cid:48) − ¯ d (cid:48) . Using Eq. (B8), we can alsoexpress λ as λ = (cid:101) c + d (cid:48) − ¯ d (cid:48) − (cid:101) d, (B11)which is the same the decomposition (138). We can also consider a similar quantity that measure the non-squarenessof the lower-right part, λ (cid:48) := v (cid:48) + c (cid:48) − e (cid:48) − ¯ d (cid:48) = (cid:101) d (cid:48) + d (cid:48) − ¯ d (cid:48) − (cid:101) c (cid:48) , (B12)where the second expression is obtained using Eq. (B9). In fact, due to the squareness of the whole matrix A , λ (cid:48) isequal to the influence index, λ = λ (cid:48) . This results in the following relation, (cid:101) c − (cid:101) d + (cid:101) c (cid:48) − (cid:101) d (cid:48) = 0 . (B13) Appendix C: Emergent conserved charges in chemical reaction networks
In this section, we discuss the role of emergent conserved charges in chemical reaction networks.
1. Systems with emergent conserved charges
As far as we observe, the chemical reaction systems with emergent conserved charges in output-complete subnetworksare pathological , in either of the following senses:(A) The steady-state condition, (cid:88) A S iA r A ( x, k ) = 0 , (cid:88) i d ¯ αi x i = (cid:96) ¯ α , (C1)does not fully determine the steady-state solution, and arbitrary parameters have to be introduced to specifythe solution.(B) No steady-state solution exists.(C) The reaction kinetics is unphysical.Below, we discuss some examples of each case. a. Pattern A : solutions have arbitrary parameters An example of pattern (A) is given by ddt (cid:18) xy (cid:19) = (cid:18) − − (cid:19) (cid:18) r r (cid:19) , (cid:18) r r (cid:19) = (cid:18) k yk y (cid:19) . (C2)The matrix A for this system is A = (cid:18) k k (cid:19) . (C3)This is not invertible. The steady state solution for the mass-action kinetics is given by (cid:18) ¯ r ¯ r (cid:19) = , (cid:18) ¯ x ¯ y (cid:19) = (cid:18) m (cid:19) , (C4)where m is an arbitrary parameter. If we choose a subnetwork γ = { x } , there is an emergent conserved charge. Letus consider the fluctuations around the steady state, ddt (cid:18) δxδy (cid:19) = (cid:18) − − (cid:19) (cid:18) k δyk δy (cid:19) = (cid:18) k δy − ( k + k ) δy (cid:19) , (C5)8where δx ( t ) := x ( t ) − ¯ x indicates the fluctuation from the steady state. The fluctuation associated with the emergentcharge, δx , is a zero mode. This means that the system is not asymptotically stable.Generically, when we have to introduce arbitrary parameters m , the matrix A has a null vector, as we see below.The steady-state condition reads r A ( x ( k , (cid:96) , m ) , k ) = − (cid:88) α µ α ( k , (cid:96) , m ) c αA , (C6) (cid:88) i d ¯ αi x i ( k , (cid:96) , m ) = (cid:96) ¯ α . (C7)By taking the derivative of those equations with respect to m , we find (cid:18) r A,i c αA d ¯ αi (cid:19) ∂∂m a (cid:18) x i µ α (cid:19) = . (C8)This means that A has a null vector and det A = 0. b. Pattern B: no steady-state solution An example of pattern (B) is given by the following reaction system with the mass-action kinetics, ddt x x x = − − − r ( x , x ) r ( x ) r r , r r r r = k x x k x k k . (C9)The steady-state solution does not exist in general. We need to fine-tune the parameters to have a solution. When k = k is satisfied, we have a steady state ¯ r = k (cid:0) (cid:1) T . (C10)The matrix A is A = ∂ r ∂ r ∂ r
10 0 0 10 0 0 1 , (C11)where ∂ j r i := ∂r i /∂x j and it is evaluated at the steady state. This matrix is not regular, det A = 0.Let us choose an output-complete subnetwork γ = ( { v , v } , { e } ). The matrix A for the subnetwork is A γ = (cid:18) ∂ r ∂ r − (cid:19) . (C12)The subnetwork has an emergent conserved charge, (cid:101) d T = (cid:0) − (cid:1) . The time derivative of this charge is ddt (cid:101) d T x γ = ddt (cid:0) − (cid:1) (cid:18) x x (cid:19) = r − r = k − k . (C13)Although (cid:101) d T S (cid:54) = , where (cid:101) d T := ( (cid:101) d T a ), for any parameter a , when the steady state exists, k = k , the combination (cid:101) d T x is in fact a conserved charge of the whole system. It is not conserved unless the parameters are fine-tuned.Let us consider the fluctuations around the steady state, ddt δx δx δx = − − − δr ( x , x ) δr ( x ) δr δr . (C14)The fluctuation associated with the emergent conserved charge leads to a zero mode, ddt δ ( x − x ) = δ ( r − r ) = 0 . (C15)9 c. Pattern C: example with emergent conserved charges and unphysical kinetics Here we discuss an example that has a subnetwork with vanishing influence index and also has an emergent conservedcharge, while the kinetics is unphysical. The rate equation of this system is ddt x x x x = − − − − − r ( x , x ) r ( x , x ) r ( x , x ) r ( x ) . (C16)We have added catalytic dependencies in the reactions r ( x , x ) and r ( x , x ). The stoichiometric matrix has atrivial kernel and ¯ r A = 0 at the steady state. The cokernel of S is also trivial.The matrix A is A = ∂ r ∂ r ∂ r ∂ r ∂ r ∂ r ∂ r . (C17)Its determinant is in general nonvanishing, det A = ∂ r ∂ r ∂ r ∂ r . (C18)Let us consider an output-complete subnetwork γ = ( { x , x } , { r , r } ). The index of γ is zero, λ ( γ ) = − − , (C19)and hence it is a buffering structure. The matrix A of the local system reads A γ = ∂ r ∂ r ∂ r − . (C20)Although γ is a buffering structure, the subnetwork γ has one emergent cycle and one emergent conserved charge.For the mass-action kinetics, r ( x , x ) r ( x , x ) r ( x , x ) r ( x ) = k x x k x x k x x k x , (C21)the steady-state concentrations are ¯ x ¯ x ¯ x ¯ x = m m or m (cid:48) m (cid:48) , (C22)where m , m , m (cid:48) , m (cid:48) are arbitrary parameters. With this kinetics, det A = 0.Let us instead employ the following kinetics, r ( x , x ) r ( x , x ) r ( x , x ) r ( x ) = k ( x + x ) k ( x + x ) k ( x + x ) k x , (C23)all the concentration vanish at the steady state, ¯ x i = 0. The matrix A is now invertible, det A = k k k k (cid:54) = 0.Although A is regular, the sensitivity is trivial, ∂ A ¯ x i = 0, since ∂ A ¯ r B = 0 at the steady state. The kinetics (C23) isnot physically sound, because the reaction r ( x , x ) can be nonzero even if the concentration of the reactant x iszero (note that x is catalytic). The same is true for r .0Let us consider the fluctuations around the steady state of the emergent conserved charge. ddt ( δx − δx ) = − δr ( x ) = − r , δx , (C24)where we denote r i,j := ∂ j r i . The time derivative of δx is ddt δx = − r , δx + r , δx ) − ( r , δx + r , δx ) − r , δx = − r , δx − r , δx + ( − r , − r , − r , ) δx . (C25)Hence, there is a zero mode when r , = 0 at the steady state, and then A is not invertible. d. Emergent conserved charges and zero modes of fluctuations We denote the matrix R whose components are given by r i,j = ∂ j r i and and separate it into block matrices, R = (cid:18) R R R R (cid:19) , (C26)according to the separation x = (cid:18) x x (cid:19) . The linear fluctuations around the steady state satisfy the following equationsof motion, ddt (cid:18) δ x δ x (cid:19) = SRδ x = (cid:18) S R δ x + ( S R + S R ) δ x S R δ x + ( S R + S R ) δ x (cid:19) , (C27)where we have also separated the stoichiometric matrix into submatrices, and we used R = , which follows fromthe output-completeness. We consider the fluctuation associated with the emergent conserved charge (cid:101) d , ddt (cid:101) d T δ x = (cid:101) d T S R δ x . (C28)Since it is an emergent charge, (cid:101) d T S (cid:54) = . The time derivative of the RHS of Eq. (C28) reads ddt (cid:101) d T S R δ x = (cid:101) d T S R S R δ x + ( · · · ) δ x . (C29)Therefore, an emergent conserved charge results in a zero mode when (cid:101) d T S R S R vanishes. We are not awareof a physical example in which (cid:101) d T S R S R does not vanish. In the example given by Eq. (C16), the first termof Eq. (C29) is computed as (cid:101) d T S R S R δ x = (cid:0) − (cid:1) (cid:18) (cid:19) (cid:18) r , r , r , (cid:19) (cid:18) − (cid:19) (cid:18) r , r , r , (cid:19) (cid:18) δx δx (cid:19) = (cid:0) r , r , (cid:1) (cid:18) δx δx (cid:19) . (C30)When this vanishes, the matrix A acquires a zero mode and is not invertible.
2. Absence of emergent conserved charges in monomolecular reaction networks
Here we consider monomolecular reaction networks and we show that, if there exists a nonzero emergent charge inan output-complete subnetwork, the index λ ( γ ) is necessarily negative. If the index is negative in an output-completesubnetwork, the matrix A is not invertible, and the response of the system to the parameter-perturbation is notwell-defined. We here show the following statement: Theorem . Suppose that γ is a connected and output-complete subnetwork of a monomolecular reaction network Γ.If (cid:101) d ( γ ) >
0, then the influence index λ ( γ ) is negative.1 Proof.
In order to have an emergent conserved charge in γ , all the boundaries of γ should be chemical species andnot reactions, in a monomolecular reaction network. Then, all the reaction in γ should end on the chemical speciesinside γ , which means S = .Recall that an emergent cycle is c ∈ ker S which is not a cycle of the whole network, S (cid:18) c (cid:19) = (cid:18) S c (cid:19) (cid:54) = . (C31)When S = , there is no such c meaning that all the local cycles are also a global cycle. Namely, for a given c ∈ ker S , (cid:18) c (cid:19) ∈ ker S always holds. Thus, we have (cid:101) c ( γ ) = 0.This also results in d l ( γ ) = d (cid:48) ( γ ) − ¯ d (cid:48) ( γ ) = 0 as follows. If d (cid:48) ( γ ) − ¯ d (cid:48) ( γ ) is nonzero, there should exist d and d such that [recall the definitions of spaces, Eqs. (130) and (128)] d T S + d T S = , (C32) d T S + d T S = , (C33) d T S (cid:54) = . (C34)However, when S = , Eqs. (C32) and (C34) are contradictory. Thus d (cid:48) ( γ ) − ¯ d (cid:48) ( γ ) = 0 and we have d l ( γ ) = 0.Therefore, we have shown that (cid:101) c ( γ ) = 0 and d l ( γ ) = 0, and the index is written as λ ( γ ) = (cid:101) c ( γ ) + d l ( γ ) − (cid:101) d ( γ ) = − (cid:101) d ( γ ) , (C35)which is negative due to the assumption (cid:101) d ( γ ) > Appendix D: Parameter values used in Figure 5
In Fig. 5, for an illustration purpose, we employ the mass-action kinetics, where the rate of the i -th reaction isgiven by the product of its substrate concentrations, r i = k i (cid:81) A ( x A ( t )) y iA (see Eq. (3) for the definition of y iA ).In the simulation, the initial concentrations and the reaction rate constants are chosen randomly: x =0 . x AcCoA = 0 . x Acetate = 0 . x CIT = 0 . x CO2 = 0 . x DHAP = 0 . x E4P = 0 . x Ethanol = 0 . x F16P = 0 . x F6P =0 . x FUM = 0 . x G3P = 0 . x G6P = 0 . x Glucose = 0 . x Glyoxylate = 0 . x ICT = 0 . x KG2 = 0 . x Lactate =1 . , x MAL = 0 . x OAA = 1 . , x PEP = 0 . x PG3 = 1 . , x PYR = 0 . x R5P = 0 . x Ru5P = 0 . x S7P = 0 . x SUC =0 . x X5P = 0 . k = 1, k = 4 . k = 7 . k = 5 . k = 3 . k = 9 . k = 5 . k = 6 . k = 3 . k = 9 . k =2 . k = 6 . k = 4 . k = 3 . k = 7 . k = 2 . k = 3 . k = 5 . k = 5 . k = 4 . k = 8 . k =7 . k = 9 . k = 1 . k = 9 . k = 7 . k = 7 . k = 8 . k = 6 . k = 6 . k = 6 . k = 7 . k =9 . k = 6 . k = 1 . k = 9 . k = 4 . k = 5 . k = 7 . k = 3 . k = 8 . k = 9 . k = 4 . k =2 . k = 8 . k = 3 . [1] N. Ishii, K. Nakahigashi, T. Baba, M. Robert, T. Soga, A. Kanai, T. Hirasawa, M. Naba, K. Hirai, A. Hoque, et al. ,“Multiple high-throughput analyses monitor the response of e. coli to perturbations,” Science no. 5824, (2007)593–597.[2] J. Munger, B. D. Bennett, A. Parikh, X.-J. Feng, J. McArdle, H. A. Rabitz, T. Shenk, and J. D. Rabinowitz,“Systems-level metabolic flux profiling identifies fatty acid synthesis as a target for antiviral therapy,”
Naturebiotechnology no. 10, (2008) 1179–1186.[3] N. Zamboni, S.-M. Fendt, M. R¨uhl, and U. Sauer, “13 c-based metabolic flux analysis,” Nature protocols no. 6, (2009)878.[4] N. Barkai and S. Leibler, “Robustness in simple biochemical networks,” Nature no. 6636, (1997) 913–917.[5] U. Alon, M. G. Surette, N. Barkai, and S. Leibler, “Robustness in bacterial chemotaxis,”
Nature no. 6715, (1999)168–171.[6] H. Kitano, “Biological robustness,”
Nature Reviews Genetics no. 11, (2004) 826–837.[7] H. Kitano, “Towards a theory of biological robustness,” Molecular systems biology no. 1, (2007) 137.[8] G. Shinar and M. Feinberg, “Structural sources of robustness in biochemical reaction networks,” Science no. 5971,(2010) 1389–1391.[9] A. Larhlimi, S. Blachon, J. Selbig, and Z. Nikoloski, “Robustness of metabolic networks: A review of existing definitions,”
Biosystems no. 1, (2011) 1–8. . [10] J. Whitacre, “Biological robustness: Paradigms, mechanisms, and systems principles,” Frontiers in Genetics (2012) 67. .[11] J. M. Eloundou-Mbebi, A. K¨uken, N. Omranian, S. Kleessen, J. Neigenfind, G. Basler, and Z. Nikoloski, “A networkproperty necessary for concentration robustness,” Nature communications no. 1, (2016) 1–7.[12] M. S. Okino and M. L. Mavrovouniotis, “Simplification of mathematical models of chemical reaction systems,” Chemicalreviews no. 2, (1998) 391–408.[13] T. J. Snowden, P. H. van der Graaf, and M. J. Tindall, “Methods of model reduction for large-scale biological systems: asurvey of current methods and trends,” Bulletin of mathematical biology no. 7, (2017) 1449–1486.[14] J. Kuo and J. Wei, “Lumping analysis in monomolecular reaction systems. analysis of approximately lumpable system,” Industrial & Engineering chemistry fundamentals no. 1, (1969) 124–133.[15] J. Wei and J. C. Kuo, “Lumping analysis in monomolecular reaction systems. analysis of the exactly lumpable system,” Industrial & Engineering chemistry fundamentals no. 1, (1969) 114–123.[16] Z. Zi, “Sensitivity analysis approaches applied to systems biology models,” IET Systems Biology (November, 2011)336–346(10). https://digital-library.theiet.org/content/journals/10.1049/iet-syb.2011.0015 .[17] D. Degenring, C. Froemel, G. Dikta, and R. Takors, “Sensitivity analysis for the reduction of complex metabolismmodels,” Journal of Process Control no. 7, (2004) 729 – 745. . Dynamics, Monitoring, Control andOptimization of Biological Systems.[18] M. Apri, M. de Gee, and J. Molenaar, “Complexity reduction preserving dynamical behavior of biochemical networks,” Journal of Theoretical Biology (2012) 16 – 26. .[19] S. Danø, M. F. Madsen, H. Schmidt, and G. Cedersund, “Reduction of a biochemical model with preservation of its basicdynamic properties,”
The FEBS journal no. 21, (2006) 4862–4877.[20] S. R. Taylor, F. J. Doyle 3rd, and L. R. Petzold, “Oscillator model reduction preserving the phase response: applicationto the circadian clock,”
Biophysical journal no. 4, (2008) 1658–1673.[21] T. Okada and A. Mochizuki, “Law of localization in chemical reaction networks,” Phys. Rev. Lett. (Jul, 2016)048101. https://link.aps.org/doi/10.1103/PhysRevLett.117.048101 .[22] T. Okada and A. Mochizuki, “Sensitivity and network topology in chemical reaction systems,”
Phys. Rev. E (Aug,2017) 022322. https://link.aps.org/doi/10.1103/PhysRevE.96.022322 .[23] T. Okada, J.-C. Tsai, and A. Mochizuki, “Structural bifurcation analysis in chemical reaction networks,” Phys. Rev. E (Jul, 2018) 012417. https://link.aps.org/doi/10.1103/PhysRevE.98.012417 .[24] B. Fong, “Decorated cospans,” Theory and Applications of Categories no. 33, (2015) 1096–1120, arXiv:1502.00872 .[25] J. C. Baez and B. S. Pollard, “A compositional framework for reaction networks,” Reviews in Mathematical Physics no. 09, (Sep, 2017) 1750028. http://dx.doi.org/10.1142/S0129055X17500283 .[26] J. C. Baez and K. Courser, “Structured cospans,” 2020.[27] R. Rao and M. Esposito, “Nonequilibrium thermodynamics of chemical reaction networks: Wisdom from stochasticthermodynamics,” Phys. Rev. X (Dec, 2016) 041064. https://link.aps.org/doi/10.1103/PhysRevX.6.041064 .[28] M. Polettini and M. Esposito, “Irreversible thermodynamics of open chemical networks. i. emergent cycles and brokenconservation laws,” The Journal of chemical physics no. 2, (2014) 07B610 1.[29] S. Klamt, U.-U. Haus, and F. Theis, “Hypergraphs and cellular networks,”
PLoS Comput Biol no. 5, (2009) e1000385.[30] M. Heiner, D. Gilbert, and R. Donaldson, “Petri nets for systems and synthetic biology,” in International school onformal methods for the design of computer, communication and software systems , pp. 215–264, Springer. 2008.[31] V. N. Reddy, M. L. Mavrovouniotis, M. N. Liebman, et al. , “Petri net representations in metabolic pathways.” in
ISMB ,vol. 93, pp. 328–336. 1993.[32] H. Kacser and J. Burns, “The control of flux.” in
Symposia of the Society for Experimental Biology , vol. 27, p. 65. 1973.[33] D. A. Fell, “Metabolic control analysis: a survey of its theoretical and experimental development,”
Biochemical Journal no. 2, (1992) 313–330.[34] E. Szathmary, “Do deleterious mutations act synergistically? metabolic control theory provides a partial answer.”
Genetics no. 1, (1993) 127–132.[35] R. MacLean, “Predicting epistasis: an experimental test of metabolic control theory with bacterial transcription andtranslation,”
Journal of evolutionary biology no. 3, (2010) 488–493.[36] H. C. Bagheri and G. P. Wagner, “Evolution of dominance in metabolic pathways,” Genetics no. 3, (2004)1713–1735.[37] M. Feinberg, “The existence and uniqueness of steady states for a class of chemical reaction networks,”
Archive forRational Mechanics and Analysis no. 4, (1995) 311–370.[38] M. Feinberg, “Multiple steady states for chemical reaction networks of deficiency one,”
Archive for Rational Mechanicsand Analysis no. 4, (1995) 371–406.[39] G. Craciun and M. Feinberg, “Multiple equilibria in complex chemical reaction networks: I. the injectivity property,”
SIAM Journal on Applied Mathematics no. 5, (2005) 1526–1546.[40] G. Craciun and M. Feinberg, “Multiple equilibria in complex chemical reaction networks: Ii. the species-reaction graph,” SIAM Journal on Applied Mathematics no. 4, (2006) 1321–1338.[41] A. Mochizuki and B. Fiedler, “Sensitivity of chemical reaction networks: A structural approach. 1. examples and thecarbon metabolic network,” Journal of Theoretical Biology (2015) 189 – 202. . [42] B. Fiedler and A. Mochizuki, “Sensitivity of chemical reaction networks: a structural approach. 2. regular monomolecularsystems,” Mathematical Methods in the Applied Sciences no. 16, (2015) 3519–3537. https://onlinelibrary.wiley.com/doi/abs/10.1002/mma.3436 .[43] B. Brehm and B. Fiedler, “Sensitivity of chemical reaction networks: A structural approach. 3. regular multimolecularsystems,” Mathematical Methods in the Applied Sciences no. 4, (2018) 1344–1376. https://onlinelibrary.wiley.com/doi/abs/10.1002/mma.4668 .[44] M. Feinberg, Foundations of chemical reaction network theory . Springer, 2019.[45] M. Koecher, “The generalized inverse of integral matrices,”
Linear Algebra and its Applications (1985) 187 – 198. .[46] A. van der Schaft, S. Rao, and B. Jayawardhana, “On the mathematical structure of balanced chemical reaction networksgoverned by mass action kinetics,” SIAM Journal on Applied Mathematics no. 2, (2013) 953–973.[47] S. Rao, A. van der Schaft, and B. Jayawardhana, “A graph-theoretical approach for the analysis and model reduction ofcomplex-balanced chemical reaction networks,” Journal of Mathematical Chemistry no. 9, (2013) 2401–2422.[48] G. Kron, “Tensor analysis of networks,” New York (1939) .[49] F. D¨orfler, J. W. Simpson-Porco, and F. Bullo, “Electrical networks and algebraic graph theory: Models, properties, andapplications,”
Proceedings of the IEEE no. 5, (2018) 977–1005.[50] F. Dorfler and F. Bullo, “Kron reduction of graphs with applications to electrical networks,”
IEEE Transactions onCircuits and Systems I: Regular Papers no. 1, (2013) 150–163.[51] D. Carlson, E. Haynsworth, and T. Markham, “A generalization of the schur complement by means of the moore–penroseinverse,” SIAM Journal on Applied Mathematics no. 1, (1974) 169–175.[52] L. Cardelli, “Morphisms of reaction networks that couple structure to function,” BMC systems biology no. 1, (2014) 84.[53] J. Jost and R. Mulas, “Hypergraph laplace operators for chemical reaction networks,” Advances in Mathematics (2019) 870 – 896. .[54] R. Mulas and D. Zhang, “Spectral theory of laplace operators on chemical hypergraphs,” arXiv preprintarXiv:2004.14671 (2020) .[55] J. Jost and R. Mulas, “Normalized laplace operators for hypergraphs with real coefficients,” arXiv preprintarXiv:2010.10100 (2020) .[56] X. Jiang, L.-H. Lim, Y. Yao, and Y. Ye, “Statistical ranking and combinatorial hodge theory,”
MathematicalProgramming127