Structural reduction of chemical reaction networks based on topology

We develop a model-independent reduction method of chemical reaction systems based on the stoichiometry, which determines their network topology. A subnetwork can be eliminated systematically to give a reduced system with fewer degrees of freedom. This subnetwork removal is accompanied by rewiring of the network, which is prescribed by the Schur complement of the stoichiometric matrix. Using homology and cohomology groups to characterize the topology of chemical reaction networks, we can track the changes of the network topology induced by the reduction through the changes in those groups. We prove that, when certain topological conditions are met, the steady-state chemical concentrations and reaction rates of the reduced system are ensured to be the same as those of the original system. This result holds regardless of the modeling of the reactions, namely chemical kinetics, since the conditions only involve topological information. This is advantageous because the details of reaction kinetics and parameter values are difficult to identify in many practical situations. The method allows us to reduce a reaction network while preserving its original steady-state properties, thereby complex reaction systems can be studied efficiently. We demonstrate the reduction method in hypothetical networks and the central carbon metabolism of Escherichia coli.


I. INTRODUCTION
Chemical reactions in living systems form complex networks [1][2][3]. They operate in a highly coordinated manner, and are responsible for various cellular functions. Experimentally, high-throughput measurements have been conducted to study cellular responses to perturbations for the purpose of elucidating underlying regulatory mechanisms (see, for example, Refs. [4][5][6]). 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 provide detailed quantitative predictions, a faithful modeling is challenging for most biological systems, because our prior knowledge about kinetics and parameter values is limited, and also because many parameters are difficult to measure experimentally. Furthermore, the complexity of models may confound model-independent features with model-dependent ones.
To address these difficulties, it is desirable to reduce complex reaction systems to simpler ones. A reduction is practically useful since it can reduce the number of variables and parameters needed to be included in the analysis, and it can also identify features 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 [7][8][9][10][11][12][13][14]. Chemical reaction networks inside living organisms are highly interconnected, and yet are robust under internal fluctuations and environmental perturbations. If a system is insensitive to the details of its substructure, it is natural to expect that a reduction is possible, in the spirit of renormalization. To the best of our knowledge, the reduction methods [15,16] of chemical reaction systems studied so far are based on timescale separation, lumping [17,18], sensitivity analysis [19][20][21] or optimizations [22,23]. To apply those approaches, we need detailed information about the reactions. For example, to exploit the timescale separation, we should know which reactions are fast and which are slow. The sensitivity analysis also requires the dependence of the system on various parameters.
In this paper, we develop a systematic method of reducing chemical reaction networks based on their topology (see Figs. 1 and 2). One motivation for the reduction method comes from the law of localization [24][25][26]; if a certain topological index, which we call the influence index, is zero for a subnetwork, perturbations inside the subnetwork do not affect the steady state of the remaining elements of the network (see Sec. III for the precise statement). This observation indicates that certain subnetworks are 'irrelevant', as far as the remaining part of the network is concerned. As we will show, a reduction can be systematically performed through the Schur complementation of the stoichiometric matrix with respect to a subset of chemical species and reactions. The well-definedness of the reduction process requires that the subnetwork should satisfy a condition called the 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 concentrations and reaction rates of the reduced system are exactly the same as those of the original system, as far as the remaining degrees of freedom are concerned. For a given subnetwork γ satisfying a condition called output-completeness, we assign a nonnegative integer that we call the influence index, λ(γ). A subnetwork with vanishing influence index is called a buffering structure. Although an elimination of a subnetwork (γ2 in the figure) generally modifies the steady state of the remaining part of the network, a buffering structure (γ1) can be reduced while preserving the original steady state of the remaining part. Example of the reduction. In a chemical reaction network with a stoichiometric matrix S, if a subnetwork γ with Sγ is output-complete, the system can be reduced to a smaller system, whose stoichiometric matrix S ′ is given by the generalized Schur complement S ′ := S/Sγ . The reduced system can reproduce the same steady-state properties of the original system, if the influence index of the subnetwork is zero. Note that S ′ is in general different from the corresponding submatrix of S (compare the lower-right block of S with S ′ , where the difference between them is indicated by the colored components in S ′ ). This alteration is pictorially represented as 'rewiring' of the network (e.g., the head of e5 is rewired to v3). See Fig. 5 for an application to the central carbon metabolism of E. coli.
We emphasize that those conditions are topological ones and determined solely by the network structure; hence, are insensitive to the details of how the reactions are modeled. Thus, the result is broadly applicable, because it holds regardless of the kinetics or parameter values. This is of practical merit since the kinetics of reactions or the values of parameters are difficult to identify in many situations. To characterize the topology of reaction networks, we introduce the homology and cohomology groups for chemical reaction networks. The change of the topology of chemical reaction networks under the reduction is captured by the change of the (co)homology groups. The tools of algebraic topology are convenient for tracking those changes. We recommend the readers who are interested in practical aspects of reduction to directly go to Sec. IV, where we discuss the reduction procedure with simple examples.
The rest of the paper is organized as follows. In Sec. II, we introduce concepts for characterizing the structure of chemical reaction networks. Further, 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. In Sec. III, we review the structural sensitivity analysis and the law of localization. We also show that the influence index is submodular as a function over output-complete subnetworks. In Sec. IV, we introduce the reduction procedure and illustrate the method with simple examples. In Sec. V, we discuss the relation between the structural sensitivity analysis and the reduction method. We show that the reduction of a buffering structure, that is an output-complete subnetwork with vanishing influence index, has a particularly nice property: The reduced system admits the same steady states as the original system. In Sec. VI, as an application to realistic networks, we demonstrate the reduction method for the metabolic pathways of Escherichia coli. Section VII is devoted to summary and outlook. In Appendix A, we discuss the Hodge decomposition and Laplace operators for chemical reaction networks. In Appendix B, we provide intuitive interpretations of the cycles and conserved charges of various types that appear in the decomposition of the influence index. We also illustrate how the decomposition of the index can be seen visually in the structure of the A-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 provide the details of the metabolic pathways of E. coli. discussed in Sec. VI.

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 1 (Chemical reaction network). A chemical reaction network (CRN) Γ is 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, 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 created in reaction e A . An element of N V will be referred to as a chemical complex. The system can be an open reaction network 1 , when there is a reaction whose source or target function is zero for 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. A reaction is usually represented in the following form, where v i ∈ V , and y iA andȳ iA are nonnegative integers. Those integers are given by the source and target functions as The stoichiometry of the reaction is specified by the stoichiometric matrix S, whose components are given by Remark 1. There are several equivalent ways to formulate a chemical reaction network such as a hypergraph [32] or a Petri net [33,34]. Remark 2. A reaction that involves at most one chemical species as reactants and products, such as v 1 → v 2 , is called monomolecular. When all the reactions in the system are monomolecular, the corresponding reaction network is a usual directed graph. In this case, the stoichiometric matrix is the incidence matrix of the graph. If we regard Γ as a directed 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 whose bases are chemical species/reactions. We denote the resulting vector spaces as Elements of those spaces are referred to as 0-chains and 1-chains. Higher (n ≥ 2) chains do not exist in the current setting. The stoichiometric matrix provides us with natural boundary operators on the spaces of chains, The action of ∂ 1 is defined by its action on the basis e A ∈ C 1 (Γ) and v i ∈ C 0 (Γ), We often use the notation of linear algebra, where an element i a i v i ∈ C 0 (Γ) is represented by the vector a = (a 1 , a 2 , · · · ) T , and we also write a ∈ C 0 (Γ). For b ∈ C 1 (Γ), the action of the boundary operator is given by the multiplication of the stoichiometric matrix, On the spaces of chains, let us define inner products by With these inner products, we can define the adjoint of the boundary operator, The action on the basis v i ∈ C 0 (Γ) is given by In the linear-algebra notation, the action of ∂ † 1 is the multiplication of the transpose of S to a ∈ C 0 (Γ), Example 1. Let us consider a reaction network Γ = ({v 1 , v 2 , v 3 , v 4 , v 5 }, {e 1 , e 2 , e 3 , e 4 , e 5 , e 6 }) given by the following set of chemical reactions, The stoichiometric matrix of the network is It can be drawn as We represent a monomolecular reaction by a single arrow, and we use a rectangle to represent a multimolecular reaction. In this network, e 3 is a multimolecular reaction and others are all monomolecular. The action of the boundary operator is, for example, and so on. Those are intuitively understood from the figure.
The network is open, since we have inputs from the outside (e 1 and e 2 ) and outputs to the external world (e 5 and e 6 ). For example, s(e 1 )(v i ) = 0 for any v i ∈ V . The action of ∂ † 1 is for example. Namely, the operator ∂ † 1 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, for n = 0, 1. Given an x ∈ C 0 (Γ), x(v i ) ∈ R represents the concentration of the chemical species v i . Similarly, for a given r ∈ C 1 (Γ), 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 0 (Γ) and r ∈ C 1 (Γ), 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, where we have used the linearity of the map x. Thus, we can identify the coboundary operator that acts on the chemical concentration x ∈ C 0 (Γ) as the multiplication of the matrix S T . We define the inner product of n-cochains as 2 x, y 0 : where x, y ∈ C 0 (Γ) and r, s ∈ C 1 (Γ). With these inner products, the adjoint of the coboundary operator d n , d † n : C n+1 (Γ) → C n (Γ), is defined by f n+1 , d n g n n+1 = d † n f n+1 , g n n .
Following the definition, we can identify d † 0 as follows, where r ∈ C 1 (Γ) and x ∈ C 0 (Γ). Thus, the action of d † 0 is given by By construction, the adjoint of coboundary operator satisfies (d † 0 r)(v i ) = r(∂ † 1 v i ).
2 More generally, one may define the inner product with a weight function as f, g n := where w is a R-valued function over Cn(Γ).

B. Homology, cohomology, and steady states
With the (co)chains and (co)boundary operators defined above, we can discuss (co)homology groups. We have the following chain complex, Noting that the action of ∂ 1 is the multiplication of the stoichiometric matrix S, we can identify the homology groups as Remark 3. Note that C 0 (Γ) is endowed with a standard inner product, with respect to which we can take the orthogonal linear subspace (im S) ⊥ . Moreover, the restriction of the quotient map C 0 (Γ) → 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 0 (Γ). Note also that the orthogonal subspace (im S) ⊥ is the same as the kernel of the transpose of S, ker S T . Combined with the above observation, this implies that we can always identify coker S with ker S T ⊂ C 0 (Γ).
Similarly, with the coboundary operator d 0 , we can define a complex of cochains as The associated cohomology groups are where (−) ⊥ denotes taking orthogonal spaces with respect to the standard inner product on C 0 (Γ) and C 1 (Γ). An Euler number for this complex can be defined as where |W | indicates the dimension of the vector space W . Several remarks on the homology and cohomology groups are in order: Remark 4. Since we consider the R coefficients, the homology and cohomology groups are the same, H n (Γ) ∼ = H n (Γ) for n = 0, 1.
Remark 5. In the chemistry literature, the elements of H 1 (Γ) are referred to as cycles, and this is consistent with the mathematical terminology.
Remark 6. When the network is monomolecular and the corresponding network is a directed graph, the dimension |H 0 (Γ)| is the number of connected components.
Remark 7. Similarly to the homology groups of topological spaces, Laplace operators can be defined and we can perform Hodge decomposition of C 1 (Γ). See Appendix A.
The cohomology groups defined above are closely related to the steady states of a reaction network as we see below. Let us consider the time evolution of spatially homogeneous chemical concentrations. The change of the chemical concentration is driven by the reactions. The time derivative of the concentration of species v i is given by the divergence of the reaction rate, which is more explicitly written as To solve the rate equations, we have to specify kinetics of chemical reactions, such as the mass-action kinetics and the Michaelis-Menten kinetics. A reaction's kinetics gives the reaction rate r A as a function of its substrate concentrations (i.e., the concentrations of species with y iA > 0) and parameters, r A = r A (x; k A ), where k A represents any one of the parameters for the A-th reaction; for example, in the Michaelis-Menten kinetics, k A represents the Michaelis constant or the maximum rate. The elements of H 0 (Γ) and H 1 (Γ) 3 characterize the steady states of chemical reaction networks. The rate equation (32) at the steady state reads which means that the steady-state reaction rate is an element of the kernel of S, r ∈ ker S ∼ = H 1 (Γ). The cokernel of S is related to conserved quantities of the system. Given d ∈ coker S ∼ = H 0 (Γ), that satisfies i d i S iA = 0, we have Thus, the linear combination i d i x i is independent of time and hence is conserved. For this reason, we refer to the elements of coker S as conserved charges 4 . To find the steady-state solutions, we have to specify the value of all the conserved charges. A steady state is specified by an element of H 0 (Γ) and H 1 (Γ), where {dᾱ} and {c α } are basis vectors of H 0 (Γ) and H 1 (Γ), respectively. The coefficients µ α (k, ℓ) depend on the parameters k and ℓ.
Example 2. We consider a network Γ = ({v 1 , v 2 , v 3 , v 4 }, {e 1 , e 2 , e 3 , e 4 , e 5 }) with the following reactions, The network structure can be drawn as We here take the mass-action kinetics, and the equations of motion are written as 3 Although the natural choice is to consider x and r as the elements of cohomology groups, we can equivalently consider them as elements of homology groups, since they are isomorphic in the current setting. 4 "Conserved moiety" may be more chemistry-oriented terminology.
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 by where span {v 1 , v 2 , · · · } indicates the vector space spanned by vectors v 1 , v 2 , · · · . The cokernel is one-dimensional and the system has one conserved charge. To find the steady states, we need to specify the value of the charge as The steady-state reaction rates and concentrations arē x = k1 k2 k1 k3 where we set K := k 4 k 2 1 /k 2 k 3 . The vectorr is spanned by the basis vectors of ker S and their coefficients are µ α .

C. Subnetworks
Let us consider a subset of chemicals and reactions, γ ⊂ Γ, which 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 where the indices are restricted to those of the subnetwork, v i ∈ V γ , e A ∈ E γ . We denote the space of relative chains by where Γ \ γ := (V \ V γ , E \ E γ ) is the complement of the subnetwork γ. The homology and cohomology groups for the subnetwork can be defined similarly. The chain complex for a subnetwork γ is where the action of the boundary operator ∂ 1 on the basis of C 1 (γ) is defined with the partial stoichiometric matrix S γ , The associated homologies with the complex (46) are The Euler number for a subnetwork is given by Note that The value of the concentrations and reaction fluxes inside a subset γ are given by R-valued functions over the space of chemicals and reactions, 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 γ 1 , γ 2 ⊂ Γ, which consist of γ 1 = (V γ1 , E γ1 ) and γ 2 = (V γ2 , E γ2 ). We can consider the intersection and union of the subnetworks, The exact sequence (56) below explains the relationship among cohomology groups of γ 1 ∪ γ 2 , γ 1 , γ 2 and γ 1 ∩ γ 2 .
Regarding the family {γ 1 , γ 2 } as a 'covering' of γ 1 ∪ γ 2 , we can think of Eq. (56) as an analogue of the Mayer-Vietoris sequence associated with an open covering of a topological space. Following the usual technique in topology, we will derive the long exact sequence from a short exact sequence of chain complexes. We have the following short exact sequence of chain complexes, where the horizontal maps are given by By applying the snake lemma to Eq. (54), we obtain In general, if there is an exact sequence of finite-dimensional vector spaces, the alternating sum of the dimensions of them is equal to zero. Therefore, the exact sequence (56) implies

III. LAW OF LOCALIZATION
A sensitivity analysis studies the response of the system to the perturbations of reaction parameters or initial conditions (conserved charges). In the context of metabolic networks, a theoretical framework called the metabolic control analysis has been developed [35][36][37][38][39]. Under the mass-action framework, biologically insightful results have been obtained regarding the sensitivity to conserved charges [11,14] as well as stability properties of stable states [40][41][42][43], although the mass-action law is not necessarily appropriate for some biological systems. Among the studies on sensitivity analysis, the structural sensitivity analysis [44][45][46] aims at constraining the response of reaction systems from the network structure alone.
In this section, we first review the structural sensitivity analysis and the law of localization [24][25][26]. For a given subnetwork, we assign a nonnegative integer, which we call the influence index. The influence index is determined from the topology of the subnetwork, and plays a decisive role in structural sensitivity. When the influence index is zero, the perturbation of the parameters and conserved charges inside the subnetwork does not affect the rest of the network. Such a structure is called a buffering structure. In Sec. III C, we prove that the influence index is submodular as a function over subnetworks. As a corollary of this property, we show that buffering structures are closed under intersection and union.

A. Structural sensitivity analysis
At the steady state, the reaction rates and the chemical concentrations satisfy where {dᾱ} is a basis of coker S and the second equation specifies the values of conserved charges. Considerable effort has been devoted to the study of the existence or uniqueness of steady states under the mass-action kinetics [47]. In the current analysis, we assume the existence of a steady state, and we focus on how it is perturbed under the change of parameters. The steady-state values of the concentrations and reaction rates are determined by the values of rate parameters and conserved charges, {k A , ℓᾱ}. The reaction rates r A (x(k, ℓ), k A ) have explicit dependence on k A , and also dependence on k and ℓ through x i (k, ℓ). 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 We are interested in the sensitivity of the reaction rates and concentrations under the perturbation of the parameters, By taking the derivative of Eqs. (59) and (60) with respect to k B and ℓβ, we obtain the following equations, Note that r A (x(k, ℓ), k A ) depends explicitly on k A and also depends implicitly on k and ℓ through x. The equations can be compactly written in the matrix form, where ∂ B := ∂/∂k B , ∂β := ∂/∂ℓβ, and we have introduced a partitioned square matrix, where the upper-left block is an |A| × |i| matrix whose (A, i)-th element is given by ∂rA ∂xi 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| matrix consisting of the basis {dᾱ} of coker S, and the lower-right one the |ᾱ| × |α| zero matrix. Here, we use the notation that index i for chemicals runs from 1 to |i|. The matrix A is square due to the identity, One can see from Eq. (66) that the response to the change of the parameter is determined by the inverse of the matrix A, We refer to A as the A-matrix ("A" indicates that it is an augmented matrix). Its inverse, −A −1 determines the sensitivity of the system and is called the sensitivity matrix. If we partition A −1 as and noting that ∂ B r A is a diagonal matrix, ∂ B r A ∝ δ BA , the responses of steady-state concentrations and reaction rates [or equivalently, the coefficients µ α in Eq. (60)] to the perturbations of k A and ℓᾱ are given by In this paper, we consider the following class of chemical reaction systems: Definition 2 (Regularity of a chemical reaction network with kinetics). A chemical reaction network with kinetics is called 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, we assume the regularity unless otherwise stated so that A is invertible and the response of the system is well-defined. The regularity implies the asymptotic stability of the steady state, through the relation between det A and the determinant of the Jacobian [26].

B. Law of localization
Definition 3 (Output-completeness). When a subnetwork γ = (V γ , E γ ) satisfies the condition that E γ includes all the chemical reactions affected by V γ , γ is called output-complete. Definition 4 (Influence index). For an output-complete subnetwork γ, the influence index is defined by The definitions of the spaces that appear in the influence index are given as follows: where S is the stoichiometric matrix, P 0 γ and P 1 γ are the projection matrices to γ in the space of chemical species and reactions, respectively. Namely, (ker S) supp γ is the space of vectors of ker S supported inside γ, and P 0 γ (coker S) is the projection of coker S to γ. Here, recall from Remark 3 that we regard coker S as a subspace of C 0 (Γ) via the identification coker S ∼ = (im S) ⊥ . We will use similar identifications throughout this paper. Remark 8. The influence index is nonnegative, λ(γ) ≥ 0, for a regular chemical reaction network. It will be shown in the proof of Theorem 1. Theorem 1 (Law of localization). 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 under the perturbation of rate parameters or conserved charges inside γ. Definition 5 (Buffering structures). For a given chemical reaction network Γ, an output-complete subnetwork γ with the vanishing influence index, λ(γ) = 0, is called a buffering structure. Example 3. The influence index of the empty subnetwork is zero. The index of the whole network Γ is also zero, This is natural in the sense that there is no "outside" of the whole network. Example 4. Let us take the same network as Example 2. The A-matrix of this system is where r A,i := ∂r A /∂x i and it is evaluated at the steady state. With this matrix A, the responses of the concentration and reaction rates to the change of parameters and the value of conserved charges can be obtained by Eq. (69). The is also a buffering structure, λ(γ 2 ) = −2 + 2 − 1 + 1 = 0, which contains a cycle supported on γ 2 . This explains the fact that x 1 and x 2 do not depend on the value of conserved charge is also a buffering structure, λ(γ 3 ) = −3 + 3 − 1 + 1 = 0, and hence x 2 does not depend on k 2 , k 4 , k 5 , and ℓ.
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. 3. All the component of the lower-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 on the upper-left corner is from a square matrix. The numbers c and d in Fig. 3 are given by c = |(ker S) supp γ | and d = |P 0 γ (coker S)|, which appear in Eq. (72). Because of the assumption of regularity, we have det A = 0, and the black 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.
When λ(γ) = 0, the black box in the upper-left corner is a square matrix. Then, A −1 inherits the same structure, Namely, if we denote the generic index of (A −1 ) as µ, ν, · · · and write the index inside and outside γ as µ ⋆ and µ ′ , respectively, we have (A −1 ) µ ′ ν ⋆ = 0. Because of this structure, which means that the concentrations out of γ do not depend on the parameter k B inside γ. Consequently, we have where we used the fact that r A ′ only depends on the concentrations outside γ because of the output-completeness. The same is true for the perturbation of the conserved charge, C. Submodularity of the influence index The influence index λ(γ) can be regarded as a function over subnetworks. We here show that the influence index satisfies 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.
• A function f (γ) over a set is called submodular, when it satisfies When ≤ is replaced with ≥, the function satisfying the replaced equation is called supermodular.
Proof. We show that We note that χ(γ) is a modular function, meaning that it satisfies which is derived from the Mayer-Vietoris exact sequence (56). Thus, it suffices to show that the last two terms on the right-hand side (RHS) of Eq. (83) are submodular. In fact, we show that each of them is submodular. Let us first look at |P 0 γ (coker S)|. If denote W := coker S, the submodularity of |P 0 γ (coker S)| reads We prove this equation just after this proof. Thus, we have shown the submodularity of |P 0 γ (coker S)|. Next, we show that |(ker S) supp γ | is supermodular. Consider the following vector space, Its dimension is given by Since any element of Z is supported in Thus, we have shown that |(ker S) supp γ | is a supermodular function, and −|(ker S) supp γ | is submodular. Therefore, |P 0 γ (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 1 , · · · , 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 , we denote 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 0 γ W | = rank w[v γ , N ] for a subnetwork γ = (v γ , e γ ). Let us pick two subnetworks γ 1 and γ 2 , and denote the sets of chemical species by v 1 and v 2 , respectively. We consider the submatrix w[v 1 ∩ v 2 , N ]. We can pick a row basis as We can form a row basis of w[v 1 , N ] by adding row vectors from w[v 1 \ v 2 , N ] to α 1∩2 . Let α 1 be the picked indices, then We can further pick row vectors from w[v 2 \ v 1 , N ] and form a basis of w[v 1 ∪ v 2 , N ]. Let us denote the added indices as α 2 , then Since the vectors specified by the indices α 1∩2 ∪α 2 are linearly independent and α 1∩2 ∪α . This can be written as This is equivalent to Eq. (85).
Corollary 1. Let Γ be a regular chemical reaction network. The union and the intersection of two buffering structures inside Γ are also buffering structures.

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 the original system. Let us here introduce a reduction method which consists of the following two steps, (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, The reduced network is characterized by a new stoichiometric matrix, Crucial points are, how to identify a subnetwork to be eliminated, and how to obtain the new stoichiometric matrix S ′ , which determines the structure of the reduced network. In this section, we mainly discuss step (2). We will discuss more on the choice of a subnetwork in Sec. (V).

A. Reduction procedure
Let us here illustrate a method of reduction based on the network topology. We denote the whole reaction network by Γ = (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 chemical species and reactions inside γ as internal, and those in Γ \ γ as boundary. For the given subnetwork γ, we separate the chemical concentrations and reaction rates as where 1 and 2 correspond to internal and boundary degrees of freedom, respectively. Accordingly, the stoichiometric matrix S can be partitioned as Note that the submatrix S 11 is the same matrix as S γ that appeared in Sec. II C. Hereafter we use S 11 for notational convenience. With the separation of internal and boundary degrees of freedom, the rate equations of the whole reaction system is written as d dt While the internal reaction rates r 1 = r 1 (x 1 , x 2 ) in general depend on both of the internal and boundary chemical concentrations, when γ is chosen to be output-complete, the boundary reaction rates r 2 = r 2 (x 2 ) do not depend on the internal chemical concentrations x 1 . The first equation of Eq. (96) can be solved for r 1 as where S + 11 is the Moore-Penrose inverse of S 11 , and c 11 ∈ ker S 11 . Substituting this to the second equation of Eq. (96), we get When the following condition is satisfied, S 21 c 11 = 0 and the second term of the RHS of Eq. (98) vanishes. 5 Then, the rate equation is written as where S ′ is the generalized Schur complement, As long as steady states are concerned, the subnetwork (x 2 , r 2 ) satisfies the rate equation whose stoichiometric matrix is S ′ . This motivates us to consider the subnetwork (x 2 , r 2 ) whose rate equation is given by Based on the considerations above, we define the reduction of a reaction system in the following way: Definition 6 (Reduction). 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 11 . We define a reduced network Γ ′ = (V \ V γ , E \ E γ ) obtained by eliminating γ from Γ, by a stoichiometric matrix S ′ given by the generalized Schur complement (101). We denote the resultant reaction network by Γ ′ = Γ/γ. Accordingly, the chemical concentrations and reaction rates of the reduced system (x ′ , r ′ ) are obtained from the original ones (x, r) as and the rate equation of the reduced system is given by Remark 9. The structure of the reduced network is determined by the generalized Schur complement (101) of the stoichiometric matrix. The second term in Eq. (101) represents the rewiring of the network associated with the elimination of γ.
Remark 10. The reduced system can be always defined if γ is output-complete; otherwise, the reduction is ill-defined since the reduced system would depend on x 1 through r 2 . We emphasize that the output-completeness is a topological condition 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 depends on further nature of γ. In the following sections, we will discuss more on the features of the subnetworks that behave nicely under reductions. In Sec. V C, we prove that, when γ has a vanishing influence index (see Sec. III), which is determined by the network topology, the steady state of the reduced system is assured to be the same as the steady state of the original system. Remark 11. In Sec. IV C, we show that the reduction we introduced here can be regarded as a morphism of reaction networks that 'shrinks' a subnetwork to a point, followed by the removal of degenerate (chemically meaningless) reactions. Remark 12. We note that the elements of the matrix S ′ are rational, since the Moore-Penrose inverse of an integral matrix is rational [48]. The matrix S ′ can be always transformed into an integral matrix by columnwise rescaling of S ′ together with the rescaling of reaction rates. Remark 13. The stoichiometric matrix given by the generalized Schur complement has appeared previously in flux balance analysis [49,50]. The current method is different from the ones discussed for reaction networks with the mass-action kinetics in detailed balanced [51] and complex balanced [52] situations, where the Schur complementation is performed for the weighted Laplacian similarly to the Kron reduction of electrical circuits [53][54][55]. In the current formulation, the Schur complementation is performed for the stoichiometric matrix.

B. Simple examples of reduction
To illustrate the reduction procedure, here we discuss simple examples. In Sec. VI, we discuss the reduction of the metabolic pathway of E. coli as a more realistic example. Example 5. We consider a monomolecular reaction network that consists of (V, E) = ({v 1 , v 2 }, {e 1 , e 2 , e 2 }). We take a subnetwork γ = ({v 1 }, {e 2 }) to be reduced. Under the reduction, the stoichiometric matrix changes as where we have brought the reduced part to the upper-left part. The reduction looks like The original rate equation is d dt where x 1 = x(v 1 ), r 2 = r(e 2 ) and so on. If we eliminate r 2 (x 1 ), .
The reduced equation of motion is obtained by replacing x 2 + x 1 with x 2 on the left-hand side.
To compute the steady-state solutions, let us for example employ the mass-action kinetics,  The steady-state reaction rates and concentrations are given by  The steady-state solutions of the reduced system are Note that this solution of the reduced system is exactly the same as the solution (111) of the original system for the boundary concentrations and rates. Indeed, this is a special property of buffering structures. In this example, the subnetwork has a vanishing influence index, λ(γ) = −1 + 1 − 0 + 0 = 0, and hence is a buffering structure. Generically, when the reduced subnetwork is a buffering structure, the steady-state solution of the reduced system is the same as the original system, and this is the content of Theorem 4. Although we used the mass-action kinetics in this example, the theorem applies to any kind of kinetics. We give a proof of the theorem in Sec. V C.
. The stoichiometric matrices of the original and reduced system are given by where we reduced the subnetwork γ = ({v 1 }, {e 1 }). The reduction is visually expressed as Suppose that we take the mass-action kinetics,  The system has one conserved charge and we specify the value as ℓ = x 1 + x 2 + x 3 . The steady-state reaction rates of the original system are  where K is defined by 1 K := 1 k1 + 1 k2 + 1 k3 . In the reduced system, ℓ ′ = x 2 + x 3 is a conserved charge. The steady-state rates in the reduced system are where 1 K ′ := 1 k2 + 1 k3 . Note that, if we want to have the same steady-state in the reduced system as the one in the original system, we have to choose the parameters so that ℓK = ℓ ′ K ′ . This is in contrast to Example 5, where no fine-tuning of the parameters is needed. The difference is attributed to the fact that the subnetwork γ is not a buffering structure and the index is nonzero, λ(γ) = −1 + 1 − 0 + 1 = 1.
e 3 , e 4 , e 5 , e 6 }). The stoichiometric matrix changes under reduction as where we chose the subnetwork γ = ({v 1 , v 2 }, {e 2 , e 3 }) to be reduced. The reduction is visually expressed as  The subnetwork γ is a buffering structure: λ(γ) = −5 + 7 − 2 + 0 = 0. The reduction is visually expressed as We note that, under the reduction, the stoichiometries for reactions e 3 , e 4 , e 10 are changed from the original ones.
In particular, e 10 , which is originally monomolecular, becomes non-monomolecular, 2v 6 → v 9 . To reproduce steady states of the original system, the rate r 10 is required to be the same as before the reduction; for example, in the mass-action kinetics, r 10 is given by r 10 (x 6 ) = k 10 x 6 rather than by r 10 (x 6 ) = k 10 x 2 6 , even after the reduction.

C. Reduction as a morphism of chemical reaction networks
The structure of a reduced network is characterized by the generalized Schur complement (101). Here, let us show that this form arises if we consider a map between chemical reactions that shrink a subnetwork to a point. The morphisms of chemical reaction networks have been discussed, for example, in Ref. [56]. Let us prepare some terminologies.
A degenerate reaction is a trivial reaction since it does not change anything, and the removal of degenerate reactions does not affect the chemical properties of the reaction network. A degenerate reaction is represented as a 0−column in the stoichiometric matrix.
Let us slightly extend the definition of CRNs for technical reasons.
Definition 8 (Generalized CRNs). A generalized chemical reaction network Γ is 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, Compared with the previous definition of a CRN, N is replaced with real numbers, R. We also call an element of R V as a chemical complex. In the remainder of this paper, we will mean a generalized CRN when we write a CRN. This extension is needed because the reductions we consider do not necessarily preserve the integrality of the source and target functions. However, we note that the integrality can be always recovered by reactionwise rescaling, if the original s and t functions are valued in integers.
Definition 9 (CRN morphisms). A CRN morphism ϕ from Γ = (V, E, s, t) to Γ ′ = (V ′ , E ′ , s ′ , t ′ ) is a pair of maps, (ϕ 0 , ϕ 1 ), where which we call a chemical complex map and a reaction map, respectively, such that the following diagrams commute, We introduce the matrix representation of a chemical complex map and a reaction map, On the spaces of chains, a CRN morphism induces the following commutative diagram, where ∂ ′ 1 is the boundary operator on C 1 (Γ ′ ). Namely, In terms of the matrix components, We write this relation in the matrix form, Now we are ready to discuss a morphism that corresponds to the reduction: Definition 10 (Reduction morphisms). We define a reduction morphism from Γ to Γ ′ , associated with a subnetwork γ ⊂ Γ, as a CRN morphism satisfying the following properties: 1. The chemical complexes and reactions in Γ \ γ are unchanged.
2. All the chemical complexes in γ are collapsed into one chemical complexc in Γ \ γ, in such a way that image of all the reactions in γ are degenerate in stoichiometry.
Let us here show that a reduction morphism gives rise to the reduced stoichiometric matrix given by the generalized Schur complement (101). We consider the matrix representation of a reduction morphism. From property 2 of reduction morphisms, the chemical complex map and the reaction map are both identity on v i ∈ V Γ\γ and e A ∈ E Γ\γ , Furthermore, we can always set ϕ 1 | γ = 1 without affecting the chemical properties, since degenerate reactions do nothing chemically. By arranging the rows and columns, the species and reaction maps of a reduction morphism can be written in the following form, where F is some matrix (see examples later in this section). The explicit form of F does not matter here 6 . By plugging Eq. (131) into the commutativity condition (129), we find that S ′ is written as 6 In the case of a directed graph (i.e. a monomolecular reaction network), F has one row whose elements are all 1 and other components are all zero, From the condition that the image of the reactions in γ under a reduction morphism is degenerate reactions, we have where D is a matrix satisfying DS 11 = 0. The stoichiometric matrix S ′ can be now written as The term DS 12 vanishes if and only if coker S 11 ⊂ coker S 12 .
Note that the combination S 21 S + 11 S 12 does not depend on the choice of the pseudo-inverse, as long as Eqs. (135) and (138) are satisfied. After removing degenerate reactions from Eq. (137), which does not change the chemical property of the system, we arrive at the generalized Schur complement (101) that we introduced earlier. In this way, a CRN morphism that shrinks a subnetwork gives rise to the reduced stoichiometric matrix given by the Schur complement.
What we have just shown can be summarized as the following statement: Theorem 3. Under a reduction morphism associated with γ ⊂ Γ, the stoichiometric matrix of Γ ′ can be written uniquely (up to the changes of rows and columns) in the form if and only if the following conditions are satisfied for the subnetwork γ, 7 ker S 11 ⊂ ker S 21 , coker S 11 ⊂ coker S 12 .
Conversely, for a given output-complete subnetwork such that Eq. (140) is satisfied, we can construct the following reduction map by where D is a matrix satisfying DS 11 = 0. The commutativity condition reads Since (1 − S + 11 S 11 ) is a projection matrix to ker S 11 and we have ker S 11 ⊂ ker S 21 by assumption, the matrix S 21 (1 − S + 11 S 11 ) is a zero matrix. Furthermore, DS 12 = 0 by the assumption coker S 11 ⊂ coker S 12 . Thus, we arrive at the reduced stoichiometric matrix of the form (139).
Below, let us illustrate reduction morphisms in simple examples. Example 9. Let us consider the following closed directed graph. We consider the morphism, which can be pictorially represented as 7 We note that the condition (135) is equivalent to the absence of "emergent cycles", which can be written as c(γ) = 0 in the notation of Sec. V. We show this equivalence in Sec. V below Eq. (182). The condition (138) implies the absence of "emergent conserved charges," which can be written as d(γ) = 0, but the converse is not true. We discuss more on the meaning of emergent cycles and conserved charges in Appendix B.
The reduction shrinks the vertices in γ to a single complex,c = v ′ 3 . The species and reactions are mapped as In the matrix form, Using the consistency condition (129), the stoichiometric matrix of Γ ′ can be written as This is indeed of the form (139).
. We consider the reduction of γ = ({v 1 }, {e 1 }). The corresponding reduction morphism is visualized as The chemical complex map and the reaction map are given by The action ϕ 0 (v 1 ) is determined so that the image of e 1 be degenerate in stoichiometry. The image of e 1 is the following degenerate reaction, In the matrix form, The stoichiometric matrix of Γ ′ is written as

V. REDUCTION AND BUFFERING STRUCTURES
We here explore the close connection between the structural sensitivity analysis and the reduction method we introduced in 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 be completed in Sec. V C.
The structure of this section is as follows: In Sec. V A, we show that the influence index allows for a decomposition in terms of the numbers of cycles and conserved charges. In Sec. V B, we construct a short exact sequence of the chain complexes for a subnetwork γ ⊂ Γ, under some conditions. This short exact sequence automatically derives a long exact sequence of homology groups. Using this exact sequence, we can describe the relationship among cycles and conserved charges of γ, Γ and Γ ′ . In Sec. V C, we show the main result, that is, that the steady state of the reduced network is the same as the one of the original network, under some conditions. In the proof, the long exact sequence prepared in subsection B plays an important role. In Sec. V D, we study the situation where we have nested subnetworks γ ′ ⊂ γ ⊂ Γ. In this case, we have a subnetwork γ/γ ′ ⊂ Γ/γ ′ . We will show that the reduced network Γ/γ is the same as (Γ/γ ′ )/(γ/γ ′ ). This ensures that the eventual network does not depend on the ordering of the reductions.

A. Decomposition of the influence index
As we detailed in Sec. II, steady-state properties are captured by cycles and conserved charges, which are the elements of homology groups. In this subsection, we study their meaning in more detail, and discuss the relation between the influence index λ(γ) and cycles/conserved charges in γ, Γ, and Γ ′ . We introduce a decomposition of the influence index in terms of the spaces of cycles/conserved charges of certain classes.
We first note that the index can be written as where we used Eq. (51). With the first two terms, we define c(γ) := | ker S 11 | − |(ker S) supp γ |.
The number c(γ) is a nonnegative integer, because there is an injective map from (ker S) supp γ to ker S 11 . Indeed, an element of (ker S) supp γ is written as c = c 1 0 satisfying the condition Consider an injective map c → c 1 . Equation (156) indicates that the image of this map is always included in ker S 11 , Thus, we have c(γ) ≥ 0. Now let us turn to the latter two terms in Eq. (154). Note that 8 whereP 0 γ := 1 − P 0 γ . The second term of the RHS of Eq. (159) is the number of the conserved charges of Γ supported in Γ \ γ,d 8 For a vector space V and a projection matrix P , whereP := 1 − P .
where the spaceD ′ (γ) is given byD 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, 9 where we defined Namely, we have the following decomposition of coker S, The dimension of coker S is written as where d(γ) := |D(γ)|, and d ′ (γ) := |D ′ (γ)|. We now have the expression To rewrite | coker S 11 |, we introduce the following spaces, The elements of D 11 (γ) are conserved charges in γ that can be extended to a global conserved charge, while those in D(γ) are emergent conserved charges that are only conserved in the subnetwork γ.
Observe that D(γ) ∼ = D 11 (γ). Indeed, there is a surjection The kernel of this map isD ′ (γ), and the induced map D(γ) = X(γ)/D ′ (γ) → D 11 (γ) is an isomorphism. Thus, |D 11 (γ)| = |D(γ)| = d(γ) and we have the decomposition, where d(γ) := | D(γ)| is the number of charges that cannot be obtained as the projections of conserved charges in Γ. Combining Eqs. (155), (167), and (171), we find that the influence index is written as where we defined The decomposition (172) is the central result of this subsection. Each term of Eq. (172) allows for the following intuitive interpretations: • The first term c(γ) = | ker S 11 /(ker S) supp γ | represents the number of emergent cycles in γ. Namely, c(γ) is the number 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 focusing on γ, namely those that are conserved in Γ but their projection to γ are not.
• The third term d(γ) = | D(γ)| is the number of emergent conserved charges in γ. It is the number of conserved charge in γ that cannot be extended to conserved charges in Γ. The meaning becomes evident if we note that D(γ) is isomorphic to the space that consists of d 1 ∈ coker S 11 that are orthogonal to the vectors that can be extended 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 (172) can be visually understood from the structure of A-matrices. An element of D ′ (γ) can be regarded as a conserved charge in Γ ′ via an injective mapφ 0 : D ′ (γ) → coker S ′ , which we will construct as follows. We defineφ 0 on each component of D ′ (γ) = (coker S)/X(γ) ⊕D ′ (γ). The map which is obviously injective, and is well-defined since d 2 belongs to coker S ′ by Note that the second equality follows from d T 2 S 22 = 0 and d T 2 S 21 = 0, which hold by the assumption Next, we construct an injection (coker S)/X(γ) → coker S ′ . For d 1 d 2 ∈ (coker S)/X(γ), we can always choose a Since (1 − S 11 S + 11 ) is a projection matrix to coker S 11 , we have d T 1 (1 − S 11 S + 11 ) = 0, and thus d 2 ∈ coker S ′ . This defines an injective map (coker S)/X(γ) → coker S ′ .
In general, the conserved charges in Γ ′ consists of those obtained from the conserved charges of Γ and emergent ones, where d ′ (γ) := |(coker S ′ )/imφ 0 | indicates the number of emergent conserved charges in Γ ′ .

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 steadystate characteristics of reaction networks. Suppose that we have a reaction network Γ, and choose a subnetwork γ ⊂ Γ, and reduce it to obtain Γ ′ = Γ/γ. The inter-relations of homologies of γ, Γ, and Γ ′ , can be systematically treated using a long exact sequence for a pair of chemical reaction networks, which we define momentarily. We consider the following short exact sequence of chain complexes, where the space of chains in Γ ′ is given by C n (Γ ′ ) := C n (Γ)/C n (γ). In the linear-algebra notations, the boundary maps are given by the following multiplications of matrices on vectors, We define the horizontal maps by The exactness of the rows of Eq. (177) can be checked easily. Note that ϕ is the reduction morphism (141) followed by the removal of degenerate reactions. One can check that the diagram (177) commutes when the following condition is satisfied: where c 1 ∈ C 1 (γ). The matrix (1 − S + 11 S 11 ) is the projection matrix to ker S 11 , and Eq. (181) is equivalent to This condition is the same as the condition that an arbitrary term in Eq. (98) vanishes. The condition (182) has a natural interpretation in terms of cycles: Eq. (182) is equivalent to c(γ) = 0, namely the absence of emergent cycles, which can be checked as follows. When c(γ) = | ker S 11 /(ker S) supp γ | = 0, any c 1 ∈ ker S 11 is a cycle in Γ by an inclusion to C 1 (Γ). Thus, c 1 satisfies This implies S 21 c 1 = 0 and we have ker S 11 ⊂ ker S 21 . Conversely, when ker S 11 ⊂ ker S 21 is true, the map ker S 11 ∋ c 1 → c 1 0 ∈ (ker S) supp γ is a bijection. This implies c(γ) = 0. Thus, we have shown that the diagram (177) commutes if and only if γ has no emergent cycle. Applying the snake lemma to Eq. (177), we obtain a long exact sequence, whereψ 0 andφ 0 are induced maps of ψ 0 and ϕ 0 . The map δ 1 : H 1 (Γ ′ ) → H 0 (γ) is called the connecting map. For a given c 2 ∈ H 1 (Γ ′ ), the connecting map is given by 10 where [...] means to identify the differences in im S 11 . Let us look at the consequences of the long exact sequence (184). Suppose that we choose γ so that its homology groups are trivial, Then, we have the isomorphisms, equivalently, Thus, the spaces of cycles and conserved charges before and after the reduction are isomorphic when γ has trivial homologies. Example 6 in Sec. IV corresponds to this situation, where the partial stoichiometric matrix is given by S 11 = −1 , whose kernel and cokernel are trivial. The exact sequence applies as long as the commutativity condition, ker S 11 ⊂ ker S 21 , is satisfied, and we can consider more general cases with ker S 11 = 0. If the connecting map δ 1 : H 1 (Γ ′ ) → H 0 (γ) is a zero map, the long exact sequence (184) results in the following two exact sequences, This implies the isomorphisms, Note that ker S 11 consists of only locally supported global cycles, due to the assumption ker S 11 ⊂ ker S 21 . The isomorphisms (191) represent equivalence of chemical reaction networks up to locally supported global cycles and locally supported global conserved charges (emergent conserved charges are also absent when δ 1 is a zero map, as we see below). Let us examine the condition when the connecting map is a zero map. δ 1 is a zero map if for any c 2 ∈ H 1 (Γ ′ ). Below we show that, if every conserved charge in γ is obtained by the projection of a global conserved charge in Γ (namely, when there is no emergent conserved charge), the connecting map δ 1 is a zero map. For a given d 1 ∈ coker S 11 , there exists an element of coker S, where we used d T 1 S 11 = 0. Let us pick c 2 ∈ H 1 (Γ ′ ) = ker S ′ . The quantity d T 1 S 12 c 2 can be shown to vanish as follows: Therefore, we have shown d T 1 S 12 c 2 = 0 for any d 1 ∈ coker S 11 and c 2 ∈ ker S ′ . This is equivalent to S 12 c 2 ∈ (coker S 11 ) ⊥ .
The relation between the long exact sequence and the numbers of cycles and conserved charges of various types is summarized in Fig. 4. The vertical lines represent the spaces, and the kernels are shown in black. Since it is an exact 10 The connecting map is identified as follows. An element c 2 ∈ H 1 (Γ ′ ), can be included in C 1 (Γ ′ ). sequence, the kernel and image coincide at each space, such as im ψ 1 = ker ϕ 1 and so on. The exactness is the key to the connections between cycles and conserved charges of particular types. Let us see an example. The image of δ 1 is the space of emergent conserved charges, im δ 1 = D(γ). They are emergent, because the image of δ 1 is the kernel ofψ 0 , and there is no counterpart in Γ. The connecting map δ 1 provides us with a one-to-one mapping between an emergent cycle in Γ ′ and an emergent conserved charge in γ (elements of ker δ 1 are not emergent, since they can be written as an image of ϕ 1 due to the exactness). The numbers d(γ), d ′ (γ) in Fig. 4 are the same as the dimensions of the spaces (162) and (163) that we defined previously. Compare Fig. 4 also with Fig. 7 in the Appendix B 2, where we discuss the relation between the numbers of cycles and conserved charges and the structure of the A-matrix. The long exact sequence is valid when c(γ) = 0 (i.e., when the diagram (177) commutes). This implies d ′ (γ) = 0 and there is not emergent conserved charge in Γ ′ , sinceφ 0 is surjective.

C. Reduction of buffering structures
Here we present the main result, regarding the reduction of buffering structures. The following theorem represents a particularly nice property of buffering structures under reductions. We show that the steady-state concentrations and rates of the network obtained by reducing a buffering structure are exactly the same as those of the network before reduction, without any modification of parameters. Thus, the reduction of a buffering structure preserves the steady-state properties of the boundary degrees of freedom. The theorem only relies on topological information of the network and is true regardless of the kinetics. Theorem 4. Let Γ be a regular chemical reaction network with kinetics r(x) and let γ be an output-complete subnetwork of Γ. We assume that the subnetwork γ does not have an emergent conserved charge. We consider a reduced network Γ ′ = Γ/γ. If γ is a buffering structure, we have the isomorphisms, ker S/ ker S 11 ∼ = ker S ′ , coker S/ coker S 11 ∼ = coker S ′ , Furthermore, when (r, x) is steady-state reaction rates concentrations of Γ, whose components we separate into those in γ and Γ \ γ as then, (r 2 , x 2 ) is a steady-state solution of Γ ′ .
Remark 14. Let us comment on the assumption of the absence of emergent conserved charges. Under the assumption of 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 d(γ) = 0 for a connected and output-complete subnetwork γ, assuming that Γ is regular (see Appendix C 2), and this condition is redundant. So far, the examples of buffering structures with nonzero emergent conserved charges are pathological in some sense. Presently, we have not been able to prove the absence of emergent conserved charges for a generic (sound) reaction network, and thus it is assumed. We have more discussions on this point in Appendix C.
Remark 15. We note that there is a possibility that the reduced system might have some solutions which are not allowed in the original system, depending on the kinetics 11 . This can occur when the reactions in a subnetwork have limitations in the values of reaction rates. When such a subnetwork is removed, the reduced system does not have the restrictions, and there may appear additional solutions. Let us illustrate this with an example. We consider a reaction network Γ = ({v 1 , v 2 , v 3 }, {e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 }) given by the following set of reactions, Let us here choose the kinetics as For reaction r 1 , we adopted the Michaelis-Menten kinetics, and we chose the mass-action kinetics for other reactions. The rate equations read d dt where we reparametrized the equation for d dt x 3 using d 1 , d 2 , and d 3 such that d 1 < d 2 < d 3 . From Eq. (202), we get two candidates of stable steady-state values,x 3 = d 1 , d 3 (note that d 2 is unstable). However, those candidates may not lead to the solutions of whole equations when the reaction rate r 1 has a bound, as in the current example. Using Eq. (200), the steady-state value of x 1 is given byx If the denominator of Eq. (203) is negative, it is not a valid solution. Thus, depending on the values of d 1 and d 3 , the original network may have no, or one, or two solutions. The subnetwork γ = {(v 1 ), (e 1 )} is a buffering structure, and we can consider the corresponding reduced network Γ ′ = Γ/γ. In the reduced network, such a restriction on the values of (internal) reaction rates is invisible. Hence, the reduced system may admit more solutions that were not possible in the original system.
Proof. The regularity of Γ requires λ(γ) ≥ 0 (Remark 8). In the absence of the emergent conserved charges, we have Since c(γ) and d l (γ) are nonnegative integers, we have c(γ) = 0 and d l (γ) = 0. Since c(γ) = 0, we can use the long exact sequence (184). Because there is no emergent conserved charge, d(γ) = 0, by assumption, the connecting map δ 1 is a zero map in the long exact sequence. This proves Eq. (196). Let us proceed to the latter part of the claim. The steady-state condition of Γ is written as dᾱ · x = ℓᾱ.
As usual, we divide the degrees of freedom to those in γ and Γ \ γ. Then Eq. (205) is written as The reactions r 2 (x 2 ) depend only on x 2 , because γ is chosen to be output-complete. The first equation can be solved for r 1 as r 1 = −S + 11 S 12 r 2 + c 11 , with c 11 ∈ ker S 11 , and we have where the last equality is due to c(γ) = 0, that is equivalent to ker S 11 ⊂ ker S 21 .
Let us turn to the conserved charges. Recall that d l (γ) is written as d l (γ) = |(coker S)/X(γ)|. Because of the decomposition (165), when d l (γ) = 0, the space coker S is written as the direct sum of D(γ) andD ′ (γ), Correspondingly, we can divide the basis vectors of coker S into two classes, {dᾱ} = {dᾱ γ , dᾱ ′ }, where {dᾱ γ } is a basis of D(γ), and {dᾱ ′ } is a basis ofD ′ (γ). The basis vectors are of the form, With this basis of coker S, Eq. (206) is written as In fact, dᾱ ′ 2 is a conserved charge in Γ ′ , dᾱ ′ 2 ∈ coker S ′ , as we see in the following. Since dᾱ ′ ∈ coker S , it satisfies This implies that dᾱ hence dᾱ ′ 2 ∈ coker S ′ . Thus we have obtained an injective map, This map is nothing but the induced mapφ 0 . It is important to note that, when c(γ) = 0, this map is a surjection, that is evident from the long exact sequence (184) 12 . The equations satisfied by the boundary part (denoted by 2) of the concentrations/rates of Γ are Eqs. (208) and (212). Since all the conserved charges in Γ ′ is given as a imagē ϕ 0 , we find that the set of Eqs. (208) and (212) are exactly the same as the steady-state condition for the reduced network Γ ′ , where x ′ = x 2 and r ′ (x ′ ) = r 2 (x 2 ). Thus, the steady-state solution of Γ ′ 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 γ ′ ⊂ γ ⊂ Γ. Given the stoichiometric matrix S of the whole network, we denote the stoichiometric matrices of the subnetworks γ and γ ′ by S γ and S γ ′ , respectively. The submatrices are included in the following form, where * indicates an arbitrary matrix. Let us consider the situation where γ and γ ′ has no emergent cycle and emergent conserved charge in Γ, namely c(γ) = d(γ) = 0, and c(γ ′ ) = d(γ ′ ) = 0. Under those assumptions, the quotient formula of the generalized Schur complement [57] holds, This indicates the isomorphisms of homology groups, for n = 0, 1. Thus, when we perform the reductions of nested subnetworks that have no emergent cycles and emergent conserved charges, the order of the reduction of them does not matter.

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 stoichiometric matrix presented in Ref. [45], which is constructed based on Ref. [4] with minor modifications. The network structure is shown in Fig. 5a, which consists of the glycolysis, the pentose phosphate pathway (PPP), and the tricarboxylic acid cycle (TCAC). The list of the reactions for this system is given in Appendix D 1. Here, we assume that H 2 O and cofactors such as ATP and NADH are abundant and do not affect the behavior of the system. Buffering structures in this network have been identified in Ref. [24] and there are in total 17 buffering structures, which we list in Appendix D 2. As we showed in Sec. III C, the intersections or unions of buffering structures are also buffering structures. They form a hierarchy, and such an architecture can be regarded as a source of robustness against perturbations, since buffering structures 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 as the original network as we showed in Sec. V C. We denote the whole network by Γ. We can pick a buffering structure γ 8 , which is a part of the pentose phosphate pathway (the yellow subnetwork in Fig. 5a) and given by 13 and perform a reduction to obtain Γ 1 := Γ/γ 8 . The stoichiometric matrix of the reduced reaction network can be computed by Eq. (101). The resulting network is shown in Fig. 5b. The reduction procedure induces rewiring of the reactions, which are colored in magenta in Fig Fig. 5a). This reduction is the same as reducing γ 5 ∪ γ 14 from Γ. The result of the reduction is shown in Fig. 5c. Again, rewiring occurs and the reactions 5,6,15, and 36 are modified from the original system. Finally, let us focus on the part colored in green in Fig. 5a, which consists of the following subsets of chemical species and reactions, The complement of the subset (222) is given by γ 5 ∪ γ 7 ∪ γ 14 , which hence is a buffering structure, and a reduction can be performed. The structure of the reduced network Γ 3 = Γ/(γ 5 ∪ γ 7 ∪ γ 14 ) is given in Fig. 5d. Compared to the original 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. 5 (the original network Γ and the reduced ones, Γ 1 , Γ 2 , Γ 3 ), using the same initial condition and reaction rate constants in all of the four cases (see Appendix D 3 for details of parameter values). The time-series of concentrations are presented in Fig. 6. After the initial transient dynamics, the original system approaches a (stable) steady state [ Fig. 6(a)]. We can see that the reduced systems can reproduce the steady-state concentrations that the original system eventually reaches, although they have distinct short-time dynamics [ Fig. 6(b-d)].
In this way, buffering structures work as a guide as to how to perform the reduction and simplify a complex reaction network. As long as the reduced part is a buffering structure and we use the generalized Schur complement (101) as a stoichiometric matrix of a reduced network, the steady-state concentrations and rates of the remaining part stay the same as the original ones regardless of the details of the kinetics, as a consequence of Theorem 4.

VII. SUMMARY AND OUTLOOK
The main focus of the present paper was the relationship between the structure and functions of the chemical reaction network. As a characterization of the structure, homology and cohomology groups for chemical reaction networks were introduced, in which the actions of boundary and coboundary operators are determined by the stoichiometry. The elements of homology groups correspond to cycles and conserved charges of chemical reaction networks, and steady states were shown to be determined by the elements of the cohomology groups. In a similar way to the homology and cohomology groups of topological spaces, the Mayer-Vietoris sequence and the long exact sequence of a pair of chemical reaction networks were introduced, the latter being particularly useful for studying the reduction of reaction networks.
We propose a method of reduction of chemical reaction networks. The reduced network is characterized by the stoichiometric matrix obtained by eliminating the chemical species and reactions of an output-complete subnetwork via the Schur complementation. The reduction relies only on the stoichiometry, which determines the topology of the reaction networks, and thus is applicable to any kind of kinetics. This represents an advantage since in many biological systems it is difficult to experimentally determine the kinetics and parameters of the reactions. For tracking the change of cycles and conserved charges under the reductions, the tools of algebraic topology, such as the long exact sequence, have been useful. We have studied how the law of localization can be understood from this perspective. We showed that the influence index is expressed in terms of the numbers of cycles/conserved charges of particular types, as in Eq. (172). We also showed that the influence index is a submodular function over output-complete subnetworks. A corollary of this is that buffering structures are closed under intersection and union, which is useful when we enumerate the buffering structures of a large reaction network. As a central result of the paper, we showed that buffering structures, which are subnetworks with vanishing influence index, behave nicely under the reduction. Namely, under the reduction of a buffering structure, the steady state of the remaining elements of the network stays the same as the original network (Theorem 4). The theorem justifies the intuition that buffering structures are regarded as 'irrelevant' substructures: they can be safely eliminated through the reduction method proposed here without changing the long-time behavior of the system. The reduction procedure introduces rewiring of reactions, which is necessary so that the steady state is not modified under the reduction. As an application of the reduction method, we discussed the reduction of the central metabolic pathway of E. coli and illustrated that reactions are rewired non-trivially under the reduction. We also demonstrated the invariance of the steady state under the reduction of buffering structures by numerically solving the rate equations before and after the reduction 14 .
Our results highlight that special care should be taken when simplifying a reaction network. A naive elimination of a subnetwork not of interest would alter steady-state properties of the original system. As long as the subnetwork has the vanishing influence index and reactions are rewired appropriately using the generalized Schur complement, it can be eliminated while keeping the steady state intact.
Another significance of our method is that it allows us to identify the modules in a complex network and facilitates the biological interpretation of the whole system. For example, the central metabolic pathway of E. coli consists of three modules; glycolysis, TCAC, and PPP. Interestingly, the reduced network in Fig. 5d roughly corresponds to the glycolysis. The fact of glycolysis being a reduced network may suggest that E. coli can control the glycolysis in an isolated manner, and the expression levels of enzymes in the TCAC and the PPP do not affect the physiological states of the glycolysis.
For practical applications, one important issue is how to find the buffering structures efficiently in large-scale reaction networks. Although we defer this as a future problem, let us make some comments on this point. One  Establishing a combinatorial method for identifying buffering structures is an amusing problem. We believe that the basic properties of buffering structures that we showed in this paper would be useful for this purpose. For example, if a network contains many small buffering structures, we can use the reduction method repeatedly and make the network smaller one we fine a small buffering structure. This procedure is possible because the order of reduction does not matter for the buffering structures, as we showed in Sec. V D. The submodular property of the influence index and the subsequent closure property of buffering structures under unions/intersections would also be useful in enumerating buffering structures.
We believe that the mathematical formulation that we used to characterize the topology of chemical reaction networks will be useful for understanding the static and dynamical properties 15 of reaction systems. The eigenvalues of the Laplacian operators entail the information of the topology of the network connectivity. Steady states correspond to the eigenvectors with zero eigenvalues and they incorporate the 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 to approximate the dynamics as well as the steady states.
for a 0 ∈ C 0 (Γ) and a 1 ∈ C 1 (Γ). Those are generalizations of the graph Laplacian to hypergraphs. The properties of hypergraph Laplacians were discussed recently in Refs. [58][59][60]. When all the reactions are monomolecular, the Laplacian reduces to the graph Laplacian of the directed graph.
The space C 1 (Γ) admits the following orthogonal decomposition, This is a natural generalization of the Hodge decomposition of flows on networks [61] to the case of a hypergraph. Thus, given a 1-cochain f ∈ C 1 (Γ), we can decompose it in a unique way as where c ∈ ker d † 0 ∩ ker d 1 is a harmonic cochain and a ∈ C 0 (Γ). This is the Hodge decomposition associated with the complex (27). By acting d † 0 on Eq. (A4), we have We can solve this for the potential a as Here, where M + indicates the Moore-Penrose inverse of a matrix M , and a 0 ∈ ker ∆ 0 . The harmonic component c can be obtained by Using the properties of the Moore-Penrose inverse, the action of the operator that appears on the RHS of Eq. (A7) is written as for an arbitrary b 1 ∈ C 1 (Γ). The matrix 1 − S + S is the projection matrix to ker S. Thus, the harmonic component can be identified by the projection to ker S, This is consistent with the fact that c ∈ H 1 (Γ) = ker S. The potential a can be obtained by the multiplication of the Moore-Penrose inverse of S to f , where a 0 ∈ ker ∆ 0 . Let us here discuss intuitive interpretations of the integers appearing in the decomposition (172). We will refer to the elements of ker S as "global cycles" and those of coker S as "global conserved charges." Similarly, the elements ker S 11 and coker S 11 are referred to as "local cycles" and "local conserved charges." With this terminology, the elements of (ker S) supp γ are called as "locally supported global cycles. " We first look at c(γ). We denote the space by C(γ) := ker S 11 /(ker S) supp γ .
Then, c(γ) = | C(γ)|. To clarify its meaning, we represent the space (B1) as follows, An element of C(γ) is a local cycle that is not a global cycle, by which we mean that c ∈ C(γ) has its boundary in Γ \ γ. For example, let us take the subnetwork γ = ({v 2 }, {e 1 , e 2 }) of a monomolecular reaction network, Although e 1 + e 2 ∈ C 1 (Γ) has its support in γ, its boundary, 16 is outside of γ. The element e 1 +e 2 is a local cycle, since ∂ 1 (e 1 +e 2 ) is zero as a relative chain in C 0 (γ) = C 0 (Γ)/C 0 (Γ\γ). Note that the network (B3) as a whole does not have a cycle and ker S = 0. Thus, we can identify c ∈ C(γ) to be a local cycle whose boundary is out of γ. When c is viewed in Γ, it may be extended to a global cycle, but it does not have to be. Considering its meaning, we will refer to the elements of C(γ) as emergent cycles, which only appear when we focus on a subnetwork. Let us illustrate the space C(γ) pictorially. The matrix S works as a boundary operator on the space of chemical reactions. Thus, the kernel of S are linear combinations of reactions without boundaries. Cycles and non-cycles can be drawn pictorially as Cycles Non-cycles , where the boundary of the box is identified.
We consider an output-complete subnetwork γ. The space ker S 11 is spanned by local cycles in γ, for example, where the inner box represents a subnetwork γ and the red lines constitute the basis of the space. Here, the symbol means that the cut ends are reactions and not chemical species. See the following two choices for example: γ 2 16 Recall that the boundary of each reaction is specified by the stoichiometric matrix as For the left one, cut ends are both reactions. For the right one, the cut ends are a species and a reaction. Both ends have to be reactions so that the cut cycle can be a local cycle. An element of ker S 11 may be extended to a global cycle, or it may be a part of a global noncycle. The space (kerS) supp γ for the same configuration takes into account only the global cycles supported on γ, Therefore, the coset space is generated by the following elements, As we see in the figure, the space ker S 11 /(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 as a boundary operator acting on C 0 (Γ), which is the space of chemical species. In this sense, an element of coker S has no 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 0 γ (coker S) is the projection of coker S to γ, Here, how the conserved charges are cut does not matter. The space coker S 11 is generated by the red and green elements in the following figure, where filled and open rectangles mean boundaries with chemical species and reactions, respectively. The parts we denoted by green lines, , are emergent conserved charges, which are conserved when it is seen in a subnetwork but 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 system to the perturbation of parameters is not well-defined. 17 For example, one can consider the following network and subnetwork, The whole network does not have a conserved charge, but the subnetwork γ has one, v 1 + v 2 + v 3 . However, such a reaction network cannot reach a steady state, since the concentration of v 3 continues to increase. When we take the difference |P 0 γ (coker S)| − | coker S 11 | , we can count the number of lost conserved charges minus the number of emergent conserved charges (if any), where the part colored in red in the first term indicates lost conserved charges, that are conserved in Γ but their projections to γ are not. This equation is equal to the latter two terms of the decomposition (172), d l (γ) − d(γ).
Example 11. Consider a monomolecular network Γ = (V, We take a subnetwork γ 1 = ({v 2 }, {e 1 , e 2 }) that is indicated by a box. The whole network does not have a cycle, and the subnetwork γ 1 has one emergent cycle given by c = e 1 + e 2 . Also, Γ has one conserved charge, d = v 1 + v 2 + v 3 . Its projection to γ 1 is given by v 2 and it is not a conserved charge in γ 1 . So we have one lost conserved charge. Each integer appearing in the decomposition of λ(γ 1 ) is and λ(γ 1 ) = 2. For the same Γ, let us consider a different choice of a subnetwork, The subnetwork γ 2 does not have a cycle, and there is one lost conserved charge, so we have Example 12. Consider a network (V, E) = ({v 1 , v 2 , v 3 }, {e 1 , e 2 , e 3 }) with the following structure, 17 We discuss more on this point in Appendix C.
The subnetwork γ has one emergent cycle and one lost conserved charge and the influence index is 2. It is useful to look at the A-matrix to visualize the relations among cycles/conserved charges of various types in subnetworks and reduced networks.

Embedding of A-matrices
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 for conserved charges. The numbers with a prime are associated with Γ ′ . The relevant numbers are listed as follows: • v, v ′ : number of chemical species in γ, Γ ′ • c, c ′ : number of cycles of Γ, whose projections to γ, Γ ′ are also cycles of γ, Γ ′ • d, d ′ : number of conserved charges of Γ, whose projections to γ, Γ ′ are also conserved in γ, Γ ′ •c ′ ,d ′ : number of cycles/conserved charges of Γ that are locally supported in Γ ′ • c, d: number of cycles/conserved charges of Γ that have nonzero support in γ In Fig. 7, we illustrate a more detailed structure of the matrix A than Fig. 3. In the center is the matrix A of the total system Γ. We choose an output-complete subnetwork γ, and bring the rows/columns related to γ to the upper-left part. Then the matrix A looks like one in the center of Fig. 7. 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. 7. The part surrounded by a pink 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 c and d. Because the matrix A γ is square, we have the relation, This equation is in fact the same as Eq. (51). Similarly, we can consider the matrix A for the network Γ ′ = Γ/γ 18 obtained by reducing γ from Γ. The numbers of the emergent cycles and emergent conserved charges in Γ ′ are denoted by c ′ and d ′ . The matrix A Γ ′ is also square and we have The influence index is given by which is the same the decomposition (172). We can also consider a similar quantity that measure the non-squareness of the lower-right part, where the second expression is obtained using Eq. (B9). In fact, due to the squareness of the whole matrix A, λ ′ is equal to the influence index, λ = λ ′ . This results in the following relation, 1. Systems with emergent conserved charges As far as we observe, the chemical reaction systems with emergent conserved charges in output-complete subnetworks are pathological, in either of the following senses: (A) The steady-state condition, does not fully determine the steady-state solution, and arbitrary parameters have to be introduced to specify the 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 d dt The matrix A for this system is This is not invertible. The steady-state solution for the mass-action kinetics is given by where m is an arbitrary parameter. If we choose a subnetwork γ = {x}, there is an emergent conserved charge. Let us consider the fluctuations around the steady state, d dt δx δy = 1 0 −1 −1 k 1 δy k 2 δy = k 1 δy −(k 1 + k 2 )δy , where δx(t) := x(t) −x indicates the fluctuation from the steady state. The fluctuation associated with the emergent charge, δ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 By taking the derivative of those equations with respect to m, we find This means that A has a null vector and det A = 0. An example of pattern (B) is given by the following reaction system with the mass-action kinetics, The steady-state solution does not exist in general. We need to fine-tune the parameters to have a solution. When k 3 = k 4 is satisfied, we have a steady stater = k 3 1 1 1 1 T .
The matrix A is 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 1 , v 2 }, {e 1 }). The matrix A for the subnetwork is The subnetwork has an emergent conserved charge, d T 1 = 1 −1 . The time derivative of this charge is Although d T S = 0, where d T := ( d T 1 a), for any parameter a, when the steady state exists, k 3 = k 4 , the combination 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, Here we discuss an example that has a subnetwork with vanishing influence index and also has an emergent conserved charge, while the kinetics is unphysical. The rate equation of this system is We have added catalytic dependencies in the reactions r 2 (x 2 , x 4 ) and r 3 (x 3 , x 4 ). The stoichiometric matrix has a trivial kernel andr A = 0 at the steady state. The cokernel of S is also trivial.
Let us consider an output-complete subnetwork γ = ({x 1 , x 2 }, {r 1 , r 2 }). The index of γ is zero, and hence it is a buffering structure. The matrix A of the local system reads Although γ is a buffering structure, the subnetwork γ has one emergent cycle and one emergent conserved charge. For the mass-action kinetics, the steady-state concentrations are where m 1 , m 2 , m ′ 1 , m ′ 2 are arbitrary parameters. With this kinetics, det A = 0. Let us instead employ the following kinetics, where all the concentrations vanish at the steady state,x i = 0. The matrix A is now invertible, det A = k 1 k 2 k 3 k 4 = 0. Although A is regular, the sensitivity is trivial, ∂ Axi = 0, since ∂ ArB = 0 at the steady state. The kinetics (C23) is not physically sound, because the reaction r 2 (x 2 , x 4 ) can be nonzero even if the concentration of the reactant x 4 is zero (note that x 2 is catalytic). The same is true for r 3 . Let us consider the fluctuations around the steady state of the emergent conserved charge.
Hence, there is a zero mode when r 2,2 = 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, according to the separation x = x 1 x 2 . The linear fluctuations around the steady state satisfy the following equations of motion, d dt δx 1 δx 2 = SRδx = S 11 R 11 δx 1 + (S 11 R 12 + S 12 R 22 )δx 2 S 21 R 11 δx 1 + (S 21 R 12 + S 22 where we have also separated the stoichiometric matrix into submatrices, and we used R 21 = 0, which follows from the output-completeness. We consider the fluctuation associated with the emergent conserved charge d 1 , Since it is an emergent charge, d T 1 S 12 = 0. The time derivative of the RHS of Eq. (C28) reads Therefore, an emergent conserved charge results in a zero mode when d T 1 S 12 R 22 S 21 R 11 vanishes. We are not aware of a physical example in which d T 1 S 12 R 22 S 21 R 11 does not vanish. In the example given by Eq. (C16), the first term of Eq. (C29) is computed as (C30) When this vanishes, the matrix A acquires a zero mode and is not invertible.

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 in an output-complete subnetwork, the index λ(γ) is necessarily negative. If the index is negative in an output-complete subnetwork, the matrix A is not invertible, and the response of the system to the parameter-perturbation is not well-defined. We here show the following statement: Theorem 5. Suppose that γ is a connected and output-complete subnetwork of a monomolecular reaction network Γ. If d(γ) > 0, then the influence index λ(γ) is negative.
Proof. To have an emergent conserved charge in γ, all the boundaries of γ should be chemical species and not reactions, in a monomolecular reaction network. Then, all the reactions in γ should end on the chemical species inside γ, which means S 21 = 0.
Recall that an emergent cycle is c 1 ∈ ker S 11 which is not a cycle of the whole network, When S 21 = 0, there is no such c 1 meaning that all the local cycles are also a global cycle. Namely, for a given c 1 ∈ ker S 11 , c 1 0 ∈ ker S always holds. Thus, we have c(γ) = 0.

Appendix D: Metabolic pathways of E. coli
We here provide the details of the metabolic pathways discussed in Sec. VI.