Sensitivity of chemical reaction networks:
a structural approach.
3. Regular multimolecular systems
Abstract
We present a systematic mathematical analysis of the qualitative steadystate response to rate perturbations in large classes of reaction networks. This includes multimolecular reactions and allows for catalysis, enzymatic reactions, multiple reaction products, nonmonotone rate functions, and nonclosed autonomous systems. Our structural sensitivity analysis is based on the stoichiometry of the reaction network, only. It does not require numerical data on reaction rates. Instead, we impose mild and generic nondegeneracy conditions of algebraic type. From the structural data, only, we derive which steadystate concentrations are sensitive to, and hence influenced by, changes of any particular reaction rate – and which are not. We also establish transitivity properties for influences involving rate perturbations. This allows us to derive an influence graph which globally summarizes the influence pattern of any given network. The influence graph allows the computational, but meaningful, automatic identification of functional subunits in general networks, which hierarchically influence each other. We illustrate our results for several variants of the glycolytic citric acid cycle. Biological applications include enzyme knockout experiments, and metabolic control.
*
Institut für Mathematik
Freie Universität Berlin
Arnimallee 7
14195 Berlin, Germany
Contents
1 Introduction
For large classes of biological, chemical or metabolic reaction networks, detailed numerical data on reaction rates are neither available, nor accessible, by parameter identification. See the large, and growing, data bases on chemical and metabolic pathways like [Le+06, KG00], for thousands of examples. One standard approach to establish, and check, the validity of such networks are knockout experiments: some reaction is obstructed, via the knockout of its catalyzing enzyme, and the response of the network is measured, e.g., in terms of concentration changes of metabolites. The large area of metabolic control studies how reaction rates steer the network to desired behaviour, or switch between different tasks; see for example [HS96, Fel92, Ste84] and the references there.
In this setting it is our goal to develop a reliable, and largely automatic, mathematical tool to aid our systematic understanding of large networks. Chemical reaction networks consist of reactions among metabolites. We consider reactions, quite generally, to cover biological and chemical reactions. Since we do not incorporate temperature dependence, explicitly, our setting is isothermal. The reacting chemical species, often called metabolites in biological settings, are denoted by labels .
More specifically, we address the response of steady states to rate perturbations in the network. Here steady state refers to any timeindependent longterm state of the system. “Longterm” refers to the relevant time scales of the model. Timeperiodic or chaotic responses are excluded, at present.
In experiments, only those steady states may be observable which are stable, or at least metastable on the relevant time scale. Our mathematical approach is not limited by any stability or hyperbolicity requirements other than some mild nondegeneracy assumption. For simplicity we present our approach in a local setting of linearized steady state response to small perturbations. In the concluding discussion we indicate how our results extend to the global setting of knockout experiments.
We derive a qualitative sensitivity matrix for the steady state response, which we encode as a flux influence relation influences ; in symbols: . Here indicates the label of the perturbed reaction rate, given by the experimental setting, and indicates a nonzero resulting flux change of the reaction with label in the network. The flux of measures the rate of conversion from input metabolites of reaction to output metabolites. The flux change measures the change of that conversion rate, at steady state, under the influence of the external rate perturbation of reaction . Note how the roles of and in the influence relation are subtly different. For example, may, or may not, influence itself.
Although flux changes are the more convenient object, mathematically, they are less accessible in experiments. We therefore include any nonzero resulting concentration change of any metabolite , at steady state, into our analysis; in symbols: . To combine both the flux and the concentration aspect, notationally, we define the influence relation
(1.1) 
to indicate that the effect of the perturbation on is a nonzero change of . Here denotes either a reaction or a metabolite . In this case of a nonzero influence we also say that the reaction or metabolite is sensitive to .
The influence relation (1.1), in itself, does not carry much information beyond a casuistic tabulation of whoinfluenceswho. It is a transitivity property, which makes flux influence a central tool for the understanding of sensitivity results in metabolic networks. Consider any 2step chain of influences , where denote reaction labels in the network. As a consequence, direct influence will be established. Only this transitivity of influence,
(1.2) 
justifies the notion of a hierarchy of influence.
It is the notion of nonzero influence, together with this central transitivity property of the flux influence relation, which will identify meaningful units in a network and will allow us to sum up all results of rate perturbation experiments in a single flux influence graph below.
One first phenomenological indication for the mathematical structure of flux influence, which motivated our detailed study, was the apparent sparsity of sensitivity: Given a rate perturbation of a specific reaction , many are not influenced at all, because the reaction flux or metabolite concentration associated to does not change. Such zero response is the counterpart of the flux influence graph . In fact, any zero response is a rational and mathematically rigorous test to any purported pathway structure: Any experimentally validated nonzero response, above error threshold, which contradicts a zero entry in the sensitivity matrix of influences, falsifies the underlying pathway.
We briefly comment on previous work in this series of papers, even though the present paper is selfcontained and does not require any but the most standard mathematical prerequisites. The present paper is a sequel to [MF15, FM15]. In [MF15] a hierarchy of influence patterns was observed in several examples, on a purely phenomenological level and without deeper mathematical justification. Most notably, a simplified version of a model by Ishii and coworkers, [Ish+07], on the tricarboxylic citric acid cycle (TCAC) was studied. The paper [MF15] concluded:
… our explanation of hierarchy patterns is still intuitive. Only in the monomolecular case, so far, are we able to prove the observed transitivity of the influence relation perturbation of reaction implies flux change in reaction , . This transitivity underlies the hierarchy of flux response patterns and their relation to directed cycles in the reaction network. See [FM15] for the mathematical details which are beyond the scope of our present exposition.
The mathematical companion paper [FM15], however, did not yet meet our mathematical objectives. It served as just a first attempt towards a comprehensive mathematical understanding of the sparse anecdotal evidence compiled in [MF15]. Mathematical results were limited to the monomolecular case, where each reaction just converts a single metabolite. Bimolecular reactions were excluded, as were reactions with more than a single metabolite as a product. Although of substantial mathematical interest, the restriction to monomolecular reactions still excluded most biologically relevant examples. In particular, the original challenge by the TCA cycles of Ishii and Nakahigashi, [Ish+07, Nak+09], had not yet been met.
Only in the present paper, for the first time, are we able to present a comprehensive and mathematically founded theory which explains, and proves, the appearance of sensitivity patterns and hierarchies of steady state responses in multimolecular metabolic networks. The full mathematical background of a theoretical example from [MF15] is developed in section 3 below. Three variants of the Ishii TCAC metabolism [Ish+07] will be compared in section 8, along with a further modification by Nakahigashi and coworkers, [Nak+09], in the light of our present mathematical understanding.
In section 2 we present our main mathematical results. To prepare, we describe our general mathematical setting next. We strongly recommend [Fei95, Pal06] for a general background.
Let the labels enumerate the total number of reactions in the network. Let denote the set of all reactions. Similarly, let the labels enumerate the total number of metabolites in the network. Let denote the set of all metabolites. Mimicking standard chemical notation, a reaction network is given by reactions
(1.3) 
The components of the stoichiometric coefficient vectors are nonnegative, and usually integer. Reversible reactions are, not required but, admitted and will be listed as two separate reactions. For notational convenience we will frequently identify with its index set . We will even write , instead of , to emphasize the supporting component set of vectors. For subsets we analogously abbreviate
(1.4) 
to denote the subspace of vectors which are supported on , only. We call an input, reactant, or educt of reaction , if ; in symbols:
(1.5) 
We use this notation even when catalyzes , i.e. when , and in presence of further inputs of reaction . Outputs, or products are given by . Feed reactions only depend on external chemical input species which are provided at constant concentration levels, unaffected by the network. Thus feed reactions have . Exit reactions, in contrast, only produce chemical output species which do not reenter the network. Thus exit reactions have . See section 3 for an easy example.
The stoichiometric matrix is essential to derive the differential equation (1.10), below, for changes of metabolite concentrations in the network (1.3). The matrix is defined by
(1.6)  
where , alias , denotes the th unit vector. In other words, the columns of the stoichiometric matrix are simply the differences of the stoichiometric output and input vectors and . The usually nonlinear reaction rate , at which reaction occurs per time unit, depends on the concentrations of the metabolites . Evidently only depends on the input metabolites of reaction . In symbols, the partial derivatives of the reaction rates satisfy
(1.7) 
In other words, the supports satisfy
(1.8) 
for each reaction , by our convenient abuse (1.4) of notation. The time evolution of the metabolite concentrations is then given by the coupled nonlinear ODE system
(1.9) 
for , under isothermal conditions for the reaction rates. In vector notation with and this reads
(1.10) 
For simplicity of presentation, we will assume
(1.11) 
throughout our paper. This excludes cokernel of , alias stoichiometric subspaces. These are affine linear subspaces which are time invariant under the flow of . In other words, we exclude trivially conserved linear combinations of the metabolite concentrations . For further comments on the case of stoichiometric subspaces we refer to the discussion in section 9.
Steady states are timeindependent solutions of (1.10), i.e.
(1.12) 
Considerable effort has gone into the study of existence, uniqueness, and possible multiplicity of steady state solutions of (1.12). See in particular Feinberg’s pioneering work [Fei95], his rather advanced recent results [SF13], and the many references there. In the present paper we do not address this question, at all. Rather, we assume existence of a steady state throughout this paper. We neither assume uniqueness, nor stability, of this steady state. We focus on the following central question:
(1.13)  
To define and study this response sensitivity of steady states with respect to small external changes of any reaction rate we introduce formal reaction rate parameters as
(1.14) 
We also write on the right hand side of (1.12). For extreme generality, we could choose to parametrize the functions , themselves. More modestly, we might think of just a rate coefficient like := .
The standard implicit function theorem in a setting immediately answers the response question – albeit rather abstractly and under mild nondegeneracy assumptions. Given a reference steady state solution , for , we obtain a unique differentiable family of solutions , for sufficiently small , such that
(1.15) 
holds for the partial derivatives. Here the rate matrix is the Jacobian matrix of the rate function vector with respect to ; see (1.7) for the definition of the partial derivatives of . For the partial derivatives of the rate function vector with respect to the parameters it is sufficient to consider unit vectors
(1.16) 
all other cases can be derived from that. We call
(1.17) 
the sensitivity of metabolite with respect to rate changes of . Similarly, we call
(1.18) 
the sensitivity of reaction flux with respect to rate changes of . The Kronecker symbol accounts for the externally forced flux change by variation of the parameter itself. The reactioninduced term may, or may not, counteract or even annihilate this forcing, depending on the effected metabolite changes .
In vector notation, (1.15) becomes the flux balance
(1.19) 
Our qualitative answer to the central sensitivity question (1.13) will distinguish those fluxes and metabolites with nonzero response to the external rate change , from those with zero response, i.e. without any response at all. Let denote any flux or metabolite. We recall from (1.1) that influences , in symbols , if the component of the response vector is nonzero.
The five theoretical results of theorems 2.1–2.5 below will investigate the influence relation (1.1) in detail. Our first three results are reminiscent of [FM15], but now hold in the general multimolecular setting. Our results are based on the notion of a child selection
(1.20) 
which we define as follows. We require injectivity of , and we require each metabolite to be an input of the selected reaction . In symbols:
(1.21) 
for all . Occasionally we call a mother, or parent, of the child if , i.e. . Note that the same child may have two or more “mothers” , i.e. input metabolites . A single child , for example, possesses at least one mother which has no other children, besides itself. Then any child selection must select that single child of the mother metabolite .
Theorem 2.1 addresses the standard requirement of a nonzero Jacobian
(1.22) 
in the standard implicit function theorem; see (1.15). We recall here that denotes the stoichiometric matrix and denotes the rate matrix. We call the reaction network (1.3) regular at the steady state if (1.22) holds there. Theorem 2.1 asserts that the nondegeneracy assumption (1.22) holds, algebraically, if and only if there exists a child selection such that the columns of the stoichiometric matrix form an invertible block. This eliminates all consideration of the particular reaction rates . Indeed the condition on the child selection only involves conditions on the arrows in the reaction network (1.3), as selected by , and on the stoichiometric coefficients (1.6) associated to these arrows.
Well, can that be true? Nowhere did we even bother to exclude the disastrous case of all constant reaction rates, which leaves undetermined. We consider as a polynomial expression in the nontrivial partial derivatives of the reaction rate functions , evaluated at the steady state . We call an expression nonzero, algebraically, if it is nonzero as a polynomial, or rational function, in the real variables , taken over all metabolite inputs of the reactions . It is in this precise sense, how influence (1.1) and nondegeneracy (1.22) are considered to hold. In particular the set of where our assertions may fail to hold is an algebraic variety of codimension at least one in the space of all real with . In this precise sense our results hold true for generic rate functions . Thus theorem 2.1 clarifies the relation between the nondegeneracy assumption (1.22) and child selections in the network. See section 3 for a simple explicit example.
As a caveat we have to issue a warning concerning mass action kinetics, defined by
(1.23) 
Here are nonnegative integers, and only a single parameter is available for each reaction. The partial derivatives
(1.24) 
are closely related to the steady state rates themselves – too closely, in fact, to be considered algebraically independent. Therefore our theory does not apply to pure mass action kinetics. Slightly richer classes like MichaelisMenten, LangmuirHinshelwood, or other kinetics where each participating species , enters with an individual kinetic coefficient, fall well within the scope of our algebraic results.
Theorem 2.2 clarifies flux influence , in the same algebraic spirit. Algebraically nonzero selfinfluence occurs, if and only if there exists a child selection such that and the columns of the stoichiometric matrix form an invertible block. Similarly, algebraically nonzero influence occurs, if and only if there exists a child selection such that and the columns of the stoichiometric matrix form an invertible block. Note how the influencing column has been swapped in, to replace the influenced column of the stoichiometric matrix .
Theorem 2.3 clarifies metabolite influence , again in the same algebraic spirit. Algebraically nonzero metabolite influence occurs, if and only if there exists a child selection , on the remaining metabolites, such that and the swapped columns of the stoichiometric matrix form an invertible block. Here the influenced metabolite has been swapped out of consideration, and the range of the remaining partial child selection has been augmented by the influencing column of the stoichiometric matrix .
Transitivity theorem 2.4(ii) considers any 2step chain of influences , with either a reaction or a metabolite as the terminating target. As a consequence, direct influence is established. Compare (1.1) and (1.2). Indeed, only this transitivity of influence justifies the notion of a hierarchy of influence. For example, only by transitivity of flux influences can we assert that any finite chain of successive flux influences amounts to a direct flux influence all along the chain. Indeed a rate perturbation of the first element , in the influence chain, changes the reaction fluxes all the way, down to the very last terminating target of the whole chain.
In particular transitivity theorem 2.4 allows us to summarize all algebraically nonzero flux responses as a directed acyclic influence graph . The vertices of are the lumped subsets of reactions which mutually influence each other; see theorem 2.5. The responses of metabolites appear as annotations of the reaction vertices in the influence graph. See also figs. 3.1 and 8.2 below for examples.
The remaining sections are organized as follows. Section 3 illustrates the main results of section 2 by two examples. First we treat the simplest case of a minimal number of reactions, given that the stoichiometric matrix possesses full rank . Any flux influences turn out to be absent in that case; see (3.8). Moreover we give a detailed mathematical account of the simple bimolecular theoretical example from [MF15] which we had not been able to treat in [FM15], due to the monomolecular restriction.
Algebraic nondegeneracy of networks, the flux influence relation , and the metabolite influence relation as stated in theorems 2.1 – 2.3 involve some analysis of child selections . The theorems are proved in section 4. The augmenticity section 5 addresses the question of the fate of influence relations, and influence regions, under modifications of the network. Specifically we augment an existing network by additional reactions among the existing metabolites and/or the addition of new metabolites. See in particular theorem 5.1. Similiarly to the first examples in section 3, this illustrates theorems 2.1 – 2.3 in a general context.
In section 6 we prove transitivity theorem 2.4, which is central to all claims on hierarchy of influences. Although transitivity involves nondegeneracy , as a prerequisite, our proof is not based on child selections directly.
Symbolic packages would typically involve terms of a complexity which grows exponentially with network size. Our practical computational approach to influence graphs, as outlined in section 7, is based on integer arithmetic modulo large primes , instead. The reaction coefficients are chosen randomly .
In section 8 our approach is brought to bear on the original Ishii TCAC example [Ish+07] and the Nakahigashi augmentation [Nak+09], in several variants. In the spirit of augmenticity section 5, it turns out that the choice of exit reactions deserves particular attention.
We conclude with a discussion, in section 9. We briefly address the proper adaptation to large perturbations of reaction rates, as required by knockout experiments and in control settings. Existence and multiplicity of steady states is a related issue. Although this case did not appear in our present examples, we also comment on the appearance of stoichiometric subspaces. Revisiting augmenticity, as discussed in section LABEL:sec:B4:aug, reveals a cautioning menetekel which indicates how monomolecular exit reactions for all metabolites may destroy all hierarchic influence structure. We also point at beautiful recent progress concerning monomolecular networks. We then turn to the elegant upper estimates [Oka+15, OM16] by Okada on the influence regions within certain subnetworks, in a multimolecular setting. A glimpse at the beautiful matroid results by Murota, as summarized in [Mur09], concludes the paper: in pioneering work, Murota has established response patterns for layered matrices more than three decades ago.
Acknowledgments. We are greatly indebted to Hiroshi Matano and Nicola Vassena for their encouraging and helpful comments on mono and multimolecular reactions. To Hiroshi Matano we owe the emphasis on the CauchyBinet formula in section 4. Marty Feinberg has patiently explained his many groundbreaking results on existence, uniqueness, and multiplicity of steady states. Takashi Okada has generously explained and shared his elegant upper estimates on influence regions, as discussed in section 9. Anna Karnauhova has drawn the TCA cycles of section 8 for us, with artistry and taste. Ulrike Geiger has tirelessly typeset and corrected several versions of the first manuscript. This work has been supported by the Deutsche Forschungsgemeinschaft, SFB 910 “Control of SelfOrganizing Nonlinear Systems”.
2 Main results
In this section we present our main results, theorems 2.1–2.5. For background notation and an outline see section 1. Throughout this paper, and in all theorems, we assume surjectivity (1.11) of the stoichiometric matrix defined in (1.6). After our analysis of this condition in theorem 2.1 we also assume the reaction network (1.3) to be regular, i.e. the nondegeneracy condition
(2.1) 
holds, algebraically, for the rate Jacobian which is the product of the stoichiometric matix with the rate matrix ; see also (1.5) and (1.22). We repeat and emphasize that, here and below, any nonzero quantities in assumptions or conclusions are understood in the algebraic, polynomial sense as explained in the introduction, section 1.
Recall that we have required full rank of in (1.11), i.e. surjectivity . Consider any subset of . We say that selects an basis of the stoichiometric matrix : , if the columns of are a basis of range . This means and for the square minor of defined by the columns of . In other words,
(2.2) 
According to our notation convention (1.4) this means that : does not possess any nontrivial kernel vector supported only in .
Theorem 2.1.
See (1.20), (1.21) for the definition of child selections . Henceforth, and throughout the rest of the paper, we will assume , i.e. the existence of such a kernelfree child selection .
The meaning of the child selection will become more evident in the proof; see section 4. We will in fact show that can be written as a polynomial
(2.4) 
Here the sum runs over all child selections , and abbreviates the monomial
(2.5) 
Up to sign, the coefficient abbreviates the determinant
(2.6) 
where is the minor of the stoichiometric columns . Of course condition (2.3) is equivalent to a (truly) nonzero coefficient of the algebraic monomial .
Let us briefly illustrate the role of single children in our nondegeneracy setting . We call a reaction a single child, if there exists an input metabolite such that implies . In other words, the single child is the only child reaction of some mother metabolite of . In such a case
(2.7) 
holds for any child selection .
Let the reaction be a single child. Our standing assumption and theorem 2.1 then imply the existence of a child selection . By definition, the child selection must be injective. Hence, by (2.7), there cannot exist any other singlechild parent of , i.e. , with the same reaction as its own single child. Therefore the mother , of which is a single child, is determined uniquely by the single child in our setting of .
The next theorem characterizes flux influences in terms of swapped child selections. We recall definition (1.18) of the flux response , and how means , algebraically.
Theorem 2.2.
Assume holds, algebraically, as asserted by theorem 2.1. Then flux influence is characterized as follows.

Selfinfluence occurs, if and only if there exists a child selection : such that selects an basis, i.e.
(2.8) 
Influence occurs, if and only if there exists a child selection : such that and the swapped set selects an basis, i.e.
(2.9)
Both cases may occur for the same perturbation .
To illustrate theorem 2.2 we observe that single children have no flux influence. Indeed let be the unique singlechild mother of the single child . Then (2.7) implies . This contradicts both (2.8) and (2.9), and hence prohibits flux influence of the single child on any , including itself. We will see many examples of this simple principle later.
To characterize metabolite influences we have to define partial child selections : . We require injectivity and the input property as in (1.21), verbatim but just for all metabolites , of course.
Theorem 2.3.
Assume holds, algebraically, as asserted by theorem 2.1. Then metabolite influence occurs, if and only if there exists a partial child selection : , such that the augmented reaction set selects an basis, i.e.
(2.10) 
Again the above single child case with unique singlechild mother of the single child is instructive. We claim that single children only influence their own mother. Indeed consider , first. We have assumed . By theorem 2.1 this provides a kernelfree child selection : , see (1.21). Moreover , by (2.7). Define the partial child selection as the restriction of to . Then (2.10) follows from (2.3). Hence influences its unique singlechild mother , by theorem 2.3. Next consider . Then implies , for any partial child selection on , again by (2.7). This prevents any metabolite influence of the single child of , other than .
In the introduction we have emphasized how transitivity (1.2) of flux influence (1.1) is at the heart of any notion like “hierarchy of influence” in reaction networks. In the monomolecular case of [FM15] we were able to show how transitivity of flux influence followed from a study of paths in the directed reaction graph with vertex set . See section 9 for further discussion. We now formulate our transitivity theorem in a multimolecular setting. Instead of any reaction digraph, our proof in section 6 will involve direct differentiation with respect to an intermediate variable .
Theorem 2.4.
Assume holds, algebraically, as asserted in theorem 2.1. Then the following two transitivity properties hold.

Assume that influences an input metabolite of reaction , and influences . In symbols,
(2.11) Then , that is, influences .

Assume that influences reaction , and influences . In symbols,
(2.12) Then , that is, influences .
Consider transitivity theorem 2.4(ii) for the special case of just reactions and . Then the theorem establishes transitivity of flux influence, as claimed in (1.2) above, for the general multimolecular case.
Although we have admitted metabolites in our transitivity theorem, we have not even defined yet what an influence or of a metabolite on a metabolite or a reaction is supposed to mean. Loosely speaking, such influences describe an algebraically nonzero sensitivity response of the steady state under the addition of an artificial constant external feed of metabolite . For a precise definition see (3.4), (3.5) below, and the discussion there.
Our next goal is a description of all flux influence sets
(2.13) 
and of all metabolite influence sets
(2.14) 
for rate perturbations of any single reaction . Of course theorems 2.2 and 2.3 characterize these sets. For single children we have seen and .
Transitivity theorem 2.4 suggests a hierarchy of influence for the influence sets. To describe this hierarchy, we first construct the directed pure flux influence graph as follows. We start from an equivalence relation on , defined by artificial reflexivity
(2.15) 
and by mutual influence between different reactions :
(2.16) 
Note that reflexivity (2.15) may, or may not, be founded on actual selfinfluence . Our definition generously glosses over this delicate point.
The equivalence classes of are the vertices of the pure flux influence graph . In general, we denote an equivalence class by brackets
(2.17) 
and call it a (mutual) flux influence class. Note how any reaction actually influences itself, by transitivity, whenever the class contains at least two elements. The subtle case of single element equivalence classes, however, comes in two flavors. We write such single element classes of as
(2.18) 
We simply omit the brackets in case does not influence itself. This is the case where reflexivity was generously decreed by mathematical fiat, in (2.15), in spite of lacking actual selfinfluence.
Transitivity defines a partial order on the equivalence classes, by actual directed influence. We emphasize this important point. Abstractly, any block triangularization of a matrix can be stylized into a purely formal “partial order” of the blocks. Only in the presence of a transitivity result, however, the formal partial order becomes one of actual influence. Mathematically, this would correspond to the additional assertion that the blocks in the triangularization are fully occupied by nonzero elements.
The above partial order of actual flux influence can be expressed, equivalently, by a minimal finite directed graph . Directed paths run in the direction of, and imply, flux influence. Since the vertices are equivalence classes, does not possess any directed cycle. More precisely, a directed path from one class vertex to another class implies that each single in influences each in – but there does not exist any single influence in the opposite direction. Also, influences all in its own class – except possibly itself; see (2.18). Therefore the pure flux influence graph , in the above notation, describes the flux influence set of a rate perturbation at any given reaction .
To include the metabolite responses to , i.e. the metabolite influence sets , we define the (full) flux influence graph , by a simple annotation at the vertices of . We keep all directed edges, as defined between the equivalence classes of the pure flux influence digraph .
To determine the metabolite influence sets , we apply transitivity theorem 2.4(ii) with reactions , and observe
(2.19) 
Together with flux transitivity (1.2) this implies
(2.20)  
Now fix any mutual flux influence class := or := , i.e. any vertex of the pure flux influence graph . By (2.20), all share the same metabolic influence set
(2.21) 
We define the, possibly empty, indirect metabolite influence set as
(2.22)  
In other words, indirect influence of on requires flux influence on an intermediary agent in another flux influence class to exert its influence on via . By transitivity (2.19) that intermediary mediates the (indirect) influence This influence does not depend on the choice of the representatives , . Conversely, the intermediary agent of does not influence .
Analogously we define the, possibly empty, complementary set of direct metabolite influence
(2.23) 
to consist of all those metabolites which are influenced by , but cannot ever be influenced by any intermediary agent . In particular we obtain the decompositions
(2.24)  
(2.25) 
for . By definition, the union in (2.24) is disjoint. We may replace the intermediary agents in (2.25) by a single representative in each of their mutual flux influence classes.
Note that the classes of mutual flux influence form a partition of . The same metabolite , in contrast, may appear in several sets of direct metabolite influence. This is the case, if and only if there exist two distinct reactions such that neither influences the other, but each influences . In particular, the union in (2.25) need not be disjoint.
We define the (full) flux influence graph as follows. The annotated vertices of are given by the pairs
(2.26) 
of mutual flux influence classes , annotated by their direct metabolite influence sets ; see (2.17), (2.18) and (2.23). Edges are directed and coincide with the directed edges of the pure flux influence graph on the vertices . Summarizing our above construction and discussion, we have proved the following theorem on the transitive hierarchic structure of flux influence.
Theorem 2.5.
Assume holds, algebraically, as asserted by theorem 2.1. Define the full flux influence graph as above, and let be the flux influence class of ; see (2.17), (2.18).
Then a rate perturbation of reaction influences the flux of , i.e. , if and only if, either , or else there exists a directed path in from the vertex to the vertex of the flux influence class of . Equivalently, the same directed path runs from to in the pure flux influence graph . Selfinfluence holds unless with removed brackets; see (2.18).
A perturbation influences the metabolite ; i.e. , if and only if,

either, is under the direct influence of the class of ; see (2.23);
In the flux influence graph this means that we pass from , either, to directly in the same vertex , or else, to in the annotation of another vertex , via some directed path in .
We recall and emphasize that these nonzero influences are all generated, in unison and simultaneously, by flux perturbations at any single one reaction of the same class .
3 Two theoretical examples
By our standing assumption (1.11), the stoichiometric matrix : has full rank . Therefore : the number of metabolites does not strictly exceed the number of reactions. Our first example addresses the simplest case where . As announced in the introduction, our second example studies the (full) flux influence graph for a simple hypothetical reaction network which is monomolecular, except for a single bimolecular reaction; here . That second example was first presented in [MF15] and had to be treated in an adhoc and casebycase fashion, because the present multimolecular mathematical framework was not available at the time. Both examples illustrate the workings of theorem 2.5 about direct and indirect flux influence.
Here, and for all our proofs in sections 4–6 below, it will be convenient to rewrite the system (1.18), (1.19) for the response vector := in block matrix form as
(3.1) 
As always, we have padded the unit vector , for , with zeros in the components. We have used the block matrix
(3.2) 
in (3.1). Again denotes the stoichiometric matrix, and is the rate matrix. We recall that, by definition, we have
(3.3) 
for any .
We digress briefly to consider metabolite “perturbations” , as well. Define as the (flux, metabolite)response to an external perturbation of the metabolite . In other words, we define
(3.4) 
From an applied point of view this means that we study the (infinitesimal) steady state response of , and the associated flux changes , to external feeds of ,
(3.5) 
as the (infinitesimal) rate of that feed varies.
After this metabolic digression we now present the case as our first example. The matrices and are square, and
(3.6) 
Note , by our full rank assumption (1.11).
The steady state equation (1.12) becomes , for . If for , componentwise, this precludes the existence of positive steady states . However, we did not impose any such positivity restrictions on . Even in the monomolecular case, for example, a feed reaction in (1.3) can be lumped with a subsequent forward reaction into a single reaction term like , at the expense of violating positivity .
The child selections : of theorem 2.1 readily appear in the evaluation of . In fact, (2.4), (2.6) hold with
(3.7) 
if we view as a permutation of elements with signature . Our nondegeneracy assumption of theorem 2.2 amounts to (algebraic) invertibility of .
Theorem 2.2 informs us that there are no influences at all, since is impossible. Of course this also follows directly, because in (1.19) and imply , for any . Thus the flux influence sets of (2.13) are all empty. The pure flux influence graph consists of the isolated vertices , sadly without any edges. The full flux influence graph , then, does not possess any edges either. Therefore all metabolite influences are direct, i.e.
(3.8) 
Theorem 2.3 implies that , if and only if , for some child selection permutation : . Indeed . Hence
(3.9) 
Clearly the metabolite influence sets will not be disjoint, here, if the reaction network admits more than one child selection permutation .
We now revisit the adhoc bimolecular example of [MF15] with 15 reactions , 10 metabolites , and a single bimolecular reaction
(3.10) 
see fig. 3.1(a). In spite of this notation inherited from [MF15], we do not plan to confuse metabolites B,F,J with the matrix , the pure flux influence graph , or child selections .
The matrix of (3.1), (3.2) can easily be inverted symbolically. In our example,
(3.11) 
We have omitted the redundant input metabolite index in , for monomolecular reactions . In (3.11) we see explicitly what it means for to be nonzero, algebraically. We have to require nonzero prefactors and the linear nondegeneracy . In this explicit sense, the bimolecular chemical reaction network of fig. 3.1(a) is regular, algebraically; see also theorem 2.1.
The feed reactions are and the exit reactions are . We have omitted the formal metabolite entry in these open system reactions. The seven single children , as defined just after theorem 2.2, are
(3.12) 
By theorem 2.2 they have no flux influence, and therefore define terminal sinks of the pure flux influence graph ; see fig. 3.1(b). The remaining terminal sink is not a single child.
Let us study (3.11) in a little more detail, in terms of child selections and (2.3)–(2.7). The seven single children in (3.12) are forced to appear in the index set of any monomial in (3.11). Indeed, their unique respective singlechild mothers have no choice but to select their only child ; see (1.21) and (2.7). Since the bimolecular reaction is the single child of metabolite H, this also forces the index of to appear in the prefactor of (3.11). The only remaining choices for a child selection are
(3.13) 
The choice immediately leads to the indicator of the network cycle as a kernel vector of supported on . Indeed is a single child, and we just saw why is also forced to hold. Therefore we must choose in theorem 2.1. This explains why must also appear in the prefactor of (3.11).
We could easily proceed with such elementary arguments to fully derive (3.11) by hand. Likewise, the modified child selections of theorems 2.2 and 2.3 would determine all flux and metabolite influences of rate perturbations , explicitly. Instead, we present a symbolic version of in figure 3.2. A black area entry for in row and column of is equivalent to an algebraically nonzero component of the response, i.e. to the influence ; see (3.1), (3.3). By the construction (2.21)–(2.23), (2.26) this immediately defines the flux influence digraphs of fig. 3.1(b),(c). Theorem 2.5 explains how these graphs determine all flux influence sets and metabolite influence sets , explicitly, by following directed paths downward in the diagrams. Note the motherchild pairs which provide terminal sink vertices in , fig. 3.1(c). The bimolecular input of reaction , interestingly, cannot respond to a perturbation of