{'fraction_non_alphanumeric': 0.04789164697571314, 'fraction_numerical': 0.015395483666658548, 'mean_word_length': 3.8565006506565718, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Computational Approach to Extinction Events in Chemical Reaction Networks with Discrete State Spaces

Recent work of M.D. Johnston et al. has produced sufficient conditions on the structure of a chemical reaction network which guarantee that the corresponding discrete state space system exhibits an extinction event. The conditions consist of a series of systems of equalities and inequalities on the edges of a modified reaction network called a domination-expanded reaction network. In this paper, we present a computational implementation of these conditions written in Python and apply the program on examples drawn from the biochemical literature, including a model of polyamine metabolism in mammals and a model of the pentose phosphate pathway in Trypanosoma brucei. We also run the program on 458 models from the European Bioinformatics Institute's BioModels Database and report our results.


Introduction
A chemical reaction network describes the conversion of chemical reactants into products as a directed graph where the vertices are aggregates of species called complexes and the edges are reactions. A dynamical system can be associated to such a network in several ways. When the counts of reacting species are high, as is typical in industrial chemistry and pharmacology, the system can be reasonably modeled with a system of differential equations over the continuous state space of reactant concentrations. When the count of reacting species is low, however, as is typical in genetic and enzymatic systems in systems biology, it is more accurate to model the system over the discrete state space of molecular counts, for instance as a continuous-time Markov chain [3,4,18] or stochastic Petri Net [6,23].
Recent work has focused on when the long-term behaviors predicted by the continuous and discrete state models are markedly different. While classical work of T. Kurtz guarantees that the probability density function of the discrete model converges in an appropriate scaling limit to the solution of the continuous state model on compact time intervals [0, T ], T < ∞ [18], on the unbounded interval [0, ∞) this convergence may fail. For example, consider the following network, which was presented as Example 2.1 in the companion paper [13]: where the numbers correspond to the enumeration of the reactions. For almost all initial conditions and parameter values, the continuous state differential equation model predicts convergence to a strictly positive steady state. For the discrete state model, however, the reaction X 2 → X 1 may irreversibly deplete the species X 2 so that the inevitable final state of the system for M = X 1 (0) + X 2 (0) is {X 1 = M, X 2 = 0}.
Extinction events in discrete state space models were studied in the context of chemical reaction network theory by D.F. Anderson et al. in [2]. In that paper, the authors showed that a large subset of networks which exhibit "absolute concentration robustness" in the continuous state model [26] exhibit an extinction event when modeled with a discrete state space. A primary extension made in that paper was generalizing the notion of "differing in one species" to a domination relationship between complexes. This domination relation, and the corresponding relationship with extinction events in discrete state space systems, were further generalized and adopted to the conventions of Petri Net Theory by R. Brijder in [7]. Finally, in this paper's companion paper [13], M.D. Johnston et al. clarify this relationship by presenting sufficient conditions for an extinction event which depend upon the evaluation of computationally-tractable systems of equalities and inequalities.
The method presented in [13] involves creating a modified reaction network called a domination-expanded network. For example, for the network above, we may correspond the following expanded network: where the dashed edges are called "domination reactions." This network is then associated by certain rules with some systems of equalities and inequalities. The main result of [13] states that, if any of these systems cannot be satisfied then the discrete state space system must exhibit an extinction event.
For example, to the domination-expanded network above we associate the following system of equalities and inequalities on the vector α = ((α R ) 1 , (α R ) 2 , (α R ) 3 , (α D ) 1 , (α D ) 2 )) ∈ Z 5 ≥0 , where the entries correspond to the reaction edges: Since this system has no nontrivial solution, we may conclude that the discrete state space system has an extinction event.
In this paper, we implement the conditions of the companion paper [13] for affirming an extinction event in discrete state space systems into a computational package written in Python. The program utilizes a series of mixed-integer linear programming (MILP) modules for verifying the technical conditions required to generate the governing systems of equalities and inequalities, and also for evaluating these resulting systems. We then run the algorithm on 458 models from the European Bioinformatics Institute's BioModels Database and report our results. As an illustration of the power of the program, we further analyze a model of polyamine metabolism in mammals and a model of the pentose phosphate pathway in Trypanosoma brucei, which were identified by the BioModels database run as exhibiting an extinction event.

Background
In this section, we outlined the definitions, terminology, and background results necessary for our study of extinction events. We follow the notation of chemical reaction network theory (CRNT) [8]. As this background largely follows that of the companion paper [13], we will be brief and use the example presented in the introduction as a running example.

Chemical Reaction Networks
The following is the basic object of study in CRNT. y ij X j . The values y ij ∈ Z ≥0 are called stoichiometric coefficients. Allowing a slight abuse of notation, we will let y i denote both the ith complex and the corresponding complex vector y i = (y i1 , y i2 , . . . , y im ) ∈ Z m ≥0 . 3. The reaction set R = {R 1 , . . . , R r } consists of ordered pairs of complexes, e.g.
(y i , y j ) ∈ R. We furthermore define the mappings ρ, ρ ′ : {1, . . . , r} → {1, . . . , n} such that ρ(k) = i if y i is the source complex (tail of arrow) of the kth reaction, and ρ ′ (k) = j if y j is the product complex (head of arrow) of the kth reaction, so that we may write R k = (y ρ(k) , y ρ ′ (k) ). We also represent reactions with directed arrows, e.g. y i → y j or y ρ(k) → y ρ ′ (k) .
We will assume that (i) every species is contained in at least one complex; and (ii) every complex is contained in at least one reaction. The interpretation of reactions as directed edges naturally gives rise to a reaction graph G = (V, E) where the set of vertices is given by the complexes (i.e. V = C) and the set of edges is given by the reactions (i.e. E = R). A maximal set of connected complexes is called a linkage class (LC) while the maximal set of strongly connected complexes is called a strong linkage class (SLC). An SLC will be called terminal if there are no outgoing edges. A complex is terminal if it is contained in a terminal SLC. A subset of complexes Y ⊆ C is called an absorbing complex set if it contains every terminal complex and has no outgoing edges, i.e. y i → y j ∈ R and y i ∈ Y implies y j ∈ Y. Complexes y i ∈ Y are called Y-interior while complexes y i ∈ Y are called Y-exterior. Similarly, reactions y i → y j ∈ R where y i ∈ Y are called Y-interior and Y-exterior otherwise.
We associate the following structural matrices to a CRN. These matrices will be useful in both defining the dynamical system we are interested in, and in the computational implementation of our conditions for an extinction event.
is the matrix which has columns Y ·,i = y i , i ∈ {1, . . . , n}, where y i is the complex vector for the ith complex.

The adjacency matrix
and [I s ] i,j = 0 otherwise.

The unweighted Laplacian
Notice that Γ = Y · I a and A = I a · I s , where Γ is the stoichiometric matrix, and that I s may be obtained by replacing the negative entries in I a with ones, setting the remaining entries to zero, and taking the transpose. It therefore suffices to know Y (stoichiometry of CRN) and I a (connectivity of CRN) in order to generate all of the required structural data of a given CRN. It is typical in CRNT to use a Laplacian where the entries are weighted by the rate constant of the corresponding reaction [8,11]. By contrast, the unweighted Laplacian simply encodes the connectivity structure of the complexes in the reaction graph. Although we will not need to consider weighted Laplacians in this paper, we will use well-known properties of them to relate ker(A) to the terminal SLCs of the reaction graph [9].
To each reaction y i → y j ∈ R we associate a reaction vector y j − y i ∈ Z m which tracks the net gain and loss of each chemical species as a result of the occurrence of this reaction. The stoichiometric subspace is defined as Example 2.1. Reconsider the CRN given in the introduction, which was taken from Example 2.1 of [13]. We have the sets S = {X 1 , X 2 }, R = {X 1 + X 2 → 2X 2 , 2X 2 → X 1 + X 2 , X 2 → X 1 }, and C = {X 1 + X 2 , 2X 2 , X 2 , X 1 }. We have the following structural matrices: The stoichiometric subspace is given by S = span{(1, −1)} and the CRN is conservative with respect to the vector c = (1, 1). This represents the observation that X 1 + X 2 is constant.

Chemical Reaction Networks with Discrete State Spaces
The evolution of a CRN on the discrete state space Z m ≥0 is given by where X(t) = (X 1 (t), . . . , X m (t)) ∈ Z m ≥0 is the discrete state of the system at time t, and N(t) = (N 1 (t), . . . , N r (t)) and N k (t) ∈ Z ≥0 is the number of times the kth reaction has occurred up to time t. Several frameworks exist for precisely modeling the stochastic evolution of the discrete state X(t) over time, including the theories of continuous-time Markov chains (CTMC) [4,20] and stochastic Petri Nets [6,23]. We will not be interested in these details; rather, we will be interested in where in the discrete state space Z m ≥0 the trajectory X(t) may travel. For a similar treatment, see [22].
We adopt the following definitions from [13].
Definition 2.2. Consider a CRN on a discrete state space. Then: . . , m}.

A reaction
9. The CRN exhibits a guaranteed extinction event on Y ⊆ C if it has an extinction event on Y from every X ∈ Z m ≥0 .
Since subconservative CRNs on discrete state spaces have a finite number of states (see Theorem 1, [19]), the notion of recurrence presented above corresponds to the notion of positive recurrence from the theory of stochastic process [4,20].
Reconsider the CRN introduced in Example 2.1 and the initial state X = {X 1 = 0, X 2 = 3}. We have that the complexes X 2 and 2X 2 are charged at X but the complex X 1 + X 2 is not charged at X. The state X reacts to X ′ = {X 1 = 1, X 2 = 2} through either the reaction 2X 2 → X 1 + X 2 or X 2 → X 1 (i.e. X → X ′ ). We can furthermore see that the states . Since no reactions are charged at X ′′′ , however, we may not return to any of X, X ′ , or X ′′ . These states are therefore transient, as are all of the source complexes and reactions. It follows that the CRN has an extinction event on In fact, repeated application of X 2 → X 1 leads to an extinction event regardless of initial state, so that the CRN has a guaranteed extinction event on Y.

Domination-Expanded Reaction Networks
We briefly restate the key extensions to CRNs made in [13]. We start with the following, which was adapted from [2] and [7]. Definition 2.3. Consider two complexes of a CRN, y i , y j ∈ C, y i = y j . We say that y i dominates y j if y j ≤ y i . We define the domination set of a CRN to be The motivation underlying the domination relations (y i , y j ) ∈ D is that, if there is sufficient molecularity in the discrete state space CRN for a reaction from y i to occur from state X ∈ Z m ≥0 , then there is necessarily sufficient molecularity for a reaction from y j to occur from X. We use the domination relations to construct the following.
In the context of dom-CRNs, we will also represent the domination relations (y i , y j ) as y i → y j . Note that reaction goes from the dominating complex to the dominated complex, i.e. y j ≤ y i implies y i → y j is added to the CRN.
On a dom-CRN, we are interested in the following subgraphs.
Y-external forests define a flow from the distal portion of the reaction graph to the absorbing complex set Y in the dom-CRN. By convention, Y-exterior forests will include all Y-interior reactions of the dom-CRN. We will be interested in the following property on Y-exterior forests.
Otherwise, the Y-exterior forest is said to be unbalanced.
Condition 3 of Definition 2.6 may be interpreted as stating that, for Y-exterior complexes, the weight of the incoming edges may not be greater than the outgoing edge. Notice that Conditions 1 − 3 generated a system of equalities and inequalities on the edges of the dom-CRN. The computational implementation of the conditions of Definition 2.6 is the primary focus of this paper. We label y 1 = X 1 + X 2 , y 2 = 2X 2 , y 3 = X 2 , and y 4 = X 1 . We have the domination relations y 3 ≤ y 1 , y 3 ≤ y 2 , and y 4 ≤ y 1 so that D * = {(y 1 , y 3 ), (y 2 , y 3 ), (y 1 , y 4 )}. The maximal dom-CRN is This dom-CRN is not admissible since the domination reaction Since this dom-CRN does not contain any domination reactions which lead directly to the terminal complex X 1 , it is admissible. There are several options for external forests on this dom-CRN, including those indicated as follows in bold red: Note that, within the external forests above (bold red), there is a unique path from every complex to the terminal complex X 1 (shaded blue). In order to be balanced, we need to find satisying the conditions of Definition 2.6. For the external forest on the left, we have the following system of equalities and inequalities: which can be satisfied by the vector (1, 0, 1, 0, 1). For the external forest on the right, we which has no nontrivial solution. It follows that the external forest on the left is balanced while the one on the right is unbalanced. The system (4) corresponds to the one in the introduction.

Conditions for Extinction Events
In [13], the authors present and prove the following main results regarding extinction events in CRNs with discrete state spaces. When algorithmically establishing that a CRN on a discrete state space has an extinction event, we will use Corollary 2.1 (the contrapositive of Theorem 2.1). Note that it is sufficient to find a single Y-exterior forest which is unbalanced, even if the dom-CRN has other balanced Y-exterior forests.

Computational Implementation
In this section, we outline a computational algorithm capable of affirming whether a given CRN exhibits an extinction event according to Corollary 2.1. The program is written in Python and utilizes mixed-integer linear program (MILP) modules. MILP algorithms have been used increasingly within CRNT in recent years to verify a variety of structural properties of CRNs [10, 12, 14-16, 25, 27].
For brevity, we have placed technical discussion of the algorithm in the Appendix. In Appendix A, we present the background theoretical results required for the program to check the technical conditions of Corollary 2.1. In Appendix B, we present details on how the algorithm cycles through and expands the absorbing complex set Y. In Appendix C, we present a detailed description of the modules outlined in Section 3.1. In Appendix D, we summarize the output of our BioModels Database run.

Description of Program
The pseudocode given in Algorithm 1 implements Corollary 2.1. The modules are briefly described below. More detailed explanations are given in Appendix C.

CreateModel(N ):
This module takes in a CRN N = (S, C, R) and generates the structural matrices Y , I a , Γ, I s , and A.

IsSubconservative(N ):
This module determines whether the CRN is subconservative. If the CRN is subconservative, it further runs the module IsConservative(N ) to determine whether it is also conservative.  if IsSubconservative(N ) then

BioModels Database
We ran Algorithm 1 on 458 curated models from the European Bioinformatics Institute's BioModels Database. In total, 86 models were found which had a guaranteed extinction event according to Corollary 2.1. Of these, 53 models were conservative and 33 were subconservative but not conservative. We limited our run to those models which contained 50 or fewer reactions. We further classify the models which exhibit an extinction event according to the following: 1. A species X i is classified as source only if, in the unbalanced Y-exterior forest which guarantees extinction, X i appears only in source complexes.

2.
A species X i is classified as product only if, in the unbalanced Y-exterior forest which guarantees extinction, X i appears only in product complexes.
These classifications identify mechanisms whose primary purpose is to convert one substrate (source) into another (product). In such mechanisms, we expect discrete extinction events to occur as a consequence of a source being used up and/or a product being formed. For example, consider the classical Michaelis-Menten mechanism [21]: where the unbalanced exterior forest is indicated in bold red, and terminal complex is shaded blue. The CRN is conservative, and has the source only species S and product only species P . Notice that S is source only even though it appears as a product in the original CRN since the reaction SE → S + E is not contained in the exterior forest (red). The remaining species E and SE appear as both a source and product in the exterior forest. For large CRNs, the property of having source only or product only species may be difficult to verify directly.
Of the 86 models classified as having a discrete extinction event from the BioModels run, 81 have source only species, 57 have product only species, and 52 have both source only and product only species. These results are consistent with the observation that many models in the database model processive biochemical mechanisms like the Michaelis-Menten mechanism. A full enumeration and classification of the BioModels networks identified by the algorithm as having a guaranteed extinction event is given in Appendix D.
The results of the BioModels run has suggested some interesting CRNs with extinction events and some unique pathways for obtaining extinction. It also suggests avenues for future work. We now further investigate two models which the program identified as having an extinction event according to Corollary 2.1: a model of polyamine metabolism in mammals (biomd0000000190) [24], and a model of the pentose phosphate pathway in Trypanosoma brucei (biomd0000000513) [1,5,17].

Polyamine Metabolism
Consider the model of polyamine metabolism given in Table 1 which corresponds to model biomd0000000190 in the BioModels database [24].  The program further identifies M et and Orn as source only species which suggests that they are inputs to the system. Indeed, the results are consistent with the observation that the mechanism requires input of M et and Orn in order to maintain function as without these species all reactions except (12) will eventually shut down.

Pentose Phosphate Pathway
Consider the model for the pentose phosphate pathway in Trypanosoma brucei given in Table 2 which corresponds to model biomd0000000513 in the BioModels Database [1,5,17]. The algorithm presented in Section 3.1 identifies this mechanism as having a guaranteed extinction event.
(1) Glce ⇄ Glcc The algorithm identifies several structurally distinct pathways by which an extinction event can occur. The first is the following subnetwork of the dom-CRN, which incorporates reactions (9), (10), (19), and (20) as indexed in Table 2: The program identifies the complexes PEP c +ADP c , Pyr c +ATP c , Pyr c , 2ADP c , AMP c +ATP c , and ATP c as transient. It follows immediately that, at any extinction state, we have the following counts: Pyr c = 0, ADP c ≤ 1, and ATP c = 0. It can be furthermore seen from reactions (10) and (20) that ADP c = 0 is not absorbing, so that ADP c ≥ 1 infinitely often. It follows that ADP = 1 at the extinction states, from which it follows from the transience of PEP c + ADP c that PEP c = 0 at the extinction states.
We now construct the pathway to extinction. It follows from the conservation on the original CRN and the observation that Pyr e only appears as a product in any reaction that any reaction which forms Pyr e must have a final occurrence. Furthermore, we can exhaust PEP c by converting it into Pyr c through the forward reaction in (19) since ADP c ≥ 1 infinitely often, as previously argued. A less trivial pathway to extinction occurs for the species ADP c and ATP c . The conservation AMP c + ADP c + ATP c = total A suggests that the system may become locked by converting ADP c into AMP c and ATP c through reaction (20), and then converting ATP c into ADP c through reaction (10), and then repeating as many times as possible. Eventually we will arrive at a state where ADP c = 0, AMP c = total A − 1, and ATP c = 1, after which reaction (10) locks both reactions. We must, however, consider the possible that ATP c is converted to ADP c by another pathway and, in fact, the backward reaction in (19) is exactly such a reaction. It is also, however, the only such pathway. Consequently reactions (9) and (19) must shut down before reactions (10) and (20) can be shut down. The complete sequence of reactions required to shut down the indicated complexes is therefore: 1. Convert all possible substrates in PEP c or Pyr c and define total P =PEP c +Pyr c .
2. Fire the forward reaction in (19) to convert all PEP c into Pyr c , replenishing ADP c through reaction (10) as required.
3. Fire reaction (9) to convert all Pyr c into Pyr e . (10) to convert all ATP c into ADP c . (20) to convert all ADP c into AMP c and ATP c .

Fire the forward reaction in
6. Repeat steps 4. and 5. until you arrive at the state ADP c = 0, AMP c = total A − 1, and ATP c = 1, and then fire (10).
The final extinction state is PEP c = 0, Pyr c = 0, Pyr e = total P , AMP c = total A , ADP c = 1, and ATP c = 0. Notice that the extinction of PEP c is only guaranteed for all trajectories by the observation that ADP c ≥ 1 infinitely often. We also require that no more PEP c and Pyr c can be produced by the remaining pathways. This example illustrates that the program is able to identify transient complexes which might be very difficult to determine by direct analysis of potential pathways to extinction. The other subnetwork which was identified as leading to extinction is the following, which incorporates reactions (4) and (14) in Table 2: DHAP c +Gly-3-P g DHAP g +Gly-3-P c Gly-3-P c DHAP c 14b 14b 4 D The program identifies the complexes DHAP c +Gly-3-P g , DHAP g + Gly-3-P c , and Gly-3-P c as transient. It follows that at the extinction state we have Gly-3-P c = 0. It furthermore follows from the observation that DHAP c > 0 at the extinction state that we must have Gly-3-P g = 0. This is consistent with the observation that, after firing the forward reaction in (14) and then reaction (4) repeatedly to exhaust Gly-3-P g and Gly-3-P c , there is no mechanism by which to convert DHAP c or any of its derivatives back into Gly-3-P or any of its derivatives. Table 2 and this complex was not identified as transient by the program. From reactions (6), (7), and (8), we can furthermore identify the complexes PEP c , 2-PGA c , 3-PGA c , and 3-PGA g as transient since they can all be transformed in Pyr e . We have that the set of complexes which is necessarily transient at the extinction states is greater than that strictly guaranteed by Condition 2.1.

Conclusions
In this paper, we have presented a computational implementation of the conditions derived in the companion paper [13] for an extinction event in a CRN with a discrete state space. We have run the program on one of the most widely studied network databases in system biology, the European Bioinformatics Institute's BioModels database. This work has yielded some mathematically and biologically interesting pathways by which extinction in biological systems may be attained, as indicated by our analysis of the model of polyamine metablism in mammals and the model of the pentose phosphate pathway in Trypanosoma brucei.
The most notable avenue for future work opened up by the study is in extending the current results to determine the maximal set of transient complexes. Currently, even when Algorithm 1 guarantees an extinction event by Corollary 2.1, it does not necessary give the maximal such set. For example, for the polyamine metabolism model studied in Section 3.3, the species S is transient but was not identified as such. For the pentose phosphate pathway model studied in Section 3.4, PEP c is transient but not identified as such. Further work will investigate methods for expanding the transient complex set given by Corollary 2.1 and Algorithm 1 into the maximal such set.

A Background Results
We require the following results in order to computationally implement the modules of Algorithm 1. The proofs of Theorem A.1 and Theorem A.2 can be found later in this Appendix. The proof of Theorem A.3 can be found in the Appendix of [9]. Consider a CRN (S, C, R) and a Y-admissible dom-CRN (S, C, R ∪ D) where Y ⊆ C is an absorbing complex set on the dom-CRN. Let I a denote the adjacency matrix of the dom-CRN. Consider a subnetwork of the dom-CRN (S, C, R F ∪ D F ) where R F ⊆ R and D F ⊆ D, and suppose that, for all y i ∈ Y, y i is the source for exactly one reaction in the subnetwork (S, C, R F ∪ D F ). Then the following are equivalent: 2. There are no cycles in the Y-exterior portion of the reaction graph of (S, C, R F ∪ D F ). Proof of Theorem A.1. We prove Condition 1 ⇔ Condition 2 ⇔ Condition 3.

There is a vector
Condition 1 ⇒ Condition 2: Suppose that (S, C, R F ∪ D F ) is a Y-exterior forest and that there is a cycle in Y-exterior portion of the reaction graph. Since every Y-exterior complex has a path to Y, there is a complex in the cycle which has a reaction leading out of the cycle. From this complex, however, there are two distinct paths which lead to Y: (i) the direct path with the reaction leading out of the cycle as its first reaction; and (ii) the path which transverses the cycle first, then the path in (i). It follows that (S, C, R F ∪ D F ) is not a Y-exterior forest, which is a contradiction. It follows that Condition 1 implies Condition 2.
Condition 2 ⇒ Condition 1: Suppose there are no cycles in the Y-exterior portion of the reaction graph. Consider a Y-exterior complex y i ∈ Y. By definition, there is a unique reaction (y i , y j ) ∈ R F ∪ D F with y i as its source complex. If y j ∈ Y, we are done; if y j ∈ Y, however, we continue the path with the unique reaction (y j , y k ) ∈ R F ∪ D F , and so on.
Since there are a finite number of Y-exterior complexes, this process may not terminate before either reaching a Y-interior complex, or a complex already in the path, which creates a cycle and is therefore a contradiction. It follows that the path reaches Y and is unique. Since y i ∈ Y was chosen arbitrarily, it follows that (S, C, R F ∪ D F ) is a Y-exterior forest. We have shown Condition 2 implies Condition 1.
Condition 2 ⇒ Condition 3: Suppose there are no cycles in the Y-exterior portion of the reaction graph. We will inductively construct the required vector v = (v R , v D ) ∈ R r+d ≥0 . We start by setting v i = 0 for y ρ(i) → y ρ ′ (i) ∈ R F ∪ D F . Next, we define the set: We set V 0 = 1 and v i = V 0 for all i ∈ I 0 . That is, reactions which are not contained in the Y-exterior forest are assigned a weight of 0 in the vector v, and those which lead directly into Y from the Y-exterior portion are assigned a weight of 1.
We now inductively work out from Y in the Y-exterior forest. We define the following: Since the set of complexes is finite, there are a finite number of such sets which are nonempty: I 0 , I 1 , . . . , I nI −1 , I nI . Further, since every Y-exterior complex is the source for exactly one reaction in the Y-exterior forest and there are no cycles, each such complex is assigned to exactly one set I k , k ∈ {1, . . . , n I }. It follows that this is a partition of the Y-exterior complexes, i.e. I k1 ∩ I k2 = Ø for all k 1 , k 2 ∈ {1, . . . , n I }, k 1 = k 2 , and nI k=1 I k = supp(Y). Now define n k = |I k | and where δ > 0 is a constant to be a determined later. For each i ∈ I k , we set v i = V k . Now consider a complex y j ∈ Y and the reaction y ρ(i) → y ρ ′ (i) ∈ R F ∪ D F where ρ(i) = j and i ∈ I k . Notice that this outgoing reaction is unique while there are at most n k+1 reactions which lead to y j . We therefore have: We now choose δ > 0 so that V k n k − δ > 0 for all k ∈ {1, . . . , n I }. We have constructed a vector v ∈ R r+d ≥0 with the desired properties, which shows Condition 2 implies Condition 3.
Condition 3 ⇒ Condition 2: Suppose also that there is a cycle on the Y-exterior portion of the reaction graph. We denote this cycle by: We define the following sets: . . , n} | µ(j) = i for some j ∈ {1, . . . , l}} and I in (j) = {k ∈ {1, . . . , n} | y i → y j ∈ R F ∪ D F where i ∈ I cycle and j ∈ I cycle }.
Notice that, since every Y-exterior complex is the source for exactly one reaction, there are no reactions which lead out of the cycle.
. It follows that, for j ∈ {1, . . . , l}, we have: Taking the summation of (5) over j ∈ {1, . . . , l} gives where the inequality follows from v ∈ R r ≥0 . This is inconsistent with the requirement that [I a v] j < 0 for j ∈ {1, . . . , n} such that y j ∈ Y, so we have shown Condition 3 implies Condition 2.

Proof of Theorem A.2. Consider a vector
where ρ(k) = j. It follows that α satisfies Condition 3 of Definition 2.6.

B Expanded Absorbing Set Y
When applying Corollary 2.1, the choice of absorbing complex set Y ⊆ C is often nontrivial. In fact, it is sometimes possible to show an extinction event occurs for one choice of Y but not for another (see Example 3.7 of the companion paper [13]).
In this Appendix, we will consider an algorithmic way for starting with a minimal absorbing complex set Y ⊆ C and systematically making it larger. We motivate the procedure with the following example.
Example B.1. Consider the following CRN: It can be seen directly that the CRN exhibits an extinction event since all of the X 1 can be converted into X 2 through reaction 1. This shuts down reactions 1 and 3 and therefore eliminates the ability to replenish X 1 . We now attempt to obtain this result by applying Corollary 2.1. We need to first determine a Y-admissible dom-CRN for some absorbing complex set Y ⊆ C on the dom-CRN. We will attempt to determine the minimal such set Y. This corresponds to the terminal complex set for the dom-CRN with the maximal set D ⊆ D * . We have the following domination set: This produces the following dom-CRN: This dom-CRN is not admissible for Y = {X 3 } since the domination reaction X 3 +X 4 → X 3 leads to the terminal complex X 3 . It can, however, be made admissible by removing this reaction. We therefore start the process by considering the dom-CRN with The corresponding admissible dom-CRN is as follows, where Y is shaded blue, and a path common to all exterior forests is shown in bold red: There are several exterior forests which need to be checked. Since the path in bold red above is common to all exterior forests, however, they are all balanced by the following vector α = (α R , α D ) ∈ R 8 ≥0 : Consequently, we cannot conclude anything by Corollary 2.1. It is suggestive, however, that the vector α which balances the CRN does not have support on the portion of the CRN which experiences an extinction event (specifically, the reactions involving X 1 ). This suggests that we redefine the set Y to correspond to the portion of the dom-CRN where the balancing vector α has support. We therefore expand Y to include the complexes in the red path above: Notice that, given this absorbing complex set Y ′ , we must modify the set D in order to allow the resulting dom-CRN to be admissible. In particular, we may not include We have the reduced set The corresponding dom-CRN is as follows, where the complexes in Y ′ are shaded blue, and a Y ′ -exterior forest is shown in red (Y ′ -interior reactions omitted): Notice that the set Y ′ is absorbing complex set but is not terminal, and that there is a unique path in the Y ′ -exterior forest from each Y ′ -exterior complex to the set Y ′ . This Y ′ -exterior forest is balanced if we can find a vector α = ( Condition 1 and the first equality in Condition 2 implies that (α R ) 1 = (α R ) 3 = 0 so that the forest in unbalanced. It follows from Corollary 2.1 that the discrete state space system has a guaranteed extinction event and that the system is transient on This example shows that it is sometimes possible to apply Corollary 2.1 to guarantee an extinction event for one absorbing complex set Y ⊆ C but not for another. While this is a limitation for our ability to apply Corollary 2.1, this example also suggests the following method for determining an initial set Y and then systematically expanding it: 1. Initialize Y * = Ø and let D * denote the domination set of the original CRN. (c) Remove from D any domination reaction with its product complex in Y.
(d) Repeat (b) and (c) until the resulting dom-CRN is admissible.
3. For the Y-exterior forests on this dom-CRN, attempt to compute balancing vectors α = (α R , α D ) ∈ Z r+d ≥0 with the maximal support (Definition 2.6). 4. For a representative vector α from this set, let Y * correspond to all the complexes corresponding to either the source or product complex of a reaction with α k > 0.
5. Repeat steps 2-4 until either an unbalanced Y-exterior forest is found, or Y = C.
Note that Condition 2(a) is only required if Y initially is not an absorbing complex set. It is natural to wonder whether it is necessary to the continually refine the sets Y and D to ensure that the resulting dom-CRN is admissible. The following example shows that an incorrect conclusion may be reached if admissibility is not taken into account.
Example B.2. Consider the following CRN, which corresponds to Example 3.10 of the companion paper [13]: We attempt to apply the steps contained above to seed the sets Y and D and then systematically refine them. We start with the set We set Y to correspond to the terminal complexes in the resulting dom-CRN, which gives Y = {X 2 }. Since the domination reaction X 2 + X 4 → X 2 leads to Y, we must remove it from D in order for the resulting dom-CRN to be admissible. The resulting dom-CRN is given by the following, where the only Y-exterior forest is indicated in red, and the set Y is indicated in blue: The Y-exterior forest indicated in red can be balanced by α = ((α R ) 1 , (α R ) 2 , (α R ) 3 , (α R ) 4 , (α D ) 1 , (α D ) 2 ) = (1, 1, 0, 0, 1, 0). While we cannot apply Corollary 2.1 directly, we notice that α does not have full support on the reactions in the Y-exterior forest above. We may therefore continue the algorithm presented earlier in this Appendix. We start by including in Y all those complexes which are contained in the reactions with α k > 0. This gives Y = {X 2 + X 4 , X 1 + X 4 , X 1 , X 2 }. Now suppose we omit steps 2(a − d) in the algorithm presented earlier in this section. The resulting dom-CRN would be the following, where the only Y-exterior forest is indicated in red (Y-interior reactions not shown), and the absorbing complex set Y is shaded blue: This Y-exterior forest cannot be balanced, which seems to suggest that the discrete state space system exhibits an extinction event by Corollary 2.1. We can, however, see that this dom-CRN is not admissible. To apply steps 2(a-d) of the algorithm presented earlier, we omit all domination reactions so that D = Ø (step (d)). We must then add X 1 + X 5 and X 3 + X 5 to Y since these complexes are now terminal in the resulting dom-CRN. Since we then have Y = C (step (e)), we terminate the procedure. In fact, this example does not exhibit an extinction event for most initial conditions. Provided X 4 > 0, X 5 > 0, and any one of X 1 , X 2 , and X 3 is positive, every complex is recurrent. This example therefore highlights the importance in guaranteeing that the resulting dom-CRN is admissible when expanding the absorbing complex set Y. Care must be taken when expanding the absorbing complex set Y. An example which is structurally identical, but which has a guaranteed extinction event, is given by Example 3.9 in the companion paper [13].

C Description of Algorithm 1 modules
In this Appendix, we present the details of the modules which are used in Algorithm 1. Several of the modules involve mixed-integer linear programming (MILP). Recall that a MILP problem in the vector of decision variables x = (x 1 , . . . , x k ) ∈ R k may be written in the form where c ∈ R k , A ∈ R p×k , B ∈ R q×k , a ∈ R p , and b ∈ R q [27].

CreateModel(N ):
This module takes the CRN N and generates the structural matrices Y , I a , Γ, I s , and A. We utilize the CRN module package created for Python by Elisa Tonello which reads in a string file of reactions, and can also interface with SBML files [28].

IsSubconservative(N ):
This module determines whether a given CRN N is subconservative or not by running the following MILP: over the vector of decision variables c = (c 1 , . . . , c m ) ∈ R m ≥0 . The first constraint set guarantees all components of c are strictly positive if they can be, while the second constraint set imposes that the CRN is subconservative. If (7) returns an optimal value less than zero, the module returns true; otherwise, it returns false. If the CRN is subconservative, it furthermore runs the module IsConservative(N ), which replaces the second constraint set in (7)

CycleForests(dom-N ):
This module cycles over every combination of reactions F such that every Y-exterior reaction in the dom-CRN is the source for exactly one reaction. The exterior forest conditions of Theorem A.1 are checked by the module IsExtForest(F ). If this module returns true, the module IsBalanced(F ) is run to determine if the forest is balanced in accordance with Definition 2.6 by Theorem A.2. If the module finds an exterior forest which is not balanced then the module confirms that the discrete state space CRN exhibits a extinction event in accordance with Corollary 2.1. If F is balanced, it runs ExpandedY(F); otherwise, the module returns false.

IsExtForest(F ):
This module checks whether F is an exterior forest (Definition 2.5) according to Theorem A.1. The module takes in the adjacency matrix I a for the dom-CRN and F . We let F i = 1 if R i ∈ R F ∪ D F and F i = 0 otherwise, and set a small parameter ǫ > 0. We then run the following linear program: over v = (v 1 , . . . , v n ) ∈ R n ≥0 . The first two constraint sets in (9) ensure that v has the same support on R F ∪ D F . The third constraint set guarantees that, if y i ∈ Y (i.e. Y i = 0) then [I d v] i ≤ −ǫ < 0 (there is no restriction if y i ∈ Y, i.e. Y i = 1). If (9) is feasible, the module returns the value true; otherwise, it returns false.

IsBalanced(F ):
This module checks whether F is balanced (Definition 2.6) according to Theorem A.2. It takes in the stoichiometric matrix Γ of the CRN, the adjacency matrix I a of the dom-CRN, and the exterior forest F . We let F i = 1 if R i ∈ R F ∪ D F and F i = 0 otherwise, and set a small parameter ǫ > 0. We then run the following linear program: over α = (α 1 , . . . , α r+d ) ∈ R r+d ≥0 . The first constraint guarantees that R i ∈ R F ∪ D F (i.e. F i = 0) then α i = 0. The second constraint set guarantees that Γ α R = 0. The third constraint set guarantees that, if y i ∈ Y (i.e. Y i = 0) then [I d a α] i ≤ 0 (there is no restriction if y i ∈ Y, i.e. Y i = 1). If (10) the optimal value zero, corresponding to the trivial vector α R = 0, the module returns false; otherwise, it returns true. It also returns the vector α.

ExpandedY(F ):
For a given vector α = (α R , α D ) ∈ Z r+d ≥0 , this module returns the set of complexes in the reactions R k = (y i , y j ) ∈ R F ∪ D F such that α k > 0. Note that if α has full support, then we set Y = C, which is one of the stopping criteria for Algorithm 1.

WriteOutput:
This module writes the output into a .dat file.

D Result of BioModels Database Search
The following 86 models from the European Bioinformatics BioModels Database were identified by the algorithm as being subconservative and having a discrete extinction event. The models are labeled with (C) if they were conservative as well as subconservative, (S) is they had a source only species, and (P) if they had a product only species.