GRAPH-THEORETIC CONDITIONS FOR ZERO-EIGENVALUE TURING INSTABILITY IN GENERAL CHEMICAL REACTION NETWORKS

We describe a necessary condition for zero-eigenvalue Turing instability, i.e., Turing instability arising from a real eigenvalue changing sign from negative to positive, for general chemical reaction networks modeled with mass-action kinetics. The reaction mechanisms are represented by the species-reaction graph (SR graph), which is a bipartite graph with different nodes representing species and reactions. If the SR graph satisfies certain conditions, similar to the conditions for ruling out multiple equilibria in spatially homogeneous differential equations systems, then the corresponding mass-action reaction-diffusion system cannot exhibit zero-eigenvalue Turing instability for any parameter values. On the other hand, if the graph-theoretic condition for ruling out zero-eigenvalue Turing instability is not satisfied, then the corresponding model may display zero-eigenvalue Turing instability for some parameter values. The technique is illustrated with a model of a bifunctional enzyme.

1. Introduction.Spatially homogeneous concentrations of chemical species are assumed in many dynamic models from cell biology.However, spatial diffusion processes are important for the proper functioning of a cell; for example, for the development, polarity and the formation of biological patterns referred to as morphogenesis [20].Therefore, cells can be viewed as reactors that are not well stirred, where chemical species diffuse within and between the different compartments of a cell.
Graph-theoretic methods are useful in biological applications since parameters such as rate constants and diffusion coefficients are not entirely known for complex chemical reaction networks.The relationship between the network's structure and different types of behavior such as multistability, oscillations or Turing instabilities is crucial for understanding the properties of reaction mechanisms.Realistic chemical reaction networks involve a large number of substances (genes, proteins, 1208 MAYA MINCHEVA AND GHEORGHE CRACIUN metabolites, signaling molecules) and reactions, creating large chemical (genetic, metabolic, signaling) networks [25].Therefore, it is necessary to develop graphtheoretic methods applicable to complex chemical reaction networks.
Studying Turing instabilities in biochemical models involves the analysis of a nonlinear system of reaction-diffusion equations.The first mathematical model for pattern formation that uses a reaction-diffusion system was proposed by Alan Turing in 1952.It was shown that patterns develop when an asymptotically stable spatially homogeneous equilibrium in the absence of diffusion, becomes unstable in the presence of diffusion [26].In this work, the diffusion-driven instability is referred to as zero-eigenvalue Turing instability, because it is due to a real eigenvalue crossing the imaginary axis from left to right.
In the absence of diffusion, a chemical reaction network with mass action kinetics gives rise to an ordinary differential equation (ODE) model where the rate constants associated with the network are the parameters.A chemical reaction network with species outflows or degradations is said to be injective if its ODE model has a negative Jacobian whose determinant is positive for all parameter values.If a chemical reaction network is injective, then it does not have the capacity for multistability, or, the existence of multiple positive equilibria [5].We define principal subnetworks (Definition 3.1) of a chemical reaction network, where the injectivity of a principal subnetwork is defined similarly as for a reaction network.We show that if there exists a non-injective principal subnetwork whose species set does not include all species, then the corresponding reaction-diffusion model (defined in Section 6) taken with mass action kinetics can exhibit zero-eigenvalue Turing instability for some parameter values (rate constants and diffusion coefficients).
Chemical reaction networks are often represented by various types of graphs [7,9,17,18,24].The mass-action kinetics reaction network studied here is represented by a bipartite graph with different types of nodes for the species and the reactions, and is referred to as the species-reaction (SR) graph [6].The bipartite graph representation is especially convenient for multi-molecular mass-action networks, since it allows for different paths between the same species nodes going through different reaction nodes.An inspection of the SR graph shows if the corresponding reaction-diffusion model can exhibit zero-eigenvalue Turing instability for some values of the parameters.
The ordinary differential equation counterpart to this problem has been studied elsewhere [5,6,18,27] and many of the results therein are essential for understanding the theory in this paper.Similar graph-theoretic conditions for zero-eigenvalue Turing instability have been obtain in [19,27].
In Section 2 some preliminaries about the ordinary differential equation model relevant to the corresponding reaction-diffusion system are explained.Principal subnetworks of a chemical reaction network are defined in Section 3 and the conditions for their injectivity are given in Section 4. In Section 5 the main ideas of the graphtheoretic analysis are introduced.In Section 6 the corresponding reaction-diffusion model of a reaction network and the graph-theoretic condition for zero-eigenvalue Turing instability are presented.A model of a bifunctional enzyme is analyzed for zero-eigenvalue Turing instability using the graph-theoretic condition in Section 6.
2. Reaction networks.We give a definition of a reaction network in terms of the set of species, the set of complexes, and the set of reactions.Then we define the corresponding ordinary differential equations model with general mass action kinetics.
Let R denote the set of real numbers, R + the set of positive numbers, and R+ the set of nonnegative numbers.Given a set I, R I denotes the vector space of linear combinations i∈I λ i i, generated by the elements of i ∈ I, with coefficients λ i ∈ R.

By RI
+ we denote the set of members of R I with λ i ≥ 0 for all i ∈ I.By R I + we denote set of the members of R I with λ i > 0 for all i ∈ I.The support of an element x ∈ R I + forms the set supp(x) = {i ∈ I : x i = 0}.The complexes of a reaction network are to be understood as the objects at the head or tail of reaction arrows.For example, the complexes of the reaction A + B → C are A + B and C. We write for each each complex y = s∈S y s s where y s ≥ 0 if s ∈ S .We will write y → y in place of (y, y ) when (y, y ) is a member of R. Also, if {y → y , y → y} ⊂ R we will denote the set {y → y , y → y} by y y and will say that y y is a reversible reaction.If y → y ∈ R and y → y / ∈ R, we say that y → y is an irreversible reaction.For example, consider the reaction network In this case S = {A, B, C}, The first reaction in (1) is reversible and the second is irreversible.Let c ∈ RS + denote the vector of species concentrations from S .Given a reaction network N next we define some special classes of dynamical systems associated with N .
Definition 2.2.A kinetics K for a reaction network N = (S , C , R) is an assignment to each reaction y → y ∈ R of a rate function K y→y : RS + → R+ .Definition 2.3.A kinetic system (S , C , R, K ) is a reaction network (S , C , R) taken together with a kinetics K for the network.
In the next definition we use the following notation: for two vectors in RS + , u = s∈S u s s and v = s∈S v s s, we denote by u v the product s∈S (u s ) vs with the assumption that 0 0 = 1.
Next we define a special type of kinetics referred to as a general mass action kinetics.Similar generalizations of mass action kinetics have been discussed in [15,23,29].Definition 2.4.A kinetics K for a reaction network (S , C , R) is general mass action if for each reaction y → y ∈ R there is a nonempty finite set V y→y ⊆ R S such that K y→y (c) = where k v,y→y ∈ R + for all v ∈ V y→y .The elements v of the set V y→y are called kinetic orders of y → y and the constants k v,y→y ∈ R + are called reaction rate constants.
If the set V y→y = {y} for all y → y ∈ R, then the rate functions K y→y (c) = k y,y→y c y and K is called mass action kinetics [10].In this case we denote k y→y = k y,y→y , and we say that a rate function K y→y (c) = k y→y c y is of mass action type.For example, the rate function K y→y = k A+B→C c A c B corresponding to the first reaction in example ( 6) is of mass action type.If instead of y = (1, 1, 0) the kinetic order v = (1, 0, 0) is used, then the rate function becomes When the use of mass action kinetics is not suitable in applications, various generalizations may be used [15,29].We will use general mass action kinetics to define the kinetics of special reaction subnetworks (see Definition 3.1 below) of a reaction network N with mass action kinetics.Definition 2.5.A general mass action system (S , C , R, K ) is a reaction network (S , C , R) taken together with general mass action kinetics K for the network.
Given a general mass action system, the differential equation system that governs the evolution of the species concentration vector c ∈ RS + is defined next.
We will denote the right-hand side of the system (3) with general mass action kinetics by In the case of mass action kinetics, the associated differential equation system (3) in vector form becomes for some vector k ∈ R R + of rate constants k y→y .
Definition 2.7.If c * ∈ R S + satisfies r(c * ) = 0 of a general mass-action system (3), then c * is a positive equilibrium of (3).Let 0 be the zero complex, which is understood to be the zero vector of RS + .As explained in [5], outflow reactions such as A → 0 model the contributions of the outflow stream or species degradation, while feed reactions such as 0 → A model the contributions of the feed stream or species production.
Similarly to [5], in this paper we study the "augmented network" that contains feed and outflow reactions (also called flow reactions) for all species.For example, the augmented network of (1) is The differential equations system with mass action associated with the network (6) is Assumption.From now on, it will be assumed that the reaction network (S , C , R) is the augmented network with flow reactions for each species s ∈ S .Then, the set of reaction R includes the set of true reactions R t , the set of outflow reactions R o and the set of feed reactions R f , i.e., 3. Principal subnetworks.In this section we define a special type of subnetworks of a reaction network N , referred to as principal subnetworks.We will show in Section 6 that the existence of a non-injective principal subnetwork (injectivity is defined in Section 4) in the SR graph of a chemical reaction network is a necessary condition for zero-eigenvalue Turing instability.
In order to form a principal subnetwork we consider a chemical reaction network such as (6) and eliminate some of the species from the true reactions by changing their corresponding coefficients from positive to zero.For example, if species A is eliminated from reaction (6) then the principal subnetwork is Note that the flow reactions in (8) remain the same.Therefore, the species set of a principal subnetwork will be the same as the species set of the corresponding reaction network N .However, only a subset of species from S will participate in the true reactions, while all species will take part in the flow reactions.For the reaction network (6) and its principal subnetwork (8), species B and C take part in the true reactions, while species A takes part only in the flow reactions.
The complexes of a principal subnetwork will be defined in the same way as the complexes of the original reaction network.For example, the complexes of the reaction B → C from ( 8) are B and C, and in vector form (0, 1, 0) and (0, 0, 1), respectively.
The true reactions of the original network N will be kept formally as such for a principal subnetwork, even if they take the form of a flow reaction or the reaction 0 → 0. For example, the reaction 0 → 2B in ( 8) is considered a true reaction of the subnetwork (8).
Next we give a formal definition of a principal subnetwork of a reaction network N .Definition 3.1.Consider a network N = (S , C , R) with species set S = {s 1 , s 2 , . . ., s n }.Let α l = {i 1 , . . ., i l } ⊆ {1, 2, . . ., n} and the set S α l = {s i1 , s i2 , . . ., s i l }, and denote its complement by S c α l .A principal subnetwork N α l = ( S , C , R) consists of: (i) the set of species S = S , (ii) the set of complexes C which consists of all complexes of flow reactions and all complexes obtained by projecting on the span of S α l the complexes in C that appear in true reactions from R t , (iii) the set R which consists of all flow reactions and all the projections on the span of S α l of true reactions from R t .
If a complex y ∈ C , then we denote its projection on S α l by ȳ ∈ C .If y → y is a reaction in R t we denote the projected reaction by ȳ → ȳ ∈ Rt .We will use general mass action kinetics (Definition 2.4) in order to define the kinetics associated with a principal subnetwork N α l .This is because we would like the differential equations for the concentrations of the species s ∈ S α l to be the same as the corresponding mass action equations ( 5) for the original reaction network N .Therefore, the rate function of a projected true reaction ȳ → ȳ ∈ Rt needs to be the same as for the corresponding reaction y → y ∈ R t .Let p be the projection map for all true reactions from R t on the span of S α l and the identity map for all flow reactions.We will denote the forward map by p(y → y ) = ȳ → ȳ and the inverse map by Then, the rate function of a reaction ȳ → ȳ from N α l using general mass action kinetics ( 2) is where the sum is indexed by reactions y → y such that ȳ → ȳ is the projection of reactions y → y on the span of S α l .The sum in ( 9) is necessary since different true reactions from R t can project onto the same reaction from Rt (see Example 3.2.below).Therefore, the general mass action system of a principal subnetwork Equations in (10) that correspond to a species s ∈ S c α l have the simple form ċs = −k s→0 c s + k 0→s since s participates in flow reactions only.Equations in (10) that correspond to a species s ∈ S α l have the form ċs = y→y ∈R k y→y c y (y s − y s ), which is exactly the same equation as the corresponding equation in (5) for the original reaction network N .
Example 3.1.Consider the reaction network from (6) and suppose that α 2 = {1, 2}.Then S α2 = {A, B} and the corresponding principal reaction subnetwork is The set of complexes of ( 11) is C = {A + B, A, 2B, 0} and the set of its true reactions is Rt = {A + B → 0, 0 → A + B, A → 2B}.We choose the rate constants for the projected reactions from (11) to be the same as for the corresponding reactions from (6).Mass action kinetics rate functions are used, except for the rate function of 0 → A + B where general mass action kinetics is used with kinetic order set V 0→A+B = {(0, 0, 1)} leading to the rate function k C→A+B c C .Therefore, the system of differential equations of the principal subnetwork (11) in vector form is and componentwise is Example 3.2.Consider the reaction network ( 6) and now suppose that α 1 = {1}.
Then S α1 = {A}, and the corresponding principal subnetwork N α1 is The set of complexes of ( 13) is {A, 0} and the set of its true reactions is Rt = {A → 0, 0 → A}.The rate constants for the projected reactions from ( 13) are chosen to be the same as for the corresponding reactions from (6), except for the true reactions A + B → C and A → 2B that project to A → 0. Let the kinetic order set V A→0 = {(1, 1, 0), (1, 0, 0)}.Then the rate constant for the kinetic order v = (1, 0, 0) is k A→2B + k A→0 and the rate function is (k A→2B + k A→0 )c A since A → 2B also projects to A → 0. Also, for the kinetic order v = (1, 1, 0) the rate function is Similarly, for the kinetic order set V 0→A = {(0, 0, 1), (0, 0, 0)}, the corresponding rate function is and by ( 9), we have The same system of differential equations, written componentwise, is 4. Injective principal subnetworks.In this section we will use properties of the Jacobian matrix of the right-hand side r(c, k) of the ODE system (5) to introduce chemical reaction networks referred to as injective networks [5].Then we will define injective principal subnetworks in a similar way.In Section 6 we will show that injective networks cannot exhibit zero-eigenvalue Turing instability for any parameter values.
Let {e 1 , . . ., e n } be the canonical basis in R n .Then, as described in [5], the Jacobian matrix associated with the system (5) can be written as is a weighted scalar product.We will deal often with the negative Jacobian (15) For example, the negative Jacobian matrix of the chemical network (6) with ODE model ( 7) is The notion of injective network was introduced in [5].One of several equivalent characterizations of injectivity described in [5] is given in the following definition.
for all values of the rate constants k ∈ R R + and the concentrations c ∈ R S + .Equivalently, it follows by Definition 4.1 that a reaction network is injective if det(J − (c, k)) = 0.
Injective networks cannot have multiple positive equilibria, i.e., cannot exhibit multistability, for any values of the reaction rate parameters k ∈ R R + [5].It is shown in [5] that any coefficient in the expansion of the determinant of the negative Jacobian (16) equals the product of two determinants det(y 1 , .., y n ) det(y 1 −y 1 , .., y n −y n ) for some choice of n reactions y 1 → y 1 , .., y n → y n , and that a reaction network is injective if and only if all such products of determinants are non-negative.In other words, another characterization of the injectivity property of a reaction network is given in the following theorem.
for all choices of n reactions (some of which could be flow reactions) y 1 → y 1 , ..., y n → y n in R.
By (16) it follows that, the determinant of the negative Jacobian det(J − (c, k)) is a polynomial over the variables given by the rate constants k yi→y i ∈ R + and the concentrations c j ∈ R + (see also (17)).In particular, since there is at least one coefficient (19) in the expansion of det(J − (c, k)) which is strictly positive (obtained by choosing all reactions that appear in (19) to be outflow reactions), it follows that if the hypothesis of Theorem 4.1 holds, then det(J − (c, k)) > 0 for any k ∈ R R + and c ∈ R S + .The negative Jacobian matrix of the general mass action model (10) of a principal subnetwork N α l can be obtained similarly to (16) and equals for some choice of n reactions y i → y i , = 1, . . ., n from R (some of which can be flow reactions) projected on corresponding reactions ȳi → ȳ i from R. We define an injective principal subnetwork N α l , similarly to an injective network N .If the determinant of its Jacobian, or equivalently its negative Jacobian (20), does not vanish for all values of the rate constants k ∈ R R + and the concentrations c ∈ R S + , then N α l is injective.The next corollary gives an alternative definition for injectivity of principal subnetworks.The proof follows the proof of [5,Theorem 3.3].Recall that if A is an (n × n) matrix and A[α l ] denotes the submatrix of A with rows and columns indices α l ⊆ {1, . . ., n}, l = 1, . . ., n, then det(A[α l ]) is the corresponding principal minor of order l [12].Next we show that the determinant of the negative Jacobian (20) for a principal subnetwork N α l is a positive multiple of the corresponding principal minor det(−J[α l ]) of the negative Jacobian (16).Lemma 4.4.Let α l be an arbitrary subset of {1, . . ., n} for some l ∈ {1, . . ., n}.Any principal minor of the negative Jacobian ( 16) can be written as where J − α l (c, k) is the negative Jacobian ( 20) associated with the model system (10) of a principal subnetwork N α l .
Proof.Let α l be an arbitrary subset of {1, . . ., n} for some l ∈ {1, . . ., n} and det(J − [α l ])(c, k) be the corresponding principal minor of the negative Jacobian (16).Consider the negative Jacobian (20) and suppose that after row and column reordering the first l rows and l columns are occupied by rows and columns with indices from the set α l (this requires an even number of switches of rows and columns, so it produces no change of sign ).The remaining (n − l) rows of J − α l (c, k) have a non-zero entry only on the diagonal, and i ∈ {n − l + 1, . . ., n}.Thus the matrix J − α l (c, k) has a block structure that implies that its determinant equals the product of the determinants of its two diagonal blocks, where its (l × l) block is exactly J − [α l ](c, k).Therefore equality (22) holds.
Theorem 4.5.Consider some reaction network N = (S , C , R).Then, the following three conditions are equivalent: (i) N is injective, (ii) the principal subnetworks N α l are injective for all α l ⊆ {1, . . .n}, l = 1, . . ., n, (iii) the negative Jacobian matrix J − (c, k) is a P -matrix for all values of the rate constants k ∈ R R + and all concentrations c ∈ R S + .

Proof. (i) ⇒ (ii):
Consider some α l ⊆ {1, . . ., n}, and n reactions Then after possibly renaming the species, we can assume that the last n−l reactions must be flow reactions (otherwise some column of the matrix above would be zero).Therefore, by using linear combinations of these last n − l columns, we obtain Since N is injective, we have The SR graph of the reaction networks ( 1) and ( 6).
Therefore the principal subnetwork N α l is injective.

5.
The SR graph of a principal subnetwork.The SR graph (species-reaction graph) of a reaction network N was first defined in [4,6] and is a bipartite graph where the nodes represent either species or reactions, and the edges encode information about which species participate in which reactions.We modify the definition of a SR graph for a principal subnetwork N α l , α l ⊆ {1, . . ., l}, l ∈ {1, . . ., n}.Definition 5.1.Consider some principal reaction subnetwork N α l = ( S , C , R), α l ⊆ {1, . . ., l}, l ∈ {1, . . ., n}.The SR graph Γ Nα l of N α l is an unoriented graph where each node of Γ Nα l is either a species node or a reaction node.There is one species node for each species in S α l , there is one reaction node for each reversible reaction in the set of true reactions Rt , and there is one reaction node for each irreversible reaction in Rt .Each edge of the graph Γ Nα l connects a species node to a reaction node, as follows.Consider a species node s and a reaction node r given by ȳ → ȳ or ȳ ȳ .If s ∈ supp(ȳ) then there is an edge between s and r and we label it with the complex ȳ.Similarly, if s ∈ supp(ȳ ) then there is an edge between s and r and we label it with the complex ȳ .
Remark 5.1.(a.)If a species s is contained in both supp(ȳ) and supp(ȳ ) (as in A + B → 2A), then the two edges joining the species node s to the reaction node ȳ → ȳ are labeled with ȳ and ȳ , respectively.(b.)If l = n, then N αn = N and its SR graph is the SR graph of the reaction network Γ N as defined in [4,6].
The SR graph of the reaction network ( 6) is shown in Figure 1.Chemical reaction network properties, such as multistability, depend on the presence of configurations of edges and cycles in the SR graph that are especially important, and they are described in the definition below.Definition 5.2.Consider the SR graph Γ Nα l of some principal reaction subnetwork N α l .A pair of edges in Γ Nα l that meet at a reaction node and have the same complex label is called a c-pair.A cycle that contains an odd number of c-pairs is called an o-cycle.A cycle that contains an even number of c-pairs is called an e-cycle.The stoichiometric coefficient of an edge is the coefficient of the species adjacent to that edge in the complex label of the edge.An s-cycle is one for which, if we alternately multiply and divide the stoichiometric coefficients of edges along the cycle we obtain for the final result 1.An S-to-R chain in an SR graph is a simple path from a species node to a reaction node.We say that two cycles in Γ Nα l have an S-to-R intersection if their common edges constitute an S-to-R chain or a disjoint union of two or more S-to-R chains.
The SR graph Γ Nα l of a principal subnetwork N α l can be obtained from the SR graph Γ N of the corresponding reaction network N by first removing all species s ∈ S c α l and the edges connected to them, and then removing all isolated reaction nodes.The reaction nodes of a principal subnetwork are relabeled to include only species s ∈ S α l .Similarly, the complexes over the edges are relabeled so that they include only species s ∈ S α l .For example, the SR graph of the principal subnetwork ( 11) can be obtained from the SR graph of the reaction network (6) shown in Figure 1 by removing the species node C and the edge attached to it, and then relabeling the reaction node A + B C by A + B 0. The main result in [6, Theorem 6.1] describes a sufficient condition for injectivity in terms of the SR graph.The statement of the theorem below is a modification for a principal subnetwork of a reaction network.The proof follows the proof of [6, Theorem 6.1].
Theorem 5.1.Consider some principal reaction subnetwork N α l , α l ⊆ {1, . . ., n}, l = 1, . . ., n of a reaction network N , such that in its SR graph Γ Nα l all cycles are o-cycles or s-cycles, and no two e-cycles have an S-to-R intersection.Then the principal reaction subnetwork N α l is injective.
The next corollary follows from Theorem 4.5.
Corollary 5.2.Consider some reaction network N such that in its SR graph Γ N all cycles are o-cycles or s-cycles, and no two e-cycles have an S-to-R intersection.Then all of its principal subnetworks N α l , α l ⊆ {1, . . ., n}, l = 1, . . ., n are injective.
In the next section we will show that the graph-theoretic condition described in Theorem 5.1 and Corollary 5.2 is also related to the capacity of a reaction network to exhibit zero-eigenvalue Turing instability.
6. Zero-eigenvalue Turing instability.In this section we show how the theory from Section 4 and Section 5 can be applied to the corresponding reaction-diffusion model of the system (5).We are particularly interested in a graph-theoretic condition that rules out zero-eigenvalue Turing instability which can arise when diffusion is included in the model.On the other hand, if the graph-theoretic condition is violated, we can infer that zero-eigenvalue Turing instability can occur for some values of the parameters.
If the concentration vector c = c(x, t) contains spatially non-homogeneous functions, the corresponding reaction-diffusion system to (5) with initial condition and no-flux boundary condition is where The initial condition functions satisfy c (0) (x) ≥ 0 for all x ∈ Ω, since they represent concentrations.
The no-flux boundary condition (27c) is considered standard for systems studied for Turing instabilities [20].However, in order to observe a Turing pattern in a chemical system, it is necessary for the system to be open to mass flow.Therefore, it is assumed that the mass flow enters through the top or bottom boundary and not through the side boundary.
We say that Turing instability occurs for the system (27a) with no-flux boundary condition (27c) if a spatially homogeneous equilibrium c * (k) > 0 is linearly asymptotically stable as a solution of the ordinary differential equation system (5) and unstable as a solution of the corresponding reaction-diffusion system [3].If, in addition, this instability occurs as a result of a real eigenvalue crossing the imaginary axis from left to right, then we say that zero-eigenvalue Turing instability occurs.
The following proposition and corollary are standard in textbooks on advanced linear algebra and can be found, for example, in [12,16].Proposition 6.1.The coefficients āl (c, k), l = 1, . . ., n, of (28) are sums of all principal minors of order l of the negative Jacobian (16), If the Jacobian matrix J(c, k) has only eigenvalues with negative real parts, then J(c, k) is said to be stable.If J(c, k) has at least one eigenvalue with a positive real part, then J(c, k) is unstable.
The next corollary is a necessary condition for a stable matrix and its proof can be found in [12].Corollary 6.2.If the Jacobian J(c, k) is stable, then all coefficients of (28) are positive, āl (c, k) > 0, l = 1, . . ., n.
The problem of Turing instability is related to the problem of stability of the matrix J(c, k)−µ j D, where µ j ≥ 0, j ∈ N is an eigenvalue of the negative Laplacian [19].The eigenvalues µ j , j ∈ N depend on the set Ω and form an increasing sequence {µ j } where µ 1 = 0 and µ j+1 > µ j for j > 1.For convenience we let µ = µ j > 0 and µ will be treated as a continuous variable.If for µ = 0, the matrix J(c, k) − µD is stable and for some eigenvalue of the negative Laplacian µ = µ j > 0, the matrix J(c, k) − µD becomes unstable, as a single real eigenvalue changes sign from negative to positive, then zero-eigenvalue Turing instability occurs [3].For large µ 0, assuming the diffusion coefficients d j , j = 1, . . ., n are kept fixed, it follows by Geršgorin's theorem [16] that the matrix J(c, k) − µD is also stable.
It follows by the assumption for zero-eigenvalue Turing instability that the Jacobian J(c, k) is stable for some equilibrium c = c * (k) ∈ R S + and some values of the rate constants k ∈ R R + .We can also assume that the diffusion coefficients d j > 0, j = 1, . . ., n are arbitrary but fixed, since this will not alter the results presented below.Therefore, we can assume that the only parameter involved in the matrix Let the characteristic polynomial of the matrix J(c, k) − µD be where the coefficients a m = a m (µ), m = 1, . . ., n.
The next proposition will be used in the proof of Theorem 6.4 below.Then all principal minors of J − (c, k)D −1 are a positive multiple of the principal minors of J − (c, k).Therefore, J − (c, k)D −1 is a P -matrix if and only if J − (c, k) is a P -matrix.
We say that a reaction network N has the capacity for zero-eigenvalue Turing instability, if the reaction-diffusion model (27a)-(27c) can display zero-eigenvalue Turing instability for some parameter values of µ > 0. have to be chosen small (so that | det(J − (c, k)D −1 )[α l ]| is large) and the remaining diffusion coefficients have to be chosen large (making the positive principal minors in b l (µ) small) so that b l (µ) < 0. This condition is another generalization of the condition for two-species systems that the ratio of the diffusion coefficients should be much larger than one [28].
Example 6.1.A model for the bifunctional enzyme phosphofructo-2-kinase:fructose-2,6-bisphosphatase which is involved in the glycolysis (gluconeogenesis) switch [13,14] was studied in [19] for zero-eigenvalue Turing instability.We show here that the chemical reaction network given in (31) below has the capacity for zeroeigenvalue Turing instability using the graph-theoretic condition from Corollary 6.5.By Proposition 6.6 it follows that if the reaction network has the capacity for zero-eigenvalue Turing instability, then a non-injective principal subnetwork that does not include all species exists.By Theorem 5.1 the SR graph of a non-injective principal subnetwork contains either an e-cycle that is not an s-cycle or contains e-cycles with a S-to-R intersection.
We write the augmented chemical network of the bifunctional enzyme with flow reactions using the notation in this article: The SR graph of the chemical reaction network (31) is shown in Figure 2.
The mass-action kinetics reaction-diffusion system taken with nonnegative initial condition and no-flux boundary condition is The SR graph of the reaction network (31) contains e-cycles with S-to-R intersections.For example, the e-cycle through species nodes A 1 and A 4 and reaction nodes A 4 + A 6 → A 1 + A 2 and A 1 + A 2 → A 1 + A 4 , and the e-cycle through the species node A 1 and the reaction node A 1 + A 2 → A 1 + A 4 intersect at the edge starting at species node A 1 and ending at reaction node A 1 +A 2 → A 1 +A 4 .Therefore, by Corollary 6.5, the reaction network (31) can display zero-eigenvalue Turing instability for some parameter values.Moreover, Proposition 6.6 tells us that the reaction network (31) should contain at least one non-injective principal subnetwork that does not include all seven species, if zero-eigenvalue Turing instability exists for some parameter values.Two non-injective principal subnetworks that can lead to zero-eigenvalue Turing instability in the reaction-diffusion model of (31) are A 4 + A 6 → 0, A 4 + A 5 0, A 5 → A 5 + A 6 , 0 → A i → 0, i = 4, 5, 6. (33b) The SR graphs of the principal subnetworks (33) are shown in Figure 3. Clearly, the principal subnetworks from (33) contain a S-to-R intersection between e-cycles, and thus are not injective.
Turing patterns are obtained as a result of a zero-eigenvalue Turing instability for some values of the parameters in [19].Also, note the similarity between the principal subnetworks in Figure 3 and the subnetworks referred to as critical fragments in Figure 3 in [19], given as the structural reason for the zero-eigenvalue Turing instability.7. Discussion.In this article we have obtained a necessary graph-theoretic condition for zero-eigenvalue Turing instability in reaction-diffusion models of mass action kinetics networks.This graph-theoretic condition is easy to use and applies to general reaction networks with any number of species.Related algebraic conditions for Turing instabilities in larger reaction-diffusion models are discussed, for example, in [21,22].Using the obtained graph-theoretic condition, we are able to find specific structures in the SR graph of a reaction network, referred to as noninjective principal subnetworks (see Corollary 4.3 and the definition above it), that are necessary for the existence of zero-eigenvalue Turing instability.The graph-theoretic criterion obtained in this article is related to the graphtheoretic condition for zero-eigenvalue Turing instability in [19,27], but is easier to

Definition 2 . 1 .
([10, 11]) A chemical reaction network N = (S , C , R) consists of three finite sets: (i) a set S of species of the network, (ii) a set C ⊂ RS + of complexes of the network, (iii) a set R ⊂ C × C of reactions, with the following properties:(a) (y, y) / ∈ R for any y ∈ C , (b) for each y ∈ C there exists y ∈ C such that (y, y ) ∈ R or such that (y , y) ∈ R.

Definition 2 . 6 .
For a general mass action system (S , C , R, K ) the associated differential equation system is ċ = y→y ∈R K y→y (c)(y − y) = y→y ∈R v∈V y→y k v,y→y c v (y − y).

k
y→y c y (y * e 1 )(y − y), . . ., y→y ∈R k y→y c y (y * e n )(y − y) where y * e i = k

Proof. Since [ 5 ,
Theorem 3.2] holds for any matrix M of the form M (c, k) =   y→y ∈R k y→y c y (y * e 1 )w y,y→y , . . ., y→y ∈R k y→y c y (y * e n )w y,y→y   , (21) where the constant vectors w y,y→y are arbitrary, replacing w y,y→y = y − y by ȳ − ȳ proves the theorem.

Figure 2 .
Figure 2. The SR graph of the reaction network (31).
n and ∆ denotes the Laplacian in (27a).The directional derivative normal to the boundary ∂Ω is denoted by ∂/∂ν in (27c).The set Ω ⊆ R n is bounded, open and connected with twice continuously differentiable boundary ∂Ω.