Good and bad children in metabolic networks

: Equilibrium bifurcations arise from sign changes of Jacobian determinants, as parameters are varied. Therefore we address the Jacobian determinant for metabolic networks with general reaction kinetics. Our approach is based on the concept of Child Selections: each (mother) metabolite is mapped, injectively, to one of those (child) reactions that it drives as an input. Our analysis distinguishes reaction network Jacobians with constant sign from the bifurcation case, where that sign depends on speciﬁc reaction rates. In particular, we distinguish “good” Child Selections, which do not a ﬀ ect the sign, from more interesting and mischievous “bad” children, which gang up towards sign changes, instability, and bifurcation.


Introduction
One of the major obstacles in the study of dynamical systems arising from chemical reaction networks is the lack of precise quantitative knowledge of the parameters involved. For this reason, a natural approach is to connect dynamical properties with the network structure, only. In fact, this idea has a long tradition, in chemical reactions settings. Pioneering works by Horn&Jackson [1] and Feinberg [2,3] have been directed towards sufficient network conditions to assert existence and uniqueness of a stable equilibrium. The concept of deficiency has been developed towards this purpose. Further network features such as injectivity [4][5][6][7] and concordance [8,9] emerged in this context to forbid multistationarity. In particular, saddle-node bifurcations are excluded.
In another direction, in 1981 Thomas [10] conjectured a necessary network condition for multistationarity. In particular, certain positive loops in the network were identified as necessary network motives to sustain multistationarity. The conjecture was then fully proven in mathematical settings by Soulé [11]. Refined conditions for restricted classes of networks have been further studied till present  [19] and the graphical representation is courtesy of Anna Karnauhova. The inflow feed reaction is named f 1 . Outflow exit reactions are labeled d1 − d6 and dd1 − dd9. Here, for image simplicity, a reversible arrow reaction m ←→m encodes two different opposite reactions. Metabolites PEP, PYR and CO2 have been graphically repeated, for sake of clarity of the picture.
An excellent compendium of the various approaches to multistationarity questions can be found in [15], where the authors also presented and clarified the important related work of Ivanova [16][17][18], independently developed in the '70s and unfortunately not easily accessible to english readers.
Focusing on metabolic networks, the novelty of our present contribution is in two regards. Firstly, we fully characterize on a graphical level the sign of the determinant of the Jacobian matrix. Secondly, we provide a recipe to detect sign changes of the Jacobian determinant in quite general examples. This implies the identification of parameter areas for equilibria bifurcations, leading to stability changes and multistationarity.
Metabolism is central to life. A metabolic process is a sequence of chemical reactions designed to transform nutrients into energy. A prominent example of a metabolic network is the central carbon metabolism of Escherichia coli, figure 1. The network depicts the metabolic transformation of glucose, which allows the production of energy. Furthermore, this network is widely used as a general model for cellular respiration, due to the fact that E. coli are genetically well-known and easy to treat. We will concentrate on this network in Section 6. However, this reference example helps us at first to focus on some mathematical features of metabolic networks.
In particular, we underline ( ): 1. The stoichiometric coefficients in metabolic networks are mostly {1,0}. In the E. coli network in figure 1 they are only {1,0}. 2. The number of metabolites involved in a reaction is considerably low compared to the size of the network. The E. coli network in figure 1 consists of 30 metabolites interacting through 58 reactions. The maximum number of metabolites involved in a reaction is four (e.g. reaction 14), and many reactions involve only two metabolites (e.g. reaction 10). 3. The reaction rate functions (kinetics) describe the mathematical form of the reactions. Metabolic networks, such as the one in figure 1, represent only the metabolic transformations without considering the enzymes network. That is, enzymes and other secondary chemicals do not appear explicitly in the network. The enzymes are usually taken in account by using enzymatic kinetics (e.g. Michaelis-Menten) instead of elementary kinetics (e.g. mass action).
Using the above outline as a guide, we can begin to discuss the details of our approach.
We consider general metabolic chemical reaction networks Γ with M metabolites and N reactions. For notation, we use labels A, B, C, D, ... for metabolites and 1, 2, 3, ... for reactions. We call M the set of metabolites and E the set of reactions, such that |M| = M and |E| = N. We use the small letter m ∈ M for a generic metabolite and the small letter j ∈ E for a generic reaction.
A chemical reaction j is represented as with nonnegative stoichiometric coefficients s j ,s j ∈ R. In a metabolic context, we repeat, these coefficients are mostly 0 or 1. A metabolite m is called an input or a reactant of the reaction j, if s j m 0. Respectively, m is called an output or a product of the reaction j, ifs j m 0. We say, conversely, that a reaction j is outgoing from the metabolite m if m is an input of reaction j. We say that a reaction j is an ingoing reaction of the metabolite m if m is an output of j.
Metabolic systems are intrinsically open systems, that is, they exchange chemicals with the outside environment by feed and exit reactions. Within our settings, the constant feed reactions, or inflows, are reactions with no inputs (s j = 0) and the exit reactions, or outflows, are reactions with no outputs (s j = 0).
Graphically, we represent a reaction where we have omitted stoichiometrically zero terms, as follows.
In (1.3), the arrow orientation is inherited from (1.2). The stoichiometric coefficient 2 of metabolite B is indicated as a weight in the lower tail of the directed arrow j, and stoichiometric coefficients 1 are omitted, as well as non-participating other reactants. In particular, this graphical representation considers the metabolites as vertices and the reactions as arrows of the network, and it is one natural representation widely used in chemistry, biology, and mathematics.
Explicit autocatalytic reactions j are defined as reactions for which a metabolite m is both an input and an output of the reaction. In symbols, s j m ,s j m 0, for at least one metabolite m. Throughout this paper, we exclude explicit autocatalytic reactions. In particular, self-loops are not allowed in the graphical representation of the network.
This assumption is only made for mathematical simplicity. Explicit autocatalysis can be treated via intermediaries. For example, the reaction is explicit autocatalytic. However, the same chemical process can be modeled without explicit autocatalytic reactions as with m 0 as an intermediary.
To construct the M × N stoichiometric matrix S , let us consider any reaction j. We associate to any stoichiometric coefficient s j m of an input metabolite m of the reaction j a negative stoichiometric entry of the stoichiometric matrix S, that is: Conversely, we associate to any stoichiometric coefficients j m of an output metabolite m of the reaction j a positive stoichiometric entry of S, that is: For example, a monomolecular reaction j is a reaction which possesses as input one single metabolite m 1 and as output one single metabolite m 2 , Such a reaction translates into the j th column of the stoichiometric matrix S as (1.9) Here we have indicated the rows by m 1 , ..., m M , explicitly. In particular, with this construction we model a reversible reaction simply as two distinct irreversible reactions Columns associated with feed reactions contain only positive entries and the ones associated with exit reactions contain only negative entries. All other columns contain both positive and negative entries.
In the E. coli picture of figure 1, only for image simplicity, a reversible arrow reaction m ←→m encodes two different opposite reactions. Stoichiometric matrices of metabolic networks, due to our observations ( ), are then sparse matrices with most entries being {-1, 0, 1}.
Let x m (t) be the time evolution of the concentration of the metabolite m. The isothermal dynamics of the vector x ∈ R M of the concentrations is described by the system of differential equationṡ (1.12) The M × N matrix S is the stoichiometric matrix constructed above. The N-dimensional vector r(x) represents the reaction rates as functions of x: the kinetics of the system. The feed reactions are represented by constant functions; i.e.: for a feed reaction j f . Throughout this paper, we pose the following assumptions on the reaction rates r(x): 1. We assume the reaction rates r j (x) to depend only on those concentrations x m such that the metabolite m is an input metabolite of reaction j. In particular, Moreover, we use the notation r jm for the nonzero partial derivatives, i.e., if m is an input of reaction j.

2.
We consider strictly positive monotone reaction rate functions r j (x) ∈ C 1 , for every j = 1, ..., N: and, for the nonzero partial derivatives r jm , strictly positive slopes This monotonicity restriction is indeed satisfied for most, but not all, chemical reaction schemes.
With these assumptions, all required information of the network is encoded in the stoichiometric matrix S , only. In particular, we do not specify here the mathematical form of the kinetics, remaining in wide generalities.
Great effort in the mathematical community has been spent in finding network characterizations of the existence and the uniqueness of equilibrium solutions of (1.12). For extensive reference, see the comprehensive book by Feinberg [20].
With our approach, we do not address this question at all. In fact, throughout this paper, we assume the existence of a dynamical equilibrium x * that solves (1.14) The assumption of the existence of a dynamical equilibrium is not smoothly untroubled. In particular, linear constraints have been implicitly imposed on the reaction rates r, because of (1.14). Note that these constraints do not necessarily fix the precise value of an equilibrium x * , and can be considered posed a priori, so that the existence of the equilibrium is an assumption on the reaction rates r, only.
Here, our analysis is based entirely on the derivatives r jm of the reaction rates and we do not want to be concerned by the equilibrium constraints (1.14). To avoid this, we must assume a certain independence of the derivatives r jm from the reaction rates themselves, at the equilibrium. In particular, locally at the fixed equilibrium, we require the possibility of freely choosing the value of any r jm independently from each other and from the constraints S r = 0. In this sense, the partial derivatives r jm can be considered positive free parameters. This requires a certain mathematical complexity of the reaction rates r j . In fact: too mathematically 'simple' kinetics fail to satisfy this assumption.
As an example, for polynomial mass action kinetics, the value of r j (x) and r jm (x) are related, a priori, at any value x, and for any j and m. In particular, the theory developed here does not fully apply to mass action kinetics. In contrast, Michaelis-Menten kinetics, more suited for metabolic networks, satisfies our independence assumption. We present here an exemplification of this fact. The complete mathematical argument can be found in [21], which is a contribution by Fiedler that shares this premise. Let us consider the following reaction, whose single input is a metabolite m: The product of the reaction is irrelevant for this discussion, and it is omitted. The rate of the reaction 1, according to the law of mass action, reads Here, k 1 is a positive coefficient. The derivative of r 1 (x m ) with respect to x m is then In particular, at any fixedx m , the value r 1 (x m ) uniquely determines the derivative r 1m (x m ). Indeed, For general nonlinearities, of course, we may choose r 1 (x m ) and r 1m (x m ) independently of each other.
In the special case of Michaelis-Menten kinetics, for example, the rate of reaction 1 is where both k 1 and a 1 are positive parameters, and its derivative reads (1. 20) This implies that at any fixedx m , the value r 1m (x m ) can be chosen as small as needed, independently from r 1 (x m ): In particular, for large a 1 → ∞ we get small r 1m → 0, and for small a 1 → 0 we get r 1m → r 1 (x m ) x m . Having fixed the valuex m requires then to pick k 1 := (1+a 1xm )r 1 (x m ) x m , for any chosen a 1 . Note that for the limit value a 1 = 0 we recover the mass action case.
Let us concentrate now on the relevant case in which the fixedx is an equilibrium. In the mass action case, it is not possible to separate the question of the existence of the equilibrium (related to the reaction rates r(x)) from the question of its stability (related to the derivatives r jm (x)), because the former determines the latter as shown in the analysis above. On the contrary, for more general kinetics as Michaelis-Menten, the two questions can be addressed separately and independently.
The Jacobian matrix of a dynamical system plays a central role in the stability analysis of equilibria. The sign of the eigenvalues is an indication of stability; therefore, a change in sign of the determinant hints at a change in stability and bifurcation phenomena. Let G be the Jacobian matrix of the equilibrium system 1.14, that is: Note that the entries G i j are linear forms in the variables r jm .
The leading question of this paper is the following: When is det G of fixed sign? (1.23) That is: When -for any choice of positive parameters r jm -does the determinant carry the same sign?
In this paper, we fully characterize, on a graphical level, the answer to the question (1.23). In particular, in Section 2, we introduce Child Selections and we use the Cauchy-Binet formula to expand the Jacobian determinant in a polynomial, in which each monomial summand is associated to a Child Selection. Depending on the sign of the coefficients of these monomials, each Child Selection is abstractly classified in good or bad. The coexistence of a good and a bad Child Selection characterizes the condition of indeterminate sign Jacobian. In Section 3, the main Theorem 3.1 abstractly characterizes whether a given Child Selection is good or bad and Section 4 translates this abstract condition into a graphical network condition. Consequently, Section 5 uses the developed concepts to find a bifurcation parameter responsible for a change of sign in the Jacobian determinant, with possible consequent bifurcation phenomena. Section 6 contains an example of an application for the central metabolism of E. coli. Section 7 and 8 conclude with the discussion and proofs.

Cauchy-binet analysis via child selections
For the metabolic chemical reaction system 1.12 the Jacobian matrix reads The reactivity matrix R of the partial derivatives is an N × M matrix, whose entries R jm are given by: The entry R jm is nonzero, i.e. R jm = r jm , if and only if the metabolite m is an input of the reaction j.
The algebraic structure of G is thus completely characterized by the network structure, only.
The following definition, originally introduced by Brehm and Fiedler [19], presents a central tool for this paper.
Definition 1 (Child Selections, mothers, children). A Child Selection is an injective map J : M −→ E, which associates to every metabolite m ∈ M a reaction j ∈ E such that m is an input metabolite of reaction j. We call the reaction j = J(m) a child of m, and the metabolite m = J −1 ( j) a mother of the reaction j. Remark 1. It is possible that a metabolite m is an input of j but not a mother of j, due to injectivity of Child Selections. Indeed, consider the following example: In this minimal case, metabolite B is an input of reaction 1, but never a mother: 1 is the only outgoing reaction from metabolite A. Therefore, due to injectivity, J −1 (1) = A.
The following Jacobian analysis, based on the Cauchy-Binet formula, is developed from a previous result in [19].
Proposition 2.1. Let G be a network Jacobian matrix, in the above settings. Then: where S J is the M × M matrix whose m th column is the J(m) th column of S .
Remark 2. Because of the exclusion of explicit autocatalytic reactions from our model, the stoichiometric coefficient S m j of an input metabolite m to reaction j is a priori negative. This implies that any diagonal element of S J is always negative: S J mm = S mJ(m) < 0 for any Child Selection J and any metabolite m.
Remark 3. If there are no Child Selections at all, then det(G) ≡ 0 for any choice of parameters r jm . For example, this is the case of a network where there are more metabolites than reactions. Another possible example is when two metabolites m 1 and m 2 are both inputs to one reaction j, and they are inputs to no other reaction. In such a situation, due to the injectivity requirement in Definition 1, we automatically have no Child Selections and det(G) ≡ 0 for all parameters. Note that, in metabolic networks, this is practically never the case due to the omnipresence of outflow reactions.
Via Proposition 2.1, the Jacobian determinant of a metabolic network is a homogeneous multilinear polynomial in the variables r jm . The possible sign of the polynomial depends on the signs of the monomial coefficients det S J . Hence, it is natural to state the following classification of Child Selections, according to the sign of the determinant of the reshuffled minor S J .
Definition 2 (Child Selection behavior). Let J be a Child Selection. We say that J is bad, or J ill-behaves, if sign(det S J ) = (−1) M−1 . We say that J is good, or J well-behaves, if sign(det S J ) = (−1) M . If det S J = 0, we say that J zero-behaves.
The choice of the terminology is consistent. In fact, in a metabolic network context, important classes of Child Selections well-behave, for example all acyclic Child Selections (see Section 4). Moreover, a 'stable' Child Selection J, for which all the eigenvalues of S J have negative real part, well-behaves. In this sense, a ill-behaving Child Selection is an indication of possible instability.

Example 1: bad child selection
This network represents a bad Child Selection: The stoichiometric matrix associated is

Example 2: good Child Selection
This network represents a good Child Selection: The stoichiometric matrix associated is At this point, the reader may wonder whether Definition 2 is well-posed in the first place; i.e., whether the behavior of a Child Selection depends on the specific labeling of the network. Section 4, Remark 8, clarifies this point, assuring the well-posedness of the definition.

Main result
In metabolic networks, stoichiometric coefficients are mostly 0 and 1. For this reason, we continue here assuming that S has entries S m j ∈ {−1, 0, 1}. By 2, Remark 2, then, the diagonal entries of the reshuffled minor S J are S J mm ≡ −1, for any m. We refer to the Phd thesis [22] for a more general version of Theorem 3.1, accounting for real stoichiometric coefficients S m j ∈ R.
Here, via a structural analysis of det S J , we characterize whether the given Child Selection J well-behaves or ill-behaves. Note, however, that the importance of the result is mainly revealed in its interpretation, see Section 4.
The Leibniz expansion formula for the determinant, applied to S J , reads Again, π indicates permutations of M elements and sgn(π) is the signature of π. Let denote the summand associated to the permutation π in the Leibniz expansion. For example, denoting as Id the identity permutation, Let π Id be a permutation such that E(π) 0. Combinatorially, the permutation π can be expressed as the product of ϑ disjoint permutation cycles c i of length l i > 1, Again, ϑ is the number of cycles in the permutation expansion. If π consists of a single cycle c, i.e. for ϑ = 1, we call the good (resp. bad)-completion a good(resp. bad)-cycle.
We clarify in the next Section 4 what a completion completes, as it requires some further arguments. Firstly, given the above definition, we state the main result of this section. For a Child Selection J, Theorem 3.1 characterizes the sign of det S J in terms of permutation cycles of the Leibniz expansion (3.1). The following Section 4 relates these permutation cycles to certain cycles in the network.

Graphical interpretation of the result
The Metabolite-Reaction graph (MR-graph) is an undirected bipartite graph with two sets of vertices given by the metabolites m 1 , ..., m M and the reactions j 1 , ..., j E , respectively. For a metabolite m participating in a reaction j, edges e = (m, j) are adjacent to a metabolite vertex m and a reaction vertex j. With this construction, edges in the MR-graph are in a one-to-one relation with the nonzero entries of stoichiometric matrix S . In particular, with S m j < 0 of in mind, we call positive the edges e = (m, j) where m is output to j. See figure 2 for a comparison between different kinds of representation graphs for the same network. Under the name Species-Reaction graph (S R-graph), this was considered by [5] and others.
We proceed with two definitions and a proposition.  Remark 6. Equivalently, because of Remark 5 above, a completion cycle is a cycle in the MR-graph of length 2l, such that J-selected edges alternate with non J-selected edges. It is now clear what the word 'completion' refers to. In fact, any completion cycle is constructed by completing J-selected elements to a cycle of length 2l in the MR-graph. In this sense, a goodcompletion π = ϑ i=1 c i can be seen, in the MR-graph, as a collection of ϑ non-intersecting completion cycles c i , such that the number of good-cycles has the same parity of ϑ. Respectively, a bad-completion is a collection of ϑ non-intersecting completion cycles, such that the number of good-cycles has opposite parity of ϑ.
Examples BIOLOGICAL Remark 7. From Definition 3 and Proposition 4.1 it follows that, for a completion cycle in the MRgraph, the parity of the negative edges not J-selected characterizes the cycle being good or bad. In fact: an odd number of negative edges not J-selected corresponds to a good completion cycle, an even number corresponds to a bad one.

MR-GRAPH Biological
Remark 8. Obviously, being a network structure, the definition of a completion cycle does not depend on the specific labeling of the network. Proposition 4.1, together with Theorem 3.1, guarantees in particular that also Definition 2 does not depend on any labeling and, thus, is well-posed.
Finally, we list some consequences of Theorem 3.1, useful for applications. well-behaves.
In particular, the bad Example 1 falls into the casuistry of the point (4) of Corollary 4.2, and the good Example 2 into the casuistry of point (2) of Corollary 4.2. Any Child Selection with a single monomolecular cycle is an example for (3). Such list can be of course enlarged for a given network of application, to easily detect good and bad Child Selections for bifurcation analysis, as discussed in the continuation of this paper.

Hunting saddle-node bifurcations
Here, we give a simple network condition under which there is the possibility, for certain parameters, of a saddle-node bifurcation of equilibria. We also identify bifurcation parameters responsible for the change of sign of the determinant and consequent change of stability of any equilibrium.
where the omitted terms can be chosen arbitrarily small, and a, b are coefficients of the same sign. In particular, the parameter may serve as a bifurcation parameter for the bifurcation of nontrivial equilibrium solutions of the system (1.12).
Theorem 5.1 does not require fixing an equilibrium, and it holds in more general settings. Nevertheless, for simplicity, we are thinking here of an equilibrium situation. The parameter ξ = ar J 1 (m b )m b − br J 2 (m b )m b is 'localized' in a single metabolite m b . In fact, the change of stability is driven by the difference between the derivatives with respect to the same metabolite m b of the reaction rates of two child reactions of m b itself. This suggests a simple biological scheme for controlling the unstable dimension of the equilibrium.

A case study: The central metabolism of E. coli
The central metabolism of E. coli (figure 1) consists of different and interconnected parts. In particular, the upper part comprises the so-called Pentose phosphate pathway and the Glycolysis. The bottom 'cyclic' part includes the Tricarboxylic acid cycle and the Glyoxylate cycle. We skip here more detailed biological explanation [23,24].
In this section, we analyze the network of the central metabolism of E. coli of figure 1. This network representation is primarily based on the original model proposed by Ishii et al. in [25] with the modifications suggested by Nakahigashi et al. in [26]. Note that, in biology papers, 'obvious' outflow exit reactions, such as d1 -d6 shown here, are frequently omitted. For our mathematical analysis, however, we are bound to include them. These reactions are the only outgoing reactions of their input metabolites. Their omission would result in an infinite production of their input metabolites and in a mathematical degeneracy of the network.
The network possesses 30 metabolites and 58 reactions. The number of Child Selections is on the order of 10 7 . Nevertheless, we can provide interesting biological insights without computing such a huge amount of Child Selections. In the same spirit as Section 5, and along its lines, we find two Child Selections J 1 and J 2 with opposite behavior, such that J 1 (m b ) J 2 (m b ) for one single metabolite m b and J 1 (m) = J 2 (m) for all other metabolites m m b . This situation, via Theorem 5.1, provides a bifurcation parameter responsible for a change of sign in the Jacobian determinant and possible consequent saddle-node bifurcations of equilibria.
To find the two Child Selections J 1 and J 2 as above, we start by imposing certain child reactions j to certain mother metabolites m. We do this arbitrarily, and only for sake of exemplification. Many other choices and analogous constructions are, of course, possible. Let us consider only Child Selections associating the metabolites PEP, PYR, and CO2 to their respective exit reactions, that is: For any Child Selection satisfying these constraints, we can consider the upper part (Pent. Phosph. Pathway -Glycolysis) and the bottom part (Tricarboxylic acid cycle -Glyoxylate cycle) as separate and independent. In fact, any Child Selection J satisfying the above constraints 1-3, identifies reshuffled minors S J , which are block diagonal. This shows that certain qualitative arguments on the dynamics of the central metabolism may be inferred, separately, from the biological components of the network. For example, for a block diagonal Jacobian matrix in our settings, indeterminate sign determinant of one block trivially implies indeterminate sign determinant for the entire matrix. In particular, we may concentrate on the bottom part of the network, only assuming that J 1 = J 2 in the upper part. By looking at the MR-graph representation, we can easily conclude that J 1 well-behaves and J 2 ill-behaves. Note indeed that J 1 does not contain any completion cycle, and therefore well-behaves. In fact, this Child Selection contains only one network cycle c = MAL − 23 − OAA − 17 − AcCoa − 27 − MAL, which is not a completion cycle as the edge AcCoa − 27 is not J 1 -selected. On the other hand, the completion cycles structure of J 2 is identical to Example 1, which had provided a simple and recognizable pattern of an bad Child Selection. In fact, this Child Selection possesses only two bad completion cycle c 1 and c 2 : Since it possesses only two intersecting bad completion cycles, the Child Selection J 26 ill-behaves. In particular, in accordance to Theorem 5.1, the parameter controls a change of sign of the Jacobian determinant of the entire system, for a certain region of parameters.
Remark 9. The choice of reaction 19 and 26 basically highlights the difference between the Tricarboxylic acid cycle (reaction 19) and the Glyoxylate cycle (reaction 26). Our analysis suggests how the control of certain dynamical properties of the metabolism of a cell is related to its network structure.

Discussion
In this paper, we have presented a new approach to address questions related to the Jacobian determinant for dynamical systems arising from metabolic networks. The idea of our approach relies on considering the partial derivatives r jm (x * ) independent, at a fixed equilibrium x * , from the value of the equilibrium x * itself and from the value r j (x * ) attained by the reaction rates at the equilibrium. Importantly, this allows us to separate the question of the existence of the equilibrium, which concern the reaction rates r j , from the questions related to the stability of the equilibrium, which concern the derivatives r jm . In particular, we interpret the Jacobian determinant as a homogeneous multilinear polynomial in the variables r jm . The coefficients of each monomial are determinants of reshuffled stoichiometric minors S J , and we have given network conditions to establish the sign of det S J .
Mass action kinetics forbids considering the partial derivatives r jm as free parameters, at an equilibrium. Still, the present structural analysis of the Jacobian determinant holds identically also for the mass action case. In particular, a Jacobian determinant of fixed sign is still a sufficient condition to exclude saddle-node bifurcations also under mass action kinetics, and the copresence of both good and bad Child Selections is therefore still a necessary condition. However, in the case where the sign of the Jacobian depends on the parameters r jm , we are not able to directly conclude on equilibria bifurcations, for mass action systems, given that the existence of the equilibrium depends as well on related parameters. Hence, the consistency between the parameter constraints coming from the existence of the equilibrium and the constraints coming from the Jacobian determinant must be additionally checked. This introduces a further difficulty in bifurcation analysis under mass action kinetics.
Most of the arguments can be lifted to any dynamical system of the forṁ where A is any real matrix (see [22]). However, the efficacy of this approach is most exploited for the case where A = S is a sparse matrix with low integer entries, as it is the case of stoichiometric matrices of metabolic networks. In this context, our approach provides bifurcation patterns and biological intuition.
Theorem 5.1 prescribes how to find parameters for a change of stability of equilibria. In particular, we shall find one good Child Selection J 1 and one bad Child Selection J 2 such that J 1 (m b ) J 2 (m b ) for a single metabolite m b , and J 1 (m) = J 2 (m) for all other metabolites m m b .
We point out that finding good Child Selections is reasonably easy, in metabolic networks. As an example, the sparsity of the stoichiometric matrix S as well as the overall presence of outflows make the task of finding acyclic Child Selections simple and, by Corollary 4.2, acyclic Child Selections are always good. An analogous argument supporting the predominance of good Child Selections follows by the fact that many reactions in a metabolic networks are monomolecular, always via 4.2.
Morally, then, the fewer bad Child Selections are the important ones to be found in the quest for parameter areas where bifurcations and changes of stability take place. This last observation underlines the importance to find simple and recognizable patterns of bad Child Selections. Example 1 of this paper is possibly one of those examples. Designed as the simplest example of a bad Child Selection on three metabolites, it turned out to possess a cycle structure that was easy to recognize in the real example of the central carbon metabolism of Escherichia coli.
Previous studies addressed equilibria bifurcations in the central carbon metabolism of E. coli. In [27] the authors have numerically studied bifurcations in the restricted E. coli model proposed by Chassagnole and coauthors [28]. They have focused on saddle-node bifurcations and Hopf bifurcations. In a similar flavor to the present paper, their bifurcation parameters are associated to individual reactions j. In particular, the rate r j is written as r j = k j g(x), where k j is a scalar and g(x) the given kinetics. They use the scalars k j as bifurcation parameters. However, changing this parameter also influences the reference equilibrium x * . For mathematical simplicity, in contrast, our approach fixes x * and changes the partial derivative r jm , only. With these adaptations in mind, our analytical results match their numerical explorations well. For example, consider the following subnetwork from their restricted model: Here reactions 1, 2, 4, 5 are outflow reactions, and PYR is not an input to any other reaction. Their bifurcation diagram for the two parameters (k 1 , k 2 ) is similar, qualitatively, as for (k 1 , k 3 ). Qualitatively similar diagrams also appear when considering the couples (k 1 , k 4 ) and (k 1 , k 5 ). This is consistent with our Child Selections point of view. Indeed: assume that J is any good Child Selection with J(PEP) = 2. Then, it is good also the Child SelectionJ such thatJ(PEP) = 3 andJ(m) = J(m) for any other metabolite m, because the cycle structure remains identical: the subnetwork (7.1) is acyclic. The same argument holds also for the other two 'sisters', reactions 4 and 5. Our analysis consistently suggests that the two parameters k 2 and k 3 (k 4 and k 5 , respectively) are bifurcation parameters that can be used in an analogous way, qualitatively. The absence of equilibria bifurcations in the E. coli network, both for mass action and Michaelis-Menten kinetics, has been claimed in [29]. Curiously, that study misses the sign changes of the Jacobian determinant, which we encounter for the Michaelis-Menten case.
Mathematically, global Hopf bifurcation of time periodic oscillations of the E. coli network, and specifically for the citric acid cycle part, has also been established along the lines of our present analysis [21].
There are several applications of the theory presented in this paper. Firstly, our theory suggests promising bifurcations parameters based on network information, only. In large systems, finding bifurcation parameters for numerical simulations is not a simple task. Our recipe of finding two Child Selections, one good and one bad, serves as an aid for this purpose. However, for realistic explanation of metabolic phenomena the mathematical approach should best be combined with deep biological insight. Further, the tools developed here also apply to sensitivity analysis of equilibria. See the thesis [22] where -in the same setting -responses to perturbations of reaction rates and metabolites concentration have been studied. The responses may have a definite sign for all parameters, or that sign may depend on the parameters choice, pretty much in the same spirit as presented in the present paper. The case in which the sign of the responses depends on the parameters suggests the possibility of controlling the sign, as the responses may be positive or negative depending on the parameters. For example, note that stable equilibria x * require the sign of the Jacobian to be (−1) M for M metabolites. Such control may therefore indicate external switches between different stable metabolic pathways, or may even induce stem cell differentiation.

Proofs
This Section is devoted to the proofs of our results. We start with Proposition 2.1.
Proof of Proposition 2.1. We apply the Cauchy-Binet formula on G = S R to obtain: Here π indicates a permutation of M elements and sgn(π) is the signature (or parity) of π. Note that m∈M r π(m)m 0 if and only if there is an associated Child Selection J such that r J(m)m = r π(m)m , for every m. In particular, the sum runs non trivially only for the selected minors S E such that the set E is the image of M through a Child Selection J. Now,  We proceed with the proof of main technical result, Theorem 3.1.

(8.4)
Let h be the number of elements m such that π(m) m. That is, h is the number of elements of π which are not fixed points of the permutation, but belong to a permutation cycle.
The steps above are made noting that (S J mm ) 2 ≡ 1, for any m and that, for a cycle c of length , sgn(c)(−1) = −1. We conclude the proof by observing that if π is a good-completion (bad-completion, respectively). This yields to the identity which proves the Theorem.
The interpretation of the above result has been discussed in the Proposition 4.1 and consequent Corollary 4.2. Therefore, for a non-zero Child Selection J possessing only monomolecular reactions and a single bimolecular reactionj of the form (4.4), either J does not contain any completion cycle and consequently well-behaves, due to point (1) of this corollary, or J contains completion cycles in the MR-graph which includej as a reaction vertex. In particular all completion cycles intersects. Without loss of generalities, assume J −1 (j) = A. We only have two possibilities for a completion cycle: a bad completion cycle containing the two adjacent edges A −j − C or a good completion cycle containing the two adjacent edges A −j − B. By assumption on the network, the completion cycles run through monomolecular reactions only, except fromj. Because they do intersect, we have at most one bad completion, G ≤ 1, and one good completion, B ≤ 1. Hence, 1 − G + B ≥ 0 and by Theorem 3.1 the Child Selection J does not ill-behave. Since, by assumption, J does not zero-behave, it must well-behave.
The last remaining proof is Theorem 5.1. We introduce some concepts, first. The set of Child Selections {J} carries a natural integer-valued distance d. If we further assume that J 1 and J 2 are such that one well-behaves and the other ill-behaves we have: with a and b constants of the same sign.
By the mere fact that d is an integer-valued distance, any other Child Selection satisfies d(J k , J 1 ), d(J k , J 2 ) ≥ 1, for any k 1, 2. (8.11) In particular, we have the following Lemma: We can consider, then, an -small choice of reaction rate parameters such that r J k (m k )m k < . The parameter ξ = ar J 1 (m b )m b − br J 2 (m b )m b becomes then a bifurcation parameter for the sign of the Jacobian determinant.