Good and bad children in metabolic networks

1.   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 days. See ^[12,13,14], among many others.

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 ($\star$):

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 $\mathbf{\Gamma}$ 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 $\textbf{M}$ the set of metabolites and $\textbf{E}$ the set of reactions, such that $|\mathbf{M}| = M$ and $|\mathbf{E}| = N$. We use the small letter $m \in \textbf{M}$ for a generic metabolite and the small letter $j\in \mathbf{E}$ for a generic reaction.

A chemical reaction $j$ is represented as

with nonnegative stoichiometric coefficients $s^{j}, \bar{s}^j \in \mathbb{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 \neq 0$. Respectively, $m$ is called an output or a product of the reaction $j$, if $\bar{s}^{j}_m\neq 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 ($\bar{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, \bar{s}^j_m \neq 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\times 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 coefficient $\bar{s}^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

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\longleftrightarrow \tilde{m}$ encodes two different opposite reactions.

Stoichiometric matrices of metabolic networks, due to our observations ($\star$), 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 \in \mathbb{R}^M$ of the concentrations is described by the system of differential equations

The $M \times N$ matrix $S$ is the stoichiometric matrix constructed above. The $N$-dimensional vector $\mathbf{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 $\mathbf{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) \in 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

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 $\mathbf{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 $\mathbf{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\mathbf{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 fixed $\bar{x}_m$, the value $r_1(\bar{x}_m)$ uniquely determines the derivative $r_{1m}(\bar{x}_m)$. Indeed,

For general nonlinearities, of course, we may choose $r_1(\bar{x}_m)$ and $r_{1m}(\bar{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

This implies that at any fixed $\bar{x}_m$, the value $r_{1m}(\bar{x}_m)$ can be chosen as small as needed, independently from $r_1(\bar{x}_m)$:

In particular, for large $a_1 \rightarrow \infty$ we get small $r_{1m}\rightarrow 0$, and for small $a_1 \rightarrow 0$ we get $r_{1m} \rightarrow \frac{r_1(\bar{x}_m)}{\bar{x}_m}$. Having fixed the value $\bar{x}_m$ requires then to pick $k_1: = \frac{(1+a_1\bar{x}_m)r_1(\bar{x}_m)}{\bar{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 fixed $\bar{x}$ 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 $\mathbf{r}(\bar{x}$)) from the question of its stability (related to the derivatives $r_{jm}(\bar{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_{ij}$ are linear forms in the variables $r_{jm}$.

The leading question of this paper is the following:

That is:

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.

2.   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 \times 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 $\mathbf{J}: \textbf{M} \longrightarrow \textbf{E}$, which associates to every metabolite $m \in \textbf{M}$ a reaction $j \in \textbf{E}$ such that $m$ is an input metabolite of reaction $j$.

We call the reaction $j = \mathbf{J}(m)$ a child of $m$, and the metabolite $m = \mathbf{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, $\mathbf{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^\mathbf{J}$ is the $M\times M$ matrix whose $m^{th}$ column is the $\mathbf{J}(m)^{th}$ column of $S$.

Remark 2. Because of the exclusion of explicit autocatalytic reactions from our model, the stoichiometric coefficient $S_{mj}$ of an input metabolite $m$ to reaction $j$ is a priori negative. This implies that any diagonal element of $S^\mathbf{J}$ is always negative: $S^\mathbf{J}_{mm} = S_{m\mathbf{J}(m)} < 0$ for any Child Selection $\mathbf{J}$ and any metabolite $m$.

Remark 3. If there are no Child Selections at all, then $\operatorname{det}(G)\equiv 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 $\operatorname{det}(G)\equiv 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 homogen{e}ous multilinear polynomial in the variables $r_{jm}$. The possible sign of the polynomial depends on the signs of the monomial coefficients $\operatorname{det}S^\mathbf{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^\mathbf{J}$.

Definition 2 (Child Selection behavior). Let $\mathbf{J}$ be a Child Selection.

We say that $\mathbf{J}$ is bad, or $\mathbf{J}$ $ill-behaves$, if $\operatorname{sign}(\operatorname{det}S^\mathbf{J}) = (-1)^{M-1}$.

We say that $\mathbf{J}$ is good, or $\mathbf{J}$ well-behaves, if $\operatorname{sign}(\operatorname{det}S^\mathbf{J}) = (-1)^M$.

If $\operatorname{det}S^\mathbf{J} = 0$, we say that $\mathbf{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 $\mathbf{J}$, for which all the eigenvalues of $S^\mathbf{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

with $\operatorname{sign}(\operatorname{det}S^\mathbf{J}) = (-1)^{M-1} = (-1)^{3-1} = +1$.

Example 2: good Child Selection

This network represents a good Child Selection:

The stoichiometric matrix associated is

with $\operatorname{sign}(\operatorname{det}S^\mathbf{J}) = (-1)^M = (-1)^3 = -1$.

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.

3.   Main result

In metabolic networks, stoichiometric coefficients are mostly 0 and 1. For this reason, we continue here assuming that $S$ has entries $S_{mj}\in \{-1, 0, 1\}$. By 2, Remark 2, then, the diagonal entries of the reshuffled minor $S^\mathbf{J}$ are $S^\mathbf{J}_{mm} \equiv -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_{mj}\in \mathbb{R}$.

Here, via a structural analysis of $\operatorname{det}S^\mathbf{J}$, we characterize whether the given Child Selection $\mathbf{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^\mathbf{J}$, reads

Again, $\pi$ indicates permutations of $M$ elements and $\operatorname{sgn}(\pi)$ is the signature of $\pi$. Let

denote the summand associated to the permutation $\pi$ in the Leibniz expansion. For example, denoting as ${{Id}}$ the identity permutation,

Let $\pi \neq {{Id}}$ be a permutation such that ${E}(\pi)\neq 0$. Combinatorially, the permutation $\pi$ can be expressed as the product of $\vartheta$ disjoint permutation cycles $c_i$ of length $l_i > 1$,

Definition 3 (good/bad-completions, good/bad-cycles). Given a Child Selection $\mathbf{J}$, we call $\pi$ a good-completion if

We call $\pi$ a bad-completion if

Again, $\vartheta$ is the number of cycles in the permutation expansion. If $\pi$ consists of a single cycle $c$, i.e. for $\vartheta = 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.

Theorem 3.1. Let $\mathbf{J}$ be a Child Selection and let $\mathbf{\mathcal{G}}$ and $\mathbf{\mathcal{B}}$ be the number of good and bad completions, respectively. Then, in the sense of Definition 2,

1. The Child Selection $\mathbf{J}$ well-behaves if $\mathbf{\mathcal{G}} > \mathbf{\mathcal{B}}-1$.

2. The Child Selection $\mathbf{J}$ ill-behaves if $\mathbf{\mathcal{G}} < \mathbf{\mathcal{B}}-1$.

3. The Child Selection $\mathbf{J}$ zero-behaves if $\mathbf{\mathcal{G}} = \mathbf{\mathcal{B}}-1$.

For a Child Selection $\mathbf{J}$, Theorem 3.1 characterizes the sign of $\operatorname{det}S^\mathbf{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.

4.   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_{mj} < 0$ of

in mind, we call negative the edges $e = (m, j)$ where $m$ is input to $j$. Conversely, with $S_{mj} > 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 ($SR$-graph), this was considered by ^[5] and others.

We proceed with two definitions and a proposition.

Definition 4 (J-selected edges). For any Child Selection $\mathbf{J}$, we call the negative edges $e = (m, \mathbf{J}(m))$ in the $MR$-graph $\mathbf{J}$-selected.

Remark 4. In particular, $\mathbf{J}$-selected edges are such that the corresponding stoichiometric entry lies on the diagonal of $S^\mathbf{J}$.

Remark 5. Injectivity of a Child Selection $\mathbf{J}$ directly implies that two $\mathbf{J}$-selected edges $e_1$ and $e_2$ never share a vertex in the $MR$-graph.

Definition 5 (Completion Cycle). For a Child Selection $\mathbf{J}$, a completion cycle in the $MR$-graph is a cycle of length $2l$, $\ell \le M$, such that $\ell$ edges are $\mathbf{J}$-selected.

Remark 6. Equivalently, because of Remark 5 above, a completion cycle is a cycle in the $MR$-graph of length $2l$, such that $\ell$ $\mathbf{J}$-selected edges alternate with $\ell$ non $\mathbf{J}$-selected edges.

Proposition 4.1. For any given Child Selection $\mathbf{J}$, there is a one-to-one correspondence between completion cycles and nonzero permutation cycles; i.e., cycles $c$ such that

It is now clear what the word 'completion' refers to. In fact, any completion cycle is constructed by completing $\ell$ $\mathbf{J}$-selected elements to a cycle of length $2l$ in the $MR$-graph. In this sense, a good-completion $\pi = \prod_{i = 1}^\vartheta c_i$ can be seen, in the $MR$-graph, as a collection of $\vartheta$ non-intersecting completion cycles $c_i$, such that the number of good-cycles has the same parity of $\vartheta$. Respectively, a bad-completion is a collection of $\vartheta$ non-intersecting completion cycles, such that the number of good-cycles has opposite parity of $\vartheta$.

Remark 7. From Definition 3 and Proposition 4.1 it follows that, for a completion cycle in the $MR$-graph, 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.

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.

Corollary 4.2 (Examples of application). The following statements hold true:

1. Acyclic Child Selections well-behave;

2. A Child Selection possessing one good-cycle, and no other completion cycles, well-behaves;

3. A Child Selection possessing one bad-cycle, and no other completion cycles, zero-behaves;

4. A Child Selection possessing two intersecting bad-cycles, and no other completion cycles, ill-behaves;

5. Any nonzero Child Selection of a network which possesses only monomolecular reactions and one single bimolecular reaction

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.

5.   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.

Theorem 5.1 (Change of Stability). Suppose there exist two Child Selections $\mathbf{J}_1$, $\mathbf{J}_2$, and a metabolite $m_b$, such that $\mathbf{J}_1(m_b)\neq \mathbf{J}_2(m_b)$ and $\mathbf{J}_1(m) = \mathbf{J}_2(m)$ for any $m\neq m_b$. Assume moreover that $\mathbf{J}_1$ well-behaves and $\mathbf{J}_2$ ill-behaves. Then the Jacobian determinant of $G$ takes the form

where the omitted terms can be chosen arbitrarily small, and $\mathit{\mathfrak{a}}, \, \mathit{\mathfrak{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 $\xi = \mathit{\mathfrak{a}} r_{\mathbf{J}_1(m_b)m_b}-\mathit{\mathfrak{b}} r_{\mathbf{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.

6.   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 $\mathbf{J}_1$ and $\mathbf{J}_2$ with opposite behavior, such that $\mathbf{J}_1(m_b) \neq \mathbf{J}_2(m_b)$ for one single metabolite $m_b$ and $\mathbf{J}_1(m) = \mathbf{J}_2(m)$ for all other metabolites $m \neq 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 $\mathbf{J}_1$ and $\mathbf{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:

1. $\mathbf{J}_1(PEP) = \mathbf{J}_2(PEP) = dd6$;

2. $\mathbf{J}_1(PYR) = \mathbf{J}_2(PYR) = dd7$;

3. $\mathbf{J}_1(CO2) = \mathbf{J}_2(CO2) = d6$.

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 $\mathbf{J}$ satisfying the above constraints 1-3, identifies reshuffled minors $S^{\mathbf{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 $\mathbf{J}_1 = \mathbf{J}_2$ in the upper part.

In figure 3, we have depicted the chosen subnetwork possessing precisely only the two chosen Child Selections $\mathbf{J}_1$ and $\mathbf{J}_2$. We identify the metabolite $m_b$, such that $\mathbf{J}_1 (m_b) \neq \mathbf{J}_2 (m_b)$, as $ICT = m_b$. Indeed, any other metabolite $m \neq ICT$ possesses a single child $\mathbf{J}_1 (m) = \mathbf{J}_2 (m)$, and only $ICT$ possesses two child reactions: reaction 19 and reaction 26. Let us call $\mathbf{J}_{1}$ the Child Selection such that $\mathbf{J}_{1}(ICT) = 19$ and $\mathbf{J}_{2}$ the Child Selection such that $\mathbf{J}_{2}(ICT) = 26$. With this choice of Child Selections, metabolites Lactate, Acetate, and Ethanol result disconnected from the rest of the network and have consequently been omitted here.

By looking at the $MR$-graph representation, we can easily conclude that $\mathbf{J}_1$ well-behaves and $\mathbf{J}_2$ ill-behaves. Note indeed that $\mathbf{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 $\mathbf{J}_{1}$-selected. On the other hand, the completion cycles structure of $\mathbf{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$:

1. $c_1 = ICT-26-Glyoxylate-27-MAL-23-OAA-17-CIT-18-ICT$;

2. $c_2 = ICT-26-SUC-21-FUM-22-MAL-23-OAA-17-CIT-18-ICT$.

Since it possesses only two intersecting bad completion cycles, the Child Selection $\mathbf{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.

7.   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 homogen{e}ous multilinear polynomial in the variables $r_{jm}$. The coefficients of each monomial are determinants of reshuffled stoichiometric minors $S^\mathbf{J}$, and we have given network conditions to establish the sign of $\operatorname{det}S^\mathbf{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 form

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 $\mathbf{J}_1$ and one bad Child Selection $\mathbf{J}_2$ such that $\mathbf{J}_1(m_b)\neq \mathbf{J}_2(m_b)$ for a single metabolite $m_b$, and $\mathbf{J}_1(m) = \mathbf{J}_2(m)$ for all other metabolites $m\neq 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 $\mathbf{J}$ is any good Child Selection with $\mathbf{J}(PEP) = 2$. Then, it is good also the Child Selection $\tilde{\mathbf{J}}$ such that $\tilde{\mathbf{J}}(PEP) = 3$ and $\tilde{\mathbf{J}}(m) = \mathbf{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.

8.   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 = SR$ to obtain:

Here $\pi$ indicates a permutation of $M$ elements and $\operatorname{sgn}(\pi)$ is the signature (or parity) of $\pi$. Note that $\prod_{m\in \mathbf{M}} r_{\pi(m)m}\neq 0$ if and only if there is an associated Child Selection $\mathbf{J}$ such that $r_{\mathbf{J}(m)m} = r_{\pi(m)m}$, for every $m$. In particular, the sum runs non trivially only for the selected minors $S^\mathcal{E}$ such that the set $\mathcal{E}$ is the image of $\mathbf{M}$ through a Child Selection $\mathbf{J}$. Now,

Last step is the observation:

We proceed with the proof of main technical result, Theorem 3.1.

Proof of Theorem 3.1. The proof follows an idea of Banaji and Craciun ^[7]. Firstly note:

Let $h$ be the number of elements $m$ such that $\pi(m)\neq m$. That is, $h$ is the number of elements of $\pi$ which are not fixed points of the permutation, but belong to a permutation cycle.

The steps above are made noting that $(S^\mathbf{J}_{mm})^2\equiv 1$, for any $m$ and that, for a cycle $c$ of length $\ell$, $\operatorname{sgn}(c) (-1)^\ell = -1$. We conclude the proof by observing that

if $\pi$ 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.

Proof of Proposition 4.1. Let us consider the computation (8.5) in the proof of Theorem 3.1. Consider, for simplicity, a single cycle permutation $\pi = c$ and concentrate on the expression

Note that the diagonal elements $S^\mathbf{J}_{mm}$ and $S^\mathbf{J}_{c(m)c(m)}$ represents $\mathbf{J}$-selected edges. $S_{c(m)m}$ shares the same column (i.e., reaction vertex) with $S^\mathbf{J}_{mm}$ and the same row (i.e., metabolite vertex) with $S^\mathbf{J}_{c(m)c(m)}$. Following the order of the cycle $c$ in the Expression (8.8) leads to the desired identification.

Proof of Corollary 4.2.

1. Acyclic Child Selections do not possess any completion cycle and consequently do not possess any completion. By Theorem 3.1, $\mathcal{G} = \mathcal{B} = 0$, and thus acyclic Child Selections well-behave.

2. Analogously, a Child Selection $\mathbf{J}$ possessing only one good cycle possesses only one good completion. By Theorem 3.1, $\mathcal{G} = 1$, $\mathcal{B} = 0$, and thus $\mathbf{J}$ well-behaves.

3. As in (2), with $\mathcal{G} = 0$, $\mathcal{B} = 1$. Via Theorem 3.1 we have the conclusion.

4. Since a completion is a collection of non-intersecting completion cycles, then two intersecting completion cycles implies two different completions. As above, both completions must be bad. Hence, again by Theorem 3.1, $\mathcal{G} = 0$, $\mathcal{B} = 2$.

5. Note that a non-zero Child Selection $\mathbf{J}$, $\operatorname{det}S^\mathbf{J} \neq 0$, a priori, cannot contain monomolecular cycles of the form:

where each reaction is monomolecular, i.e.

Therefore, for a non-zero Child Selection $\mathbf{J}$ possessing only monomolecular reactions and a single bimolecular reaction $\tilde{j}$ of the form (4.4), either $\mathbf{J}$ does not contain any completion cycle and consequently well-behaves, due to point (1) of this corollary, or $\mathbf{J}$ contains completion cycles in the $MR$-graph which include $\tilde{j}$ as a reaction vertex. In particular all completion cycles intersects.

Without loss of generalities, assume $\mathbf{J}^{-1}(\tilde{j}) = A$. We only have two possibilities for a completion cycle: a bad completion cycle containing the two adjacent edges $A-\tilde{j}-C$ or a good completion cycle containing the two adjacent edges $A-\tilde{j}-B$. By assumption on the network, the completion cycles run through monomolecular reactions only, except from $\tilde{j}$. Because they do intersect, we have at most one bad completion, $\mathcal{G}\le 1$, and one good completion, $\mathcal{B} \le 1$. Hence, $1-\mathcal{G}+\mathcal{B} \ge 0$ and by Theorem 3.1 the Child Selection $\mathbf{J}$ does not ill-behave. Since, by assumption, $\mathbf{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 $\{\mathbf{J}\}$ carries a natural integer-valued distance $d$.

Definition 6. Let $\mathbf{J}_1$, $\mathbf{J}_2$ be two Child Selections. We define the distance $d(\mathbf{J}_1, \mathbf{J}_2)$ as the number of metabolites $m \in \mathbf{M}$ such that $\mathbf{J}_1(m) \neq \mathbf{J}_2(m)$.

It is straightforward to verify that $d$ is a distance on the set of Child Selections. We consider now Child Selections at distance $d = 1$. These are Child Selections $\mathbf{J}_1$, $\mathbf{J}_2$ such that $\mathbf{J}_1(m_b)\neq \mathbf{J}_2(m_b)$ for a single metabolite $m_b$ and $\mathbf{J}_1(m) = \mathbf{J}_2(m)$ for any $m\neq m_b$ different from $m_b$, as in Theorem 5.1. Clearly:

If we further assume that $\mathbf{J}_1$ and $\mathbf{J}_2$ are such that one well-behaves and the other ill-behaves we have:

with $\mathit{\mathfrak{a}}$ and $\mathit{\mathfrak{b}}$ constants of the same sign.

By the mere fact that $d$ is an integer-valued distance, any other Child Selection satisfies

In particular, we have the following Lemma:

Lemma 8.1. Let $\mathbf{J}_1$ and $\mathbf{J}_2$ be Child Selections at distance $d = 1$, that is, $\mathbf{J}_1(m_b)\neq \mathbf{J}_2(m_b)$ and $\mathbf{J}_1(m) = \mathbf{J}_2(m)$ for any $m\neq m_b$. For any other Child Selection $\mathbf{J}_k$, there is a metabolite $m_k$ such that $\mathbf{J}_k(m_k)\neq \mathbf{J}_1(m_k)$ and $\mathbf{J}_k(m_k)\neq \mathbf{J}_2(m_k)$.

Moreover if $d(\mathbf{J}_k, \mathbf{J}_1) = d(\mathbf{J}_k, \mathbf{J}_2) = 1$, then $m_k = m_b$.

Proof. Let us consider any $m_k$ such that $\mathbf{J}_1(m_k)\neq \mathbf{J}_k(m_k)$. If $\mathbf{J}_2(m_k)\neq \mathbf{J}_k(m_k)$, we are done. Assume then that $\mathbf{J}_2(m_k) = \mathbf{J}_k(m_k)$. By construction, $m_k = m_b$. Consider now any $\tilde{m}_k$ such that $\mathbf{J}_2(\tilde{m}_k)\neq \mathbf{J}_k(\tilde{m}_k)$ and remember that $\mathbf{J}_1(m) = \mathbf{J}_2(m)$ for any $m\neq m_b$. We conclude that $\mathbf{J}_1(\tilde{m}_k)\neq \mathbf{J}_k(\tilde{m}_k)$. Otherwise, we would have found two metabolites $m_k$ and $\tilde{m}_k$ such that $\mathbf{J}_1(m_k)\neq \mathbf{J}_2(m_k)$ and $\mathbf{J}_1(\tilde{m}_k)\neq \mathbf{J}_2(\tilde{m}_k)$, contradicting $d(\mathbf{J}_1, \mathbf{J}_2) = 1$.

In the above argument, note that if $\mathbf{J}_2(m_k) = \mathbf{J}_k(m_k)$, then $d(\mathbf{J}_1, \mathbf{J}_k)\ge2$. Hence, if $d(\mathbf{J}_1, \mathbf{J}_k) = d(\mathbf{J}_2, \mathbf{J}_k) = 1$ we conclude that $\mathbf{J}_1(m_b)\neq \mathbf{J}_2(m_b)\neq \mathbf{J}_k(m_b)$ and hence $m_k = m_b$.

We are now ready to prove Theorem 5.1.

Proof of Theorem 5.1. Let $\mathbf{J}_1$ and $\mathbf{J}_2$ be Child Selections, as above. In particular, distance $d(\mathbf{J}_1, \mathbf{J}_2) = 1$. By Lemma 8.1, for any other Child Selection $\mathbf{J}_k \neq \mathbf{J}_1, \mathbf{J}_2$ we can find $m_k$ such that

We can consider, then, an $\epsilon$-small choice of reaction rate parameters such that

Then, for this $\epsilon$-choice of reaction rates,

The parameter $\xi = \mathit{\mathfrak{a}} r_{\mathbf{J}_1(m_b)m_b}-\mathit{\mathfrak{b}} r_{\mathbf{J}_2(m_b)m_b}$ becomes then a bifurcation parameter for the sign of the Jacobian determinant.

Acknowledgments

We thank the anonymous referees for their helpful and valuable comments. We are very grateful to Bernold Fiedler for his enduring support. Our thanks go also to the present and past group of Nonlinear Dynamics of Freie Universität Berlin. Anna Karnauhova with everlasting gentleness allowed us to use her masterful picture of the central metabolism of E. coli. Meghan Kane greatly helped in fixing non-native english typos. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-project number 163436311-SFB 910, and supported by the Berlin Mathematical School.

Conflict of interest

The authors have no conflict of interest.