Simplifying the Jacobian Criterion for precluding multistationarity in chemical reaction networks

Simplifying the Jacobian Criterion for precluding multistationarity in chemical reaction networks

February 14, 2012
Abstract

Chemical reaction networks taken with mass-action kinetics are dynamical systems that arise in chemical engineering and systems biology. In general, determining whether a chemical reaction network admits multiple steady states is difficult, as this requires determining existence of multiple positive solutions to a large system of polynomials with unknown coefficients. However, in certain cases, various easy criteria can be applied. One such test is the Jacobian Criterion, due to Craciun and Feinberg, which gives sufficient conditions for ruling out the possibility of multiple steady states. A chemical reaction network is said to pass the Jacobian Criterion if all terms in the determinant expansion of its parametrized Jacobian matrix have the same sign. In this article, we present a procedure which simplifies the application of the Jacobian Criterion, and as a result, we identify a new class of networks for which multiple steady states is precluded: those in which all chemical species have total molecularity of at most two. The total molecularity of a species refers to the sum of all of its stoichiometric coefficients in the network. We illustrate our results by examining enzyme catalysis networks.
Keywords: chemical reaction networks, mass-action kinetics, multiple steady states, Jacobian Criterion, injectivity, species-reaction graph

1 Introduction

A test for answering the question of multistationarity, which is implemented in some of the above-mentioned software programs, is the Jacobian Criterion, which is the focus of this article. The criterion, which is due to Craciun and Feinberg, gives sufficient conditions for precluding multiple steady states [5, 6, 7, 8]. This work has been extended by Banaji and Craciun [1, 2] and by Feliu and Wiuf [16]. A chemical reaction network is said to pass the Jacobian Criterion (or to be ‘injective’) if all terms in the determinant expansion of its parametrized Jacobian matrix have the same sign (see Lemma 3.4). We prove that for a given network the Jacobian Criterion can be checked by examining certain subnetworks (Theorem 4.9).

Our main result (Algorithm 6.1) is a procedure which simplifies the application of the criterion in two steps: first, the network is reduced to a simpler network, and then the number of subnetworks that must be examined is decreased. As an important consequence, we identify a new class of networks for which multiple steady states is precluded: those networks in which each chemical species has total molecularity of less than or equal to two (Theorem 5.1). The total molecularity of a species refers to how many non-flow reactions it appears in, where reversible reactions are counted only once and each occurrence of the species is counted with its stoichiometric coefficient; see Definition 2.2. More generally, Algorithm 6.1 allows us in certain cases to rule out the possibility of multiple steady states simply by inspecting the network. In addition, the reduction in the number of subnetworks that must be examined lends support to the observation that typically very few terms (if any) in the Jacobian determinant have anomalous sign [5]; this topic is explored by Helton, Klep, and Gomez [18] and Craciun, Helton, and Williams [9]. A proposed application of our work is to the problem of enumerating networks that admit multiple steady states: our algorithm quickly removes from consideration many small networks that are computationally expensive to enumerate (see Remark 6.4). Finally, we provide examples which demonstrate that our approach is complementary to that of the Species-Reaction Graph Theorem of Craciun and Feinberg [7, 10].

We note that our results are stated for continuous-flow stirred-tank reactors (CFSTRs), chemical reaction networks in which all chemical species are removed at rates proportional to their concentrations (see Definition 2.2). However, even for more general networks in which not all species take part in the outflows, our results often nevertheless apply in one of the following two ways. First, one can augment a network by any missing outflows so that it becomes a CFSTR, and then apply our results to determine whether it passes the Jacobian Criterion. Then, under mild assumptions on the original network (such as weak-reversibility; see Definition 2.2), the recent work of Craciun and Feinberg allows us to conclude that the original network also precludes multiple steady states [8, Corollary 8.3]. A second approach is due to work of Feliu and Wiuf. Their extended definition of the Jacobian Criterion applies to all networks [16], thereby bypassing the need to augment a network by missing outflows. Finally, we emphasize that the Jacobian Criterion applies for more general dynamics than mass-action kinetics [1, 9]; our results also may be used in those settings.

This article is organized as follows. Section 2 provides an introduction to chemical reaction systems. Section 3 describes the Jacobian Criterion, and Section 4 lays the groundwork for simplifying its implementation. The main result of Section 5 (Theorem 5.1) precludes multistationarity from any network in which each chemical species has total molecularity at most two. Section 6 presents Algorithm 6.1, our procedure which simplifies the application of the Jacobian Criterion. Finally, Section 7 contains examples such as enzyme catalysis networks which illustrate our main results.

2 Chemical reaction network theory

In this section we review the standard notation and define ‘total molecularity.’ An example of a chemical reaction is denoted by the following:

 2X1+X2 → X3 .

The are called chemical species, and and are called chemical complexes. Assigning the reactant complex to the vector and the product complex to the vector , we rewrite the reaction as . In general we let denote the total number of species , and we consider a set of reactions, each denoted by

 yk→y′k ,

for , and , with . We write each such vector as , and we call the stoichiometric coefficient of species in complex . For ease of notation, when there is no need for enumeration we will drop the subscript from the notation for the complexes and reactions.

A self-catalyzing reaction will refer to a reaction, such as , in which there is a species (here, ) that is a self-catalyst: it appears with a higher stoichiometric coefficient in the product complex than in the reactant complex.

Definition 2.1.

Let , and denote finite sets of species, complexes, and reactions, respectively. The triple is called a chemical reaction network if it satisfies the following:

1. for each complex , there exists a reaction in for which is the reactant complex or is the product complex, and

2. for each species , there exists a complex that contains .

We say that a species is a reactant species (respectively, product species) if it appears in a reactant (respectively, product) complex of some reaction. A subset of the reactions defines the subnetwork , where denotes the set of complexes that appear in the reactions , and denotes the set of species that appear in those complexes. A network decouples if there exist nonempty sets of reactions and that contain disjoint sets of species and such that . Otherwise, we say that the network is coupled.

Definition 2.2.
1. A chemical reaction network is weakly-reversible if for each reaction , there exists a sequence of reactions in from to (that is, reactions ). A reaction is reversible if its reverse reaction is also a reaction in .

2. A flow reaction contains only one molecule; such a reaction is either an inflow reaction or an outflow reaction . A non-flow reaction is any reaction that is not a flow reaction.

3. A chemical reaction network is a continuous-flow stirred-tank reactor network (CFSTR network) if it contains all outflow reactions (for all ).

4. Consider a chemical reaction network whose set of non-flow reactions is given by:

 y1 →y′1 y2 →y′2 … yl →y′l yl+1 ⇆y′l+1 yl+2 ⇆y′l+2 … yl+k ⇆y′l+k ,

where none of the first reactions is reversible. The total molecularity of a species in the network, denoted by , is the following integer:

 TM(Xi)=l+k∑j=1(yji+y′ji) ,

where is the stoichiometric coefficient of species in the complex .

Note that Craciun and Feinberg use the term “feed reactions” for inflow reactions and “true reactions” for non-flow reactions. In chemical engineering, a CFSTR refers to a tank in which reactions occur. An inflow reaction represents the flow of species (at a constant rate) into the tank in which the non-flow reactions take place, and an outflow reaction represent the removal of a species (at rate proportional to its concentration).

2.1 Dynamics

The concentration vector will track the concentration of the species at time . For a chemical reaction network , mass-action kinetics defines the following system of ordinary differential equations:

 ˙x(t)=|R|∑k=1κkx(t)yk(y′k−yk)=:f(x(t)) , (1)

where the last equality is a definition. Here, we make use of multi-index notation: . (By convention, ). Also, the ’s denote unknown (positive) reaction rate constants. A concentration vector is a (positive) steady state of the mass-action system (1) if .

Remark 2.3.

We can scale each concentration so that each outflow reaction (a reaction with and ) has outflow rate equal to one ( for such a reaction) [5, Remark 2.9]. We will make this assumption beginning in the next example.

In the following example and all others in this article, we will label species by distinct letters such as rather than .

Example 2.4.

Consider the following CFSTR:

 0[k2=1]k1⇄A0[k4=1]k3⇄BA+Bk5→2A

The total molecularity of species and are three and one, respectively. The mass-action differential equations (1) for this network are the following:

 ddtxA =k1−xA+k5xAxB ddtxB =k3−xB−k5xAxB .
Definition 2.5.

A CFSTR is said to admit multiple steady states or is multistationary if there exist reaction rate constants whose resulting system (1) admits two or more positive steady states.

2.2 Square networks

Throughout this article, we will be interested in subnetworks in which the number of species and the number of reactions are equal. We will call such subnetworks ‘square.’

Definition 2.6.
1. A chemical reaction network is said to be -square if it contains exactly species and exactly directed reactions. A network is said to be square if it is -square for some . (A -square network is the empty network.)

2. For an -square chemical reaction network with reaction set , the reactant matrix and reaction matrix are the integer -matrices

3. For a square chemical reaction network , we define its orientation as follows:

 Or(N):=sign(det(MN)det(RN)) .

As defined, our reactant and reaction matrices are the transposes of those of Craciun and Feinberg [5]. However, the orientation of a network is unaffected by transposing either matrix. Note that a reaction is encoded by the vector in the reaction matrix , which is the negative of its reaction vector which appears in the mass-action differential equations (1).

3 The Jacobian Criterion

The Jacobian Criterion concerns the determinant expansion of the Jacobian matrix of the differential equations (1) of a chemical reaction network.

Definition 3.1.

Consider a chemical reaction network with species. The Jacobian matrix associated to is the -matrix whose -th entry is the partial derivative , where is the -th differential equation in (1). The Jacobian determinant of , denoted by , refers to the determinant of the negative of the Jacobian matrix, which is a polynomial in the variables (for ) and .

Note that Craciun and Feinberg use the notation for [5]; here denotes the negative of the mass-action differential equations for the network obtained from by removing all inflow reactions.

Example 3.2.

We now return to the CFSTR examined in Example 2.4, which we denote by . Its Jacobian matrix is:

 (−1+k5xBk5xA−k5xB−1−k5xA) .

Therefore, the Jacobian determinant of is the following polynomial:

 (2)

The following result, which is due to Craciun and Feinberg [5], states that the determinant expansion of the Jacobian matrix has at most one term in the case of a square chemical reaction network.

Lemma 3.3 (Theorem 3.2 of [5]).

Let be a chemical reaction network that contains species and directed reactions. Then the Jacobian determinant of , , is either zero or a monomial (in the variables and ) whose coefficient is given by the following product:

 det(MN)⋅det(RN) ,

where is the reactant matrix and is the reaction matrix.

The following result is due to Craciun and Feinberg.

Lemma 3.4 (Theorems 3.1 and 3.3 [5]).

Let be a chemical reaction network that contains species, non-flow reactions, outflow reactions, and inflow reactions. Let denote the Jacobian determinant associated with , which is a polynomial in variables . Let enumerate the subnetworks of that contain exactly reactions, such that the first do not contain inflow reactions.

1. Each term in the expansion of is the Jacobian determinant of a unique subnetwork . Conversely, the Jacobian determinant of each subnetwork is either zero or is a monomial in the Jacobian determinant of . In other words,

 Jac(G)=(r+ro+ris)∑i=1Jac(Gi)=(r+ros)∑i=1Jac(Gi) .
2. Moreover, if is a CFSTR, then the following are equivalent:

• The polynomial is positive for all and .

• Each term is either zero or a positive monomial.

• Each subnetwork has non-negative orientation: .

If these equivalent conditions hold, then does not admit multiple steady states.

One of the main results in the next section is Theorem 4.9, which will extend the equivalent conditions that appear in part 2 of Lemma 3.4 so that instead of checking whether the subnetworks have non-negative orientation, we need only check the orientations of certain smaller ‘embedded’ networks (Definition 4.6). If any of these conditions hold for a CFSTR , for instance if has only positive terms, then we say that passes the Jacobian Criterion (or is ‘injective’).

Remark 3.5.

always has at least one positive term when contains all outflow reactions (such as for a CFSTR). This arises because one of the subnetworks consists of the outflow reactions. Clearly, the reactant and reaction matrices of this subnetwork are both the identity matrix, so its orientation is indeed positive.

The aim of this article is to simplify the computation that must be performed in determining whether a CFSTR passes the Jacobian Criterion, that is, whether all other terms in the determinant expansion are positive. In fact, there are typically few negative terms [5]; this topic is explored in [9, 18].

Notation 3.6.

Although Feliu and Wiuf have extended the definition of the Jacobian Criterion to apply also to non-CFSTR networks [16], the main focus of this article is the case when the chemical reaction network is a CFSTR. The subnetwork of containing all the non-flow reactions but none of the flow reactions will be called the non-flow subnetwork of and denoted by . Conversely, any reaction network is contained in a CFSTR obtained by including all outflow reactions: we denote this CFSTR by

 ~N:={\SS, C∪\SS∪{0}, R∪{Xi→0}Xi∈\SS} .
Example 3.7.

We now return to the CFSTR examined in Examples 2.4 and 3.2, which has species and non-flow reaction. Its non-flow subnetwork is:

 G:={A+B→2A} .

The -square subnetworks of that do not contain inflow reactions are the following:

 G1 ={A→0, B→0} , G2 ={A→0, A+B→2A},and G3 ={B→0, A+B→2A} .

As noted earlier, the subnetwork comprised of all outflow reactions, , has positive orientation: . Now we examine the reactant and reaction matrices of the remaining subnetworks and :

 MG2 =(1011) RG2 =(10−11) (3) MG3 =(0111) RG3 =(01−11) .

We see that the orientations of the two networks are and . Therefore, Lemma 3.4 implies that the CFSTR fails the Jacobian Criterion. This is consistent with the computation of the Jacobian determinant (2) in Example 3.2, which we saw contains one negative term.

3.1 Reactant and reaction rank

The Jacobian Criterion requires that each -square subnetwork of has either a singular reactant matrix, a singular reaction matrix, or both matrices have nonzero determinants with the same sign. We now give a name to these cases:

Definition 3.8.

A square chemical reaction network is said to have full reactant-rank if , and is said to have full reaction-rank if .

If a square network has both full reactant-rank and full reaction-rank, then it clearly has nonzero orientation.

Definition 3.9.

In a square network , we say that species is changed in reaction if .

In other words, species is changed in reaction if the concentration of species is changed when reaction takes place. For example, species is changed in the reaction and in the reaction . On the other hand, species and are changed only in the first reaction.

Lemma 3.10.

If is a square network with full reactant-rank, then:

1. Each species of is a reactant species.

2. No reaction of is an inflow reaction or a generalized inflow reaction of the form (where ).

3. No two reactions of share the same reactant complex.

If is an -square network with full reaction-rank, then:

1. Each species of is changed in some reaction.

2. does not contain both a reaction and its reverse.

Proof.

This result follows because nonsingular square matrices can not have a zero column, a zero row, two identical rows, or a row and its negative. ∎

3.2 Known tests for passing the Jacobian Criterion

Known sufficient conditions on the reaction network for passing the Jacobian Criterion arise from the following results: the SCL Graph Theorem of Schlosser and Feinberg [25], the Species-Reaction Graph Theorem of Craciun and Feinberg [7, 8] and the related work of Banaji and Craciun [1, 2], and results of Helton, Klep, and Gomez [18]. In this article, we present a complementary test (Algorithm 6.1).

4 Analysis of subnetworks

Theorems 4.9 and 4.11 in this subsection form the basis of the procedure (Algorithm 6.1) that we will introduce to simplify the implementation of the Jacobian Criterion. The following example motivates our results.

Example 4.1.

We return to the network in Example 3.7. Note that in both the reactant and product matrices (3) of the subnetwork , the first row is the vector . So, the orientation of depends only on the entries in the lower-right of the two matrices. In other words, by letting denote the 1-square network obtained from by ‘removing’ the outflow reaction and ‘removing’ the species from the non-flow reaction :

 N2={B→0} ,

it follows that . Similarly, for the network obtained from by removing and the species from the non-flow reaction:

 N3={A→2A} ,

we have .

In this section, we make the simplifications performed in Example 4.1 formal: the first rows in the matrices and (3) will be called ‘pivotal rows’ (Definition 4.2), Theorem 4.9 will permit us to restrict our attention to the square ‘embedded’ networks rather than the larger -square subnetworks , and Theorem 4.11 will allow us to remove the species from the network. In general the networks will be smaller and fewer in number, thereby simplifying the Jacobian Criterion.

4.1 Analysis of row and column pivots

Definition 4.2.
1. For any matrix , if there exists a row with one positive entry, (respectively, ) for some , and all other entries zero, for , we say that is a positive (respectively, negative) row pivot and the -th row is a pivotal row.

2. For any matrix , if (respectively, ) for some and for all , we say that is a positive (respectively, negative) column pivot and the -th column is a pivotal column.

Notation 4.3.

For an -matrix and indices , we define the associated matrix which is obtained from by removing the -th row and the -th column.

In the following theorem, we translate pivot conditions for reactant and reaction matrices into the language of reaction networks. Note that reactant matrices are non-negative matrices, so their pivots are necessarily positive pivots.

Proposition 4.4.

Let be a square chemical reaction network with reactant matrix and reaction matrix . Then:

1. is a row pivot if and only if the -th reactant contains only species (in other words, the -th reaction is ).

2. is a column pivot if and only if species appears in the reactant complex of the -th reaction and in no other reactant complex.

3. is a row pivot if and only if only species is changed in reaction .

4. is a column pivot if and only if species is changed only in reaction .

5. Both and are row pivots if and only if the -th reaction is of the form

 ajXj→bjXj         where   0≠aj≠bj

Moreover, is a positive pivot if and a negative pivot if . If and are positive row pivots then we will say that is a pivotal reaction.

6. Both and are column pivots if and only if species appears only in reaction . Moreover, is a negative column pivot if and only if is a self-catalyst in reaction .

7. Suppose that for some , both and are positive pivots. Then

 Or(N)=sign(det(M)det(R))=sign(det(Mˆij)det(Rˆij)) .

If is a (positive) pivot and is a negative pivot, then

 Or(N)=sign(det(M)det(R))=−sign(det(Mˆij)det(Rˆij)) .
Proof.

This proposition follows easily from definitions and elementary linear algebra. ∎

4.2 Analysis of embedded networks

Theorem 4.9 in this section will tell us that instead of analyzing all -square subnetworks of a network , we can remove the outflow reactions from and then remove the outflow species from the remaining reactions, and analyze those simpler ‘embedded networks.’ First we must define what is meant by removing a reaction or a species.

Definition 4.5.

Let be a chemical reaction network.

1. Consider a subset of the species , a subset of the complexes , and a subset of the reactions .

• The restriction of to , denoted by , is the multiset of reactions obtained by taking the reactions in and removing all species not in .

• The restriction of to , denoted by , is the set of (reactant and product) complexes of the reactions in .

• The restriction of to , denoted by , is the set of species that are in the complexes in .

2. The network obtained from by removing a subset of species is the network

 {\SS∖{Xi}, C|R|\SS∖{Xi}, R|\SS∖{Xi}} .
3. The network obtained from by removing a set of reactions is the subnetwork

 {\SS|C|R∖{y→y′}, C|R∖{y→y′}, R∖{y→y′}} .

Note that is viewed as a multiset (a generalized set in which an element may appear more than once), so that in the following definition, embedded networks are indeed square.

Definition 4.6.

Let be a chemical reaction network with species, and let be its non-flow subnetwork.  For , a -square embedded network of (or of ), which is defined by two sets, a -subset of the species, , and a -subset of the non-flow reactions, , that involve all species of , is the -square network consisting of the reactions .

The -square embedded networks with negative orientation are precisely the reactions of the form with . Therefore, a network contains a negatively oriented -square embedded reaction if and only if it contains a self-catalyzing reaction.

The following proposition will be used in the proofs of Theorems 4.9 and 4.11.

Proposition 4.7.

Let be a square network that contains one of the following:

1. an outflow reaction or a generalized outflow reaction (where ), or

2. a species that appears in exactly one reaction, and the reaction is not self-catalyzing.

Then , or the network obtained from by removing the aforementioned reaction and species is a square network and .

Proof.

Part 1 follows easily from Proposition 4.4, as such a reaction is a pivotal reaction. Now, part 2 is the case when some species appears in exactly one reaction in . If is in the product complex , then by Lemma 3.10. In the remaining case ( is in the reactant complex ), then is a pivotal reaction, so . ∎

The following lemma generalizes what we saw in Example 4.1 in the CFSTR setting: each -square subnetwork of defines an embedded network of , and they have the same orientation.

Lemma 4.8.

Let be a chemical reaction network with species, and let denote the subnetwork of non-flow reactions. Consider the map from -square subnetworks of that do not contain inflow reactions to square embedded networks of , which maps to the network obtained from by first removing all outflow reactions (for ) and then removing those outflow species from the remaining non-flow reactions of . Then this map is surjective, and

 Or(Gi)=Or(Ni) . (4)

(If is empty, its orientation is defined to be +1.)

Proof.

The surjection is clear: both objects are defined by a subset of and a subset of whose cardinalities sum to . Equation (4) follows from noting that the outflow reactions are pivotal reactions, and then applying Proposition 4.7. ∎

From Lemma 4.8, we obtain the following theorem, which will be used in Step III in Algorithm 6.1.

Theorem 4.9.

For a CFSTR , with its subnetwork of non-flow reactions, the following are equivalent:

1. passes the Jacobian Criterion.

2. Each square embedded network of has non-negative orientation.

3. has no self-catalyzing reaction, and each square embedded network of which satisfies the following properties has non-negative orientation:

• has no outflows, inflows, generalized inflow reactions , or generalized outflows (where ).

• Each species appears in at least two reactions (where a pair of reversible reactions is counted only once).

Proof.

Parts 1 and 2 are equivalent by Lemma 3.4 and Lemma 4.8. Clearly, part 2 implies part 3, so it remains only to show the converse. It suffices to show that if is a square embedded network which violates one of the conditions of part 3, then either it has zero orientation or we can find a ‘smaller’ embedded network which satisfies part 3 and has the same orientation as . In the case of inflows and generalized inflows, the orientation of is zero. All remaining cases are handled by Proposition 4.7. ∎

The following corollary follows easily from Theorem 4.9.

Corollary 4.10.

If a CFSTR contains a self-catalyzing reaction, then it fails the Jacobian Criterion. Conversely, if a CFSTR does not contain a self-catalyzing reaction, then the only square embedded networks in Theorem 4.9 that must be checked are those of size at least two.

Corollary 4.10 allows us to see quickly by inspection that the network introduced in Example 2.4 (and analyzed in Examples 3.2, 3.7, and 4.1) fails the Jacobian Criterion.

4.3 A network reduction theorem

This subsection concerns those reactions or species that can be removed from a CFSTR without affecting whether the network passes the Jacobian Criterion. The following theorem will be the foundation of Step II of Algorithm 6.1.

Theorem 4.11.

Let be a CFSTR, and let be a network obtained from by performing any of the following operations:

1. Remove any inflow or outflow reactions.

2. Remove any generalized inflow reactions .

3. Remove reactions containing only one species for which .

4. Remove any species that is not a reactant species.

5. Remove any species that appears in exactly one reaction (or a pair of reversible reactions), unless the species is a self-catalyst of the reaction.

Then passes the Jacobian Criterion if and only if all square embedded networks of have nonnegative orientation.

Proof.

Let (respectively, ) denote the subnetwork of non-flow reactions of (respectively, ). By Theorem 4.9, we need only show that all square embedded networks of have non-negative orientation if and only if all square embedded networks of do. First, we consider the operations through , in which is obtained from by removing a reaction . In this case,

 {SENs of G} = { SENs of Gred} ∪ { SENs of G that involve y→y′} , (5)

where ‘SEN’ is shorthand for ‘square embedded network.’ So, we need only show that if is a square embedded network of , then either has zero orientation, or there exists a square embedded network of which has the same orientation as . In the case that the reaction is an inflow or a generalized inflow, . So, the remaining cases are when is an outflow reaction or is of the form (where ); this is handled by Proposition 4.7.

Now, only parts 4 and 5 remain to be proven. Here, is obtained from by removing a species . The analogue to (5) now is that the set of square embedded networks of are precisely the square embedded networks of that do not involve species . Again, we must show that a square embedded network of either has zero orientation, or there exists a square embedded network of with the same orientation as . If is not a reactant species of (in the case of part 4), then has zero orientation by Lemma 3.10. Finally, if satisfies the condition in part 5, then appears in exactly one reaction in . This case is completed by applying Proposition 4.7. This completes the proof. ∎

5 Networks with total molecularity at most two do not admit multiple steady states

In this section, we prove the following theorem, which concerns networks in which the maximum total molecularity (Definition 2.2) is at most two.

Theorem 5.1.

If all species of a CFSTR have total molecularity less than or equal to two, then it passes the Jacobian Criterion, and therefore does not admit multiple steady states.

Theorem 5.1 has the following corollary.

Corollary 5.2.

For any weakly-reversible chemical reaction network (not necessarily a CFSTR) in which all species have total molecularity less than or equal to two, the network does not admit multiple steady states.

Remark 5.3.

Theorem 5.1 cannot be improved to higher total molecularity. For example, a CFSTR that contains the reaction fails the Jacobian Criterion (by Corollary 4.10).

Remark 5.4.

Theorem 5.1 can be proven by applying the Species-Reaction Graph Theorem of Craciun and Feinberg [7, Theorem 1.1]. Rather than introducing the terminology and results of species-reaction graphs, we give a direct proof here.

Our proof of Theorem 5.1 will proceed by showing that for a chemical reaction network that has total molecularity at most two, if a square embedded network has nonzero orientation, then it in fact has positive orientation. The proof requires the following lemmata, which will be used to reduce to studying square embedded networks, and then to elucidate what values can be taken by the determinants of the reactant and reaction matrices of such networks.

Lemma 5.5.

Let be a chemical reaction network whose species all have total molecularity less than or equal to some . Then the same holds for any embedded network of .

The proof of Lemma 5.5 follows immediately from the definitions of total molecularity and embedded networks.

Lemma 5.6.

Let be a coupled -square network with at least two reactions (none of which are inflows or outflows) such that all species have total molecularity at most two. Assume that has nonzero orientation. Then each reaction of has the form of or .

Proof.

We first claim that if is a coupled, -square network in which all species have total molecularity two, then each reaction of must contain exactly two molecules. This is because no reaction can have only one molecule (because inflows and outflows are not allowed). So, each of the reactions must have exactly two molecules, and it follows that the total molecularity of each species must be two. Now, a reaction of the form would yield zero orientation for the network (because the corresponding row of the reactant matrix would be zero), and a reaction is not allowed, because then the network would decouple into this reaction and the remaining subnetwork that does not contain . So, only two types of reactions are possible: and . ∎

Lemma 5.7.

Consider an -matrix such that

1. all diagonal entries are ,

2. the entries above the diagonal (mod ), , are or , and

3. all other entries are 0.

Then . (In fact, the determinant is , , or .)

Proof.

The determinant is . ∎

We now can prove Theorem 5.1 and Corollary 5.2.

Proof of Theorem 5.1.

By Theorem 4.9, Lemma 4.8, and Lemma 5.5, we need only show that if is a square network with nonzero orientation and with maximum total molecularity at most two, then . Let denote the number of species of . Now, if decouples, then we could reduce to analyzing the two subnetworks separately both of which have fewer than species; the network passes the Jacobian Criterion if and only if the two subnetworks do. The case is trivial (only is allowed, and in this case both reaction and reactant matrices are the same, so their determinants are the same), so we may proceed by induction on in the case of a coupled network .

Now Lemma 5.6 applies, so we know that the reactions of are of the form and . By relabelling the species, we can ensure that the first reaction contains species and , the second reaction contains and , and so on. Therefore, both the reactant matrix and the reaction matrix now have the form of the matrix in Lemma 5.7. By the lemma, both determinants are non-negative, so the orientation of the network is either zero or positive.

Finally, Lemma 3.4 implies that a CFSTR network that passes the Jacobian Criterion does not admit multiple steady states. ∎

Proof of Corollary 5.2.

For a weakly-reversible network that has total molecularity at most two, we consider the CFSTR obtained by augmenting the network by all outflows. Theorem 5.1 implies that passes the Jacobian Criterion. Finally, results due to Craciun and Feinberg allow us to conclude that the such a network precludes multiple steady states [8, Theorem 8.2 and Corollary 8.3]. ∎

6 Procedure for simplifying the Jacobian Criterion

Algorithm 6.1 is the main result of this article. It combines the results in previous sections into a procedure that simplifies the implementation of the Jacobian Criterion. After pre-processing in Step I, the algorithm reduces a chemical reaction network to a simpler network in Step II, and then decreases the number of embedded networks of the reduced network that must be examined in Step III.

Algorithm 6.1 (Simplifying the Jacobian Criterion).

Step I (Pre-processing step).
Input:
a CFSTR .
Output: either whether passes the Jacobian Criterion, or (to be passed to Step II).

1. If has a self-catalyzing reaction (an embedded reaction with ), then automatically fails the Jacobian Criterion.

2. If each species of has total molecularity at most two, then automatically passes the Jacobian Criterion.

If neither applies, then pass to Step II.
Step II (Reduction step).
Input: a CFSTR that passes Step I.
Output: a non-flow chemical reaction network such that the CFSTR passes the Jacobian Criterion if and only if does.
Perform the following operations on the network , in any order, as long as any one can be completed.

1. Remove any inflow or outflow reaction.

2. Remove any generalized inflow reaction .

3. Remove any reaction that contains only one species.

4. Remove any species that is not a reactant species.

5. Remove any species that appears in only one reaction (or in a pair of reversible reactions and in no other reactions). In particular, remove any species whose total molecularity is one.

6. Remove one copy of any repeated reactions.

If the reduced network is the empty network or all species in the reduced network have total molecularity at most two, then passes the Jacobian Criterion. Otherwise, pass the reduced network to Step III.
Step III (Analysis of embedded networks step).
Input: a chemical reaction network (such as the output of Step II) which passes Step I.
Output: whether passes the Jacobian Criterion (if is the output of Step II, then this step will determine whether the input of Step II passes the Jacobian Criterion).
Consider all coupled square embedded networks of with the following properties:

1. contains none of the following reactions:

• inflow or outflow reactions,

• generalized inflow reactions ,

• reactions containing only one species,

• both a reaction and its reverse reaction,

• two reactions that share the same reactant complex.

2. has at least two species.

3. Each species of is a reactant species.

4. Each species appears in at least two reactions (where pairs of reversible reactions are counted only once).

5. At least one species has total molecularity three or more.

6. has nonzero orientation.

Output: The network passes the Jacobian Criterion if and only if all such square embedded networks have nonnegative orientation.

Proof of correctness of Algorithm 6.1.

Part 1 of Step I follows from Corollary 4.10, and part 2 follows from Theorem 5.1. Step II is valid by Theorem 4.11. Parts through of Step III follow from Lemma 3.10, Corollary 4.10, and Theorem 4.9. Part 5 of the Step III follows from Theorem 5.1, and part 6 simply avoids embedded networks with zero orientation. ∎

Remark 6.2.

Instead of applying Step III to a reduced network that is the output of Step II, one can apply other tests to check the Jacobian Criterion. For instance, one can apply the Species-Reaction Graph Criterion [7] to the reduced network or directly compute the expansion of the Jacobian determinant using available software [13, 21, 22]. The reduced network typically contains fewer species and fewer reactions than the original network, so computer computation will be more tractable.

Remark 6.3.

In practice, the number of square embedded networks that must be examined in Step III of Algorithm 6.1 is relatively small compared to the number of -square subnetworks in Lemma 3.4. Comparisons of the number of square subnetworks versus the number of square embedded networks to be examined appear in the last two columns of Table 1 and in the comments of Table 2. In addition, the determinants of such square embedded networks (as given in Lemma 3.3) that must be computed are of smaller matrices than those of .

Remark 6.4.

We propose that Algorithm 6.1 could aid in enumerating chemical reaction networks that admit multiple steady states. Recall that failing the Jacobian Criterion is a necessary condition for admitting multiple steady states (Lemma 3.4), so let us consider the problem of identifying the set of networks that fail the Jacobian Criterion. More specifically, we now define a bimolecular network to be a network in which all complexes contain at most two molecules, and then pose the following question: which bimolecular CFSTRs fail the Jacobian Criterion?

Accordingly, we enumerated all bimolecular CFSTRs that contain two pairs of reversible non-flow reactions and we found that there are 386 such networks. Of these 386 networks, 142 have a maximum total molecularity of two or fewer, and thus immediately pass the Jacobian Criterion. It is interesting to note that for these 142 networks, several distinct species typically have the same total molecularity. Networks that have such symmetry are computationally the most expensive to enumerate [11]. As we are interested in listing only those networks that fail the Jacobian Criterion, we can avoid these computationally expensive networks.

Moreover, of the remaining 244 networks, applying Step II of the Algorithm 6.1 immediately reduces the number of networks that need to be examined. In other words, many networks are related in whether or not they pass the Jacobian Criterion and thus studying all 244 networks is redundant. When we discard networks in which some species has total molecularity equal to one, we are left with only 55 networks whose multistationarity needs to be investigated. Details of our enumeration of bimolecular CFSTRs appear in [20]. Note that a similar enumeration of small bimolecular reaction networks was undertaken by Deckard, Bergmann, and Sauro [11], from which Pantea and Craciun sampled networks to compute the fraction of such networks that pass the Jacobian Criterion [23, Figure 1].

7 Examples

The examples in this section illustrate that Algorithm 6.1 in some cases can determine whether a network passes the Jacobian Criterion more quickly than applying the Species-Reaction Graph Criterion of Craciun and Feinberg [7], Deficiency Theory [14, 15], or the results of Banaji and Craciun [2]. Moreover, we are able to analyze examples in [5, 7, 10, 25], which appear in Tables 1 and 2. In addition, we demonstrate that instead of checking the orientations of all the -square networks from Lemma 3.4 (which entails computing determinants of many -matrices), we need only compute the orientations of certain square embedded networks, which are both smaller and fewer in number. The final columns in both Tables 1 and 2 illustrate the great reduction in the number of subnetworks that must be examined. We begin with examples in which no orientations must be computed, as checking the Jacobian Criterion sometimes can be performed by inspection.

Example 7.1 (Networks containing a self-catalyzing reaction).

Recall that Corollary 4.10 states that if a chemical reaction network contains a self-catalyzing reaction, then it fails the Jacobian Criterion. One such network is the one examined in Examples 3.7 and 4.1. Other networks that are resolved easily by Step I of Algorithm 6.1 include networks (vi), (vii), and (viii) in Table 1.

In addition to the self-catalyzing networks, the other networks that are resolved in Step I of Algorithm 6.1 are those with total molecularity at most two.

Example 7.2 (Networks with total molecularity at most two).

Consider the following network examined by Gnacadja et al. [17]:

 L+R ↔LRA+R↔AR (6) L+T ↔LTA+T↔AT .

This ligand-receptor-antagonist-trap network summarizes the four pairs of binding and dissociating reactions that occur in the presence of two drugs, a ‘receptor antagonist’ and a ‘trap,’ which are used to treat patients with rheumatoid arthritis. In their work, Gnacadja et al. apply Deficiency Theory to determine the existence and uniqueness of a steady state [17, Proposition 2]. If one is interested only in precluding multiple steady states, then the approach advocated here can be performed more quickly. It is easily seen that the total molecularity of each of the eight species, and , is no more than two. Therefore, Theorem 5.1 implies that a CFSTR version of this network passes the Jacobian Criterion, and does not admit multiple steady states. Moreover, as the underlying network (6) is weakly-reversible, we can conclude that multiple steady states is also ruled out in the non-CFSTR setting [8].

Other networks with total molecularity at most two include the network (iv) in Table 1 and the enzyme catalysis networks (i) and (iv) in Table 2. All of these networks can be seen by inspection to pass the Jacobian Criterion, without appealing to the SR Graph Criterion.

The networks in Example 7.2 had total molecularity at most two, and so could be seen immediately to pass the Jacobian Criterion. The networks in the next example do not have this property until after applying Step II of Algorithm 6.1.

Example 7.3 (Networks with total molecularity at most two after reduction).

Consider the CFSTR introduced by Banaji and Craciun whose non-flow subnetwork is:

 D⇆A+B+CE⇆A+B+CF⇆A+B (7)

Banaji and Craciun showed that the CFSTR fails the Species-Reaction Graph Criterion, but nevertheless passes the Jacobian Criterion, because the stoichiometric matrix is ‘strongly sign determined’  [2, §6.2]. We now show that Step II of Algorithm 6.1 can easily be performed on this network, to conclude that it passes the Jacobian Criterion.

In Step I of Algorithm 6.1, we note that the total molecularity of species and in network (7) are both three, and no self-catalyzing reactions exist. So, we begin Step II by removing species , , and , each of which has total molecularity equal to one. The resulting network is:

 0⇆A+B+C0⇆A+B+C0⇆A+B

Now the first two pairs of reversible reactions are the same, so we remove one copy. Additionally, we remove the remaining two generalized inflow reactions, to obtain:

 0←A+B+C0←A+B

Clearly, each species in this reduced network has total molecularity no more than two, so we conclude that the original CFSTR passes the Jacobian Criterion. Two other networks whose reduced networks can be obtained easily and have total molecularity at most two are the enzyme catalysis networks with inhibition (ii) and (iii) in Table 2.

While the above examples did not need to proceed beyond Step II of Algorithm 6.1, the networks in the following example do require this step. We will see that the algorithm greatly reduces the number of subnetworks that must be checked.

Example 7.4 (Networks requiring the full algorithm).

Consider the reversible network (i) in Table 1; its non-flow subnetwork is:

 A+B⇆PB+C⇆QC⇆2A (8)

The output of Step II of Algorithm 6.1 is the following network:

 A+B→0B+C→0C⇆2A

As only species has total molecularity greater than two, each square embedded network we consider must contain both and one of the reversible reactions (because we need not consider networks with both directions of a reversible reaction). However, such a -square embedded network would have a species ( or ) that appears in only one complex. Therefore, we need only consider the two -square embedded networks, which are the following:

 {A+B→0 B+C→0C→2A} ,and (9) {A+B→0 B+C→0C←2A} (10)

The first embedded network (9) has negative orientation, and the second network (10) has positive orientation, so a CFSTR version of the original network (8) fails the Jacobian Criterion. Note that our algorithm reduced the Jacobian Criterion from examining the -square subnetworks (the in Lemma 3.4) to two embedded networks. Networks (ii), (iii), and (v) in Table 1 can be analyzed similarly: they require examining only two, two, and three square embedded networks, respectively.

Our final example network requires more complicated analysis, but the number of square embedded ne