Reachability Analysis of Low-Order Discrete State Reaction Networks Obeying Conservation Laws

In this paper we study the reachability problem of suband superconservative discrete state chemical reaction networks (d-CRNs). It is known that a subconservative network has bounded reachable state space, while that of a superconservative one is unbounded. The reachability problem of superconservative reaction networks is traced back to the reachability of subconservative ones. We consider network structures composed of reactions having at most one input and one output species beyond the possible catalyzers. We give a proof that, assuming all the reactions are charged in the initial and target states, the reachability problems of suband superconservative reaction networks are equivalent to the existence of nonnegative integer solution of the corresponding d-CRN state equations. Using this result, the reachability problem is reformulated as an Integer Linear Programming (ILP) feasibility problem. Therefore, the number of feasible trajectories satisfying the reachability relation can be counted in polynomial time in the number of species and in the distance of initial and target states, assuming fixed number of reactions in the system.


Introduction
Employing deterministic ordinary differential equation systems to characterize the dynamical behavior of complex networks of chemically interacting components (species) is a widely used approach in mathematical and computational systems biology [1][2][3].Such a continuous state modeling approach assumes high molecular count of species and their homogeneous (well-mixed) distribution in the surrounding media [4].However, in several (bio)chemically interesting systems, such as some enzymatic and genetic networks, the molecular count of different species is relatively low (e.g., < 100 molecules) [4][5][6] implying that the assumption of homogeneous species distribution does not hold [7,8].Hence it is necessary to introduce a discrete state model capable of keeping track of the individual molecular counts in order to properly characterize the qualitative dynamical behavior of (bio)chemical networks of species with low number of molecules [9,10].There exist several mathematical models describing the state evolution of discrete state chemical reactions networks, such as Markov chain models [8,10] and stochastic Petri nets [11].
In the context of chemical reaction networks of several interacting components, in order to completely characterize the system it is needed to simultaneously study the dynamical behavior and the underlying network structure as well.Moreover, it is also important to examine how the dynamical behavior and the network structure are related to each other, and how we can predict the dynamical behavior (e.g., in the form of possible state space trajectories) and be aware of the underlying network structure.For continuous state reaction networks obeying the law of mass action, it is recognized that the network structure (i.e., topology) is not necessarily unique; i.e., the same system of differential equations can be generated by different network topologies (different sets of interactions among the given species) [12][13][14][15].the set of non-negative integer numbers T × the set of ( × )-dimensional vectors over the set T 0 × a zero matrix of dimension  ×  1 ×  a matrix of dimension  ×  for which all the entries are equal to 1 {0, 1} ×  the set of ( × )-dimensional binary vectors (all the entries are equal to 0 or 1) {−1, 0, 1} ×  the set of ( × )-dimensional vectors composed of the entries −1, 0, 1 [] ,: the th row of the matrix   ≺  for ,  ∈ R  ,   <   for  = 1, . . .,   ⪯  for ,  ∈ R  ,   ≤   for  = 1, . . .,    an ordered sequence of states   an ordered sequence of reaction vectors   an ordered sequence of species   an ordered sequence of complexes In the case of discrete state reaction networks the socalled reachability is a strictly related problem to the dynamical behavior; namely, is it possible to reach a prescribed target state from a given initial one through a finite sequence of transition (reactions)?It is known that the reachability relation between any pair of nonnegative initial and target states is determined by the network structure itself.Through the reachability analysis several problems of great importance can be analyzed; one of them having high interest is the existence of so-called extinction events: the existence of trajectories resulting in the irreversible extinction of some species from the system.It has been shown that under some conditions on the network structure a discrete state chemical reaction network exhibits an extinction event from any point of its state space [9,16,17].The properties of recurrence (the ability of returning to any initial state) and irreducibility (the ability of reaching any state from any other one) are also examined in the context of discrete state reaction networks [18,19].
The mathematical model of discrete state chemical reaction networks is equivalent to an important model of theoretical computer science, namely, the so-called vector addition systems with states (VASS) or equivalently Petri nets [20,21].Hence the discrete chemical reaction network reachability problem is equivalent to the extensively studied problem of vector addition system (VAS) reachability.The VAS reachability problem is known to be decidable [22][23][24][25], and for the space complexity we have EXSPACE lower bound [26].Unfortunately, contrary to the proven polynomial time complexity of reachability of rate independent continuous state chemical reaction networks [21], in the case of discrete state reaction networks it is not known whether there exists an algorithm of primitive-recursive time complexity deciding this problem [27].
The aim of this paper is to study of the reachability problem of sub-and superconservative d-CRNs.We make use of the relation between the sub-and superconservative properties.In Propositions 15 and 17, we give necessary and sufficient conditions on the network structure and the initial and target states under which the reachability is equivalent to the nonnegative integer solution of the d-CRN state equation.Then these results in Corollaries 16 and 18 are extended to a subclass of superconservative d-CRNs.
The paper is organized as follows.In Section 2 the necessary mathematical notations and concepts of Chemical Reaction Network Theory (CRNT) are introduced.Section 3 discusses the classes of sub-and superconservative d-CRNs and their duality as well.In Section 4 the reachability problem of sub-and superconservative d-CRNs is examined.Firstly the case of low state space-dimensional d-CRNs is discussed, followed by the extension to the general case when the dimension of the state space is arbitrarily high.In Section 5 our findings are illustrated in a representative example.

Notations and Mathematical Background
In Table 1 we summarize the notations and concepts of discrete chemical reaction networks which will be extensively used later.

Discrete State Chemical Reaction Networks.
A discrete state Chemical Reaction Network (d-CRN) with  species,  complexes, and  reactions is a triple N = (S, C, R) so that where   is the 'th species,   is the 'th complex, and  V is the V'th reaction of the network.Moreover,   is the stoichiometric coefficient of the 'th species in the 'th complex.For a reaction and  ( V ) are the source complex and the product complex, respectively.For each complex   ∈ C,  ∈ {1, . . ., }, the stoichiometric coefficients of the species can be represented as a vector of the following form: For each  ∈ R, a reaction vector   ∈ Z  can be associated with the track of the net molecular count changes of the species upon firing the reaction: so that   and   are the corresponding source and product complexes of .In the sequel the notation   will be used for denoting both the 'th reaction of the d-CRN and the associated reaction vector as well.We will also assume that for all the examined d-CRNs a fixed order of the reaction vectors is given; i.e., an order  1 ,  2 , . . .,   is fixed and  = |R|.A d-CRN can also be represented by a directed graph  = (, ) such that the vertices and edges correspond to the complexes and the reactions, respectively, i.e.,

𝑉 = C
(4) The direction of the edges is determined by the reactions of R, so that if   →   ∈ R, then there exists an edge  ∈  from the vertex representing   to the vertex of   .
For each edge a weight corresponding to the reaction rate constant (also called intensity or propensity) corresponding to the respective reaction can also be associated.
Beyond the above representations it is also possible to describe a d-CRN in an algebraic way by means of its unique stoichiometric matrix.Definition 1.Let us consider a d-CRN N = (S, C, R).The stoichiometric matrix Γ ∈ Z × of N is defined as follows: Note that [Γ]  encodes the net molecule count change on species   upon the occurrence of reaction   .Besides Γ we also define the following matrices: where  +   denotes the vector form of the product complex belonging to reaction   while  −   represents the vector of the source complex associated with reaction   .The relation among the above defined matrices is as follows: We fix the order of species and reactions as they are listed in the above sets.
N has no information on the probabilities of the reactions, but at any given time instant at most one reaction can occur.
The molecular count of each species of a d-CRN at any time  ≥ 0 is given by its state vector () ∈ Z  ≥0 and the time evolution of the system is characterized by the following discrete state equation: where (0) is the state vector belonging to the initial time instant and () = [ 1 (),  2 (), . . .,   ()] ⊤ ∈ Z  ≥0 such that   () ∈ Z ≥0 stores the number of occurrences of the 'th reaction up to time .We note that () is typically modeled as some point process [8,10].For our further analysis the time instants when the reactions have occurred are not of interest, but only the order of reactions; therefore we abandon the notation of time  in the formulas.
The condition that a reaction  ∈ R is charged at state  ∈ Z  ≥0 can be expressed by the inequality  ⪰  −  .For a reaction sequence   a state transition sequence   =  0  1 . . . V can be uniquely associated so that where the initial state  0 is assumed to be given.A state transition sequence   is said to be admissible if   ⪰  +1 for   ∈   ,  ∈ {0, . . ., V − 1}; moreover, we say that a reaction sequence   is admissible if the corresponding state transition sequence is admissible.
From the reachability of a state   ∈ Z  ≥0 from an initial state  0 ∈ Z  ≥0 , it follows that the following equation has a nonnegative integer solution  ∈ Z  ≥0 : where []  encodes the number of occurrences for reaction   ∈ R for  ∈ {1, . . ., }.However, it is important to note that from the existence of a nonnegative integer solution  of ( 13), the reachability relation  0  N   does not necessary follow.
We note that  of ( 13) corresponds to () of ( 11).Since a solution  ∈ Z  ≥0 of ( 13) encodes the number of occurrences for each reaction in a fixed order, the following equality is fulfilled: where ℎ = ∑  =1 []  and   ∈ R for  ∈ {1, . . ., ℎ}.When we want to emphasize that a reaction vector sequence is encoded by a particular  ∈ Z  ≥0 , we will use the notation    =  1 , . . .,  ℎ and the state transition sequence determined by    will be denoted by    .
Definition 4. Let us consider a d-CRN N with stoichiometric matrix Γ ∈ Z × and an initial state  0 ∈ Z  ≥0 .The reachable state space ℎ(N,  0 ) of N with initial state  0 is the set of nonnegative discrete states reachable from  0 .

Integer Linear Programming
In this section some relevant concepts of mathematical programming that will be extensively employed later are briefly reviewed.An Integer Linear Programming (ILP) instance can be formulated as follows: where  is the -dimensional vector of decision variables while  ∈ Z  ,  ∈ Z × , and  ∈ Z  are fixed coefficients.Generally, the above ILP computational problem is known to be NP-hard, which may highly confine our ability to efficiently solve problems of integers in high dimension.However, if the value of the decision vector that minimizes (or maximizes) the prescribed objective function is not important for us, but only the existence of a  ∈ Z  vector satisfying the set of specified constraints, then the problem is called ILP feasibility problem.
An ILP feasibility problem, as a decision problem, addresses the question of whether the polytope  contains an integer lattice point, formally  ∩ Z  ?= 0.While a FP is also known to be NP-hard, it has well-decoupled time complexity with respect to the number of variables, the number of constraints, and the maximum of the absolute values of the entries of  and .Therefore, a feasibility problem of the form (17), assuming fixed dimension , can be decided in polynomial time in the number of constraints  and the maximum of the absolute values of the coefficients  and  by means of the Lenstra algorithm [28,29].Moreover, the number of integer lattice points in  can also be numerated in polynomial time in  and the maximum of the absolute value of the coefficients using Barvinok's integer lattice point counting algorithm [30][31][32][33].We note that for the Barvinok algorithm there exists an effective implementation called LattE [34].

Sub-and Superconservative d-CRNs
We define conservativity and subconservativity in the same way as they were introduced, e.g. in [4,16].
An important property related to subconservativity is the strong boundedness which is defined as follows.Definition 6.A d-CRN N is said to be strongly bounded if, for any  0 ∈ Z  ≥0 initial state, the reachable state space ℎ(N,  0 ) is bounded.
The subconservative property of the reaction network structure is a necessary and sufficient condition of strong boundedness [16,35].
As a special case covered by the intersection of sub-and superconservativity, we can define the conservative property as well.Definition 8. Let us consider a d-CRN N = (S, C, R) with stoichiometric matrix Γ ∈ Z × .The d-CRN N is said to be conservative if there exists a vector  ∈ R  >0 satisfying the matrix equation  ⊤ Γ = 0 1× .We note that the above structural properties can be easily decided in polynomial time by means of an LP of the following form: The relationship between sub-and superconservativity can be expressed by the following proposition.
We note that −Γ N means the change of the direction of each reaction in the d-CRN N of stoichiometric matrix Γ N .
Example 10. Figure 2 depicts two d-CRNs: a subconservative and a superconservative reaction network structure.Indeed, these networks are counterparts that can be easily transformed to each other by changing the sign of the entries in the stoichiometric matrices.Such a transformation results in the change of the direction of the edges in the reaction network.
From Proposition 9 it follows that, instead of the reachability problem of a superconservative network structure, one can consider an equivalent subconservative d-CRN reachability problem as is discussed in Proposition 11.The importance of Proposition 11 is that the reachability problem of a superconservative d-CRN of unbounded reachable state space can be easily traced back to the reachability problem of a d-CRN of bounded reachable state space which can make the original decision problem computationally tractable.

Low-Dimensional Case.
In this section the case of lowdimensional (rank(Γ) ≤ 2) subconservative d-CRNs is considered.We state a modified version of Proposition 5 of [36] where the conditions on the initial and target states are less strict.Then we extend the result to superconservative d-CRNs.
In order to discuss low-dimensional reachability problems, we introduce a distinguished state  = (Γ − ) as follows: Here Γ − is defined by (8).Note that the set { |  ∈ Z  ≥0 ,  ⪰ } contains all the states where each reaction is charged.
Let us assume that there exists a transition state   ,   ⪰ , so that the forthcoming state  +1 satisfies the inequality [ +1 ]  < []  for some  ∈ {1, 2}.For the target state   to be reached the inequality   ⪰  holds; hence there exists a reaction increasing the state variable along the coordinate .Let us assume that all the reactions increasing the state variable along   decrease the other coordinate   so that the resulting forthcoming state  +1 satisfies the inequality [ +1 ]   < []   .Then   =  holds.Now there are two different cases: ( 1 ) If   = , then Algorithm 1 terminates, and the correctness follows.( 2 ) If   ̸ = , then the subconservativity of N implies that it is not possible to reach a state ,  ⪰ ,  ̸ = ; i.e.,   is not reachable from   .This is contradiction, since arbitrary permutation of the initial ordering   results in the same target state   , given the initial state  0 .Then the correctness of Algorithm 1 follows.Algorithm 1 can be easily extended to the class of superconservative reaction networks.

Corollary 13. Let us consider a superconservative d-CRN N
with stoichiometric matrix Γ ∈ {−1, 0, 1} × and Γ − ∈ 0, 1 × .Assume that rank(Γ) ≤ 2 holds and consider an initial state  0 ∈ Z  ≥0 and a target state   ∈ Z  ≥0 for which  0 ⪰  and   ⪰  hold where  is defined by (20).Then the state   ∈ Z  ≥0 is reachable from  0 if and only if the equation has a nonnegative integer solution .
Proof.According to Proposition 11 we can consider a subconservative d-CRN N  of stoichiometric matrix −Γ and take the reachability problem   ? N   0 .Then Proposition 7 can be applied.

Sub-and Superconservative d-CRNs of Arbitrary High
State Space Dimension.In this section the reachability problem of arbitrary high-dimensional sub-and superconservative d-CRNs is considered.Firstly we examine network structures composed of reactions having at most one input and one output species.It is shown by an inductive proof that, under some auxiliary condition, the reachability relation  0  N   is equivalent to the existence of a  ∈ Z  ≥0 solution of the d-CRN state equation  0 + Γ =   .Then, according to the relation between sub-and superconservative reaction network structures, this result is generalized to a subclass of superconservative d-CRNs as well.We also extend the results to d-CRNs containing second-order reactions by allowing catalyzer species.
Firstly, we adopt the following necessary and sufficient condition of reachability from the theory of Petri nets (see Theorem 16, [37]) which will be extensively used in the sequel.≥0 so that  0 ⪰  and   ⪰  hold where  = (Γ − ) is defined by (20).Then the reachability relation  0  N   holds if and only if there exists a vector  ∈ Z  ≥0 satisfying the state equation  0 + Γ =   . Proof.

. 𝑠 ](ℎ) ). Let us consider an arbitrary
Now we construct a d-CRN N  = (S  , C  , R  ) from the stoichiometric matrix Γ ∈ {−1, 0, 1} (−1)× and Γ − ∈ {0, 1} (−1)× as follows: and Here   = min{ 1 ,  2 } and   = max{ 1 ,  2 }.This way we obtained a d-CRN N  satisfying the assumptions of the proposition.Figure 3 gives an illustrative example of how N  is constructed.Now we assign to each   ∈ R  the ordered pair of source complex and product complex of  ∈ R from which it is obtained.
In such a way every reaction of N  is uniquely described by an ordered pair (  , ) ∈ R  × R.
Then by the mapping ((  , )) =  one can uniquely determine the reaction  ∈ R from which   ∈ R  is derived.
consider an admissible reaction vector sequence    associated with relation (28).Since for each   ∈ R  we associated the reaction  ∈ R from which   is obtained, making use of the mapping  : R  ×R → R, we can consider the reaction vector sequence   ( ∈ R ∀ ∈   ) uniquely determined by    .We start from  0 and modify the state variable  ∈ Z  ≥0 according to the reaction vector sequence   .We may get to two invalid cases: ( 1 ) []  2 = 0, but the source complex of the forthcoming reaction   ∈   is   2 .Then, according to the ( − 1)-dimensional reachability, it is guaranteed that   1 is charged at the current state .Let us insert   1  2 into   before the current reaction   .( 2 ) []  1 = 0, but the source complex of the forthcoming reaction   ∈   is   1 .Then, according to the ( − 1)-dimensional reachability, it is guaranteed that   2 is charged at the current state .It is known that   1 can be reached from   2 along a reaction vector sequence  *  in the reaction network of N. Let us insert  *  into   before the current reaction   .By modifying   according to the above discussed cases ( 1 ) and ( 2 ), we obtain an admissible reaction vector sequence    with respect to the reachability relation where  * ⪰ 0  , [ * ]  = [  ]  for  ∈ {1, . . .},  ̸ =  1 and  ̸ =  2 ; moreover, [ * ]  1 + [ * ]  2 =   .According to the assumptions N contains directed paths both from   1 to   2 and from   2 to   1 ; hence the reachability relation  *  N   follows.Then, due to the transitivity of the relation  N , we have that  0  N   also holds.
Proposition 15 can be extended to the case of superconservative d-CRNs.

Corollary 16. Let us consider a superconservative d-CRN
Let us consider two states  0 ,   ∈ Z  ≥0 so that  0 ⪰  and   ⪰  hold where  = (Γ − ) is defined by (20).Then the reachability relation  0  N   holds if and only if there exists a vector  ∈ Z  ≥0 satisfying the state equation  0 + Γ =   .
Proof.By changing the sign of the entries in the stoichiometric matrix Γ, we get a subconservative d-CRN N  of stoichiometric matrix −Γ.Then we can consider the reachability problem   ? N   0 .
We can extend Proposition 15 by allowing the restricted application of catalyzer species as follows.
( According to the duality of the sub-and superconservativity properties, we can extend Proposition 17 to the case of superconservative d-CRNs.Let us consider two states  0 ,   ∈ Z  ≥0 for which  0 ⪰  and   ⪰  where  = (Γ − ) is defined by (20)  By the above corollary, any reachability problem on a superconservative d-CRN satisfying the conditions of Corollary 18 can be easily traced back to that of a subconservative network; hence the problem is equivalent to finding a  ∈ Z  ≥0 solution for the respective d-CRN state equation.
The reaction network class covered by the above statements might be beneficial in modeling first-and secondorder (bio)chemical reaction networks.For a representative example, see Example 19 below.We also note that any mass action type chemical reaction network can be dynamically described by an appropriately constructed reaction network containing at most second-order reactions [39].Moreover, the hypergraph representation of chemical reaction networks (see, e.g., [40]) is helpful for checking the conditions of Proposition 17.
Example 19.Nuclear factors of activated T-cells (NFAT) are proteins that can exist in highly phosphorylated states [38].They act as transcription factors; i.e., they have regulatory role in transcription.NFAT1, which is a member of the NFAT family, has 13 residues that can be dephosphorylated upon stimulation.NFAT1 has two different states: active and inactive.The transition between active and inactive states of the protein is regulated by the level of phosphorylation such that the higher the level of phosphorylation is, the lower the rate of transition becomes from inactive state to the active one and vice versa.Phosphorylation and dephosphorylation are achieved by a kinase and calcineurin, respectively.In the mathematical model the activities of kinase and calcineurin are modeled as rate constants; hence the respective reactions can be considered as first-order ones.The protein might be located in the cytoplasm or the nucleus of the cell.Cytoplasmic active NFAT1 is imported to the nucleus, while inactive NFAT1 of the nucleus is exported back to the cytoplasm.
The reaction network structure is depicted in Figure 4.It is visible that each reaction is of first order and there is no degradation and synthesis; hence the reaction network structure is conservative with a particular conservativity vector  = 1 56 and Proposition 15 can be applied.
We note that a reachability problem of the discussed reaction network class without additional constraints may be determined in polynomial time [41].However, by using an ILP feasibility approach, the number of all distinct trajectories satisfying a prescribed reachability relation can be determined efficiently (see Remark 20), assuming the fixed number of reactions in the network.In addition, the ILP formulation can also be equipped with further linear constraints. N   is equivalent to the existence of a nonnegative integer solution  ∈ Z  ≥0 of the state equation  0 + Γ =   .In this way the reachability problem can be reformulated as an ILP feasibility problem in terms of , and the Barvinok algorithm can be applied.Using the Barvinok algorithm in this particular case, the following complexity bounds are obtained: (1) exponential in the dimension of the decision variables, that is, in the number of different reactions , (2) polynomial in the number of constraints, that is, in the number of species , (3) polynomial in the maximum of the absolute values of the coefficients Γ,   −  0 .
The particular importance of Remark 20 is that the time complexity of the trajectory counting problem between a prescribed pair of states is polynomial in the number of constraints and in the distance of the initial and target states even in the case of superconservative d-CRNs for which the associated reachable state space can be unbounded for any  0 initial state.

Computational Example: A Superconservative d-CRN of First-Order Reactions
Let us consider the d-CRN depicted in Figure 5.This system is superconservative with a particular conservation vector  = 1 21×1 implying the unboundedness of its reachable state space regardless of the initial state  0 .Making use of the above results, the reachability problem of  0 ?   for any  0 ,   ∈ Z 21 ≥0 can be reformulated as a subconservative d-CRN reachability problem for which the boundedness of the reachable state space, i.e. structural boundedness, is guaranteed and is equivalent to the existence of a nonnegative integer solution of the respective subconservative d-CRN state equation.
As initial state we consider  0 given by (31) that was randomly generated from [10,100]  21 .In order to find a target state   satisfying the reachability relation  0    we randomly generated target states so that the number of each species was uniformly sampled from the interval [40,100].In the choice of the intervals from which we sample, it was taken into consideration that the discrete state model of reaction networks is typically employed in the case of low molecular counts [4,5].In order to decide the reachability relation between a pair of particular states  0 and   , we need to solve the following decision problem.

Complexity
We found the target state   given by (32) the reachability relation holds.To solve the decision problems of the form (30) the LattE [34] software was used.Now, let us examine the reachability from  0 to   with additional constraints.One can observe that   results in a significant increase in the number of molecules in the species  5 ,  6 ,  20 , and  21 and any trajectory from  0 to   results in a net increase in the number of molecules.These together imply the flow of molecules from the zero complex (environment).The flow of molecules over the network from the zero complex to  5 ,  6 ,  20 , and  21 can take place through different paths.We assume that the directed paths 2 =  3  16  17  4 (34) are slow compared to the other ones; hence we wish to minimize the flow through them in order to lower their effect in .This can be easily expressed by posing addition linear constraints on  as is done in the decision problem (35).
We also determined a particular solution  by equipping (35) with the objective function ∑ For implementation purposes we employed Python 2.7 programming language and the Gurobi mathematical optimization solver [42].A Lenovo P51s workstation with two 2.70 GHz i7-7500U CPUs and 32GB RAM (DDR4 2133 MHz) was used for all the computations.

Conclusion
In this paper the reachability problems of sub-and superconservative discrete state chemical reaction networks are considered.It is shown that the reachability problem of a superconservative reaction network of unbounded reachable state space can be transformed to that of a subconservative network for which the boundedness of the reachable state space is always guaranteed.Using an inductive proof we provided a set of necessary and sufficient conditions under which the equivalence between a d-CRN reachability problem and existence of nonnegative integer solution of the corresponding state equations is guaranteed.In such a way the reachability problem can be traced back to an IP-feasibility (decision) problem for which the number of decision variables is significantly lower than that employed in the literature [36].Moreover the number of trajectories satisfying the reachability relation can also be enumerated efficiently, assuming a fixed reaction network structure.The applicability of our approach is illustrated on a highdimensional superconservative d-CRN.

1 Figure 1 :
Figure 1: A discrete state chemical reaction network.Left: reaction network structure.The nodes and directed edges represent the complexes and the reactions, respectively.The numbers on the edges denote a fixed ordering of the reactions.Right: the stoichiometric matrix associated with the system, i.e. [Γ]  is the net change in the number of the 'th species upon occurring the 'th reaction.

Figure 2 :
Figure 2: A pair of sub-and superconservative reaction network structures denoted by N and N  , respectively.The ordering of the reactions are denoted by the numbers on the edges of the graphs.The two networks can be transformed to each other by changing the sign of the entries in their stoichiometric matrices.(a) Subconservative d-CRN.(b) Superconservative d-CRN.

(a) 𝑘 = 2 Figure 3 :
Figure 3: Graphical explanation of how the reaction network structure of N  in the proof of Proposition 15 is constructed.(a) Reaction network structure of an -dimensional d-CRN N. (b) Reaction network structure of N  resulting from merging the species   1 and   2 of N along their shared reaction   1  2 (and reverse counterpart reaction   2  1).Note that by merging   1 and   2 we obtain a stoichiometric matrix Γ  having redundant reactions (e.g., ( 1 ,   1 ), ( 1 ,   2 ) resulting in ( 1 ,    1 ), ( 1 ,    1 )) and zero reaction vectors (i.e., self-loops on    1 ), but they are omitted in (b).A directed cycle on which the chosen reaction   1  2 lies is depicted in gray.

Figure 4 :
Figure 4: Conformational switch model of NFAT1[38].Lower case letters denote the protein located in the cytoplasm while upper case letters refer to the protein in the nucleus.  ,   and   ,   for  = 0, . . .13 denote the active and inactive proteins, respectively.Lower indices denote the number of phosphorylated residues.

Remark 20 .
Let us consider a subconservative (superconservative) d-CRN N = (S, C, R) of  species,  complexes, and  reactions.Assume that N satisfies the conditions of Proposition 17 (Corollary 18).Then for any  0 ,   ∈ Z  ≥0 initial and target states for which  0 ⪰ (Γ − ),   ⪰ (Γ − ) hold we have that the number of distinct trajectories   satisfying the reachability relation  0 ? N   can be determined in polynomial time in the distance of  0 and   , given the fixed number of reactions  in the d-CRN.The explanation of this is the following.According to Proposition 17 (Corollary 18) the reachability problem  0 ?

Figure 5 :
Figure 5: A superconservative d-CRN.0 indicates the zero complex and the numbers denote the indices of the reactions on which they are located.Due to the superconservativity of the network structure, the above d-CRN is unbounded for any initial state  0 ∈ Z  ≥0 .
)  →   ,   →  +   ,   →  + } ) if there exists a reaction  ∈ R such that  is charged at state  and  +  =   , (5) a reaction (vector) sequence   is an ordered set of reaction vectors  =  1 . ..V where   ∈ R,  = 1, ..., V,(6)a state transition sequence   is an ordered set states  0 ,  1 , . . .,   so that  1 → 2 → . . .→  −1 →   , (7) a state   ∈ Z  ≥0 is reachable from a state  ∈ Z  ≥0 (denoted by  N   ) if there exists a directed path in the state space so that us take an initial state  0 ∈ Z  ≥0 and a target state   ∈ Z  ≥0 .Then the reachability  0  N   holds if and only if    N   0 also holds. 0  N   ⇒ ∃ ∈ Z  ≥0 such that  0 +Γ =   which is equivalent to   + (−Γ) =  0 .From  0  N   it follows that the solution  ∈ Z   ; i.e., all the states of    determined by    are composed of nonnegative entries.Then, by reversing    , we obtain a nonnegative state transition sequence σ  from   to  0 which is uniquely determined by means of the reaction vector sequence σ  = −  ℎ . . .−   1 .It is also needed to show that σ  is an admissible reaction sequence.This can be done as follows: for each state  ∈    \ 0 there exists a reaction  ∈    so that upon firing  the resulting state is , from which it follows that  ⪰  +  ; moreover, considering the reversed reaction sequence σ  , the reaction vector that will occur at state  is − ∈ σ  which is charged at  even if  ⪰  +  .Then the admissibility of σ  follows.(2) The proof for the other direction    N   0 works analogously as above.

1 :
procedure Reorder( 0 ](1)  ](2) . . . ](ℎ) ], ) ).Note that by merging   1 and   2 we obtain a stoichiometric matrix Γ  having redundant reactions (e.g., ( 1 ,   1 ), ( 1 ,   2 ) resulting in ( 1 ,    1 ), ( 1 ,    1 )) and zero reaction vectors (i.e., self-loops on    1 ), but they are omitted in (b).A directed cycle on which the chosen reaction   1  2 lies is depicted in gray.For  =  − 1 we assume that the reachability relation  0  N   holds.(c)  =  We have two different cases with respect to the existence of directed cycles.If the reaction network has no directed cycle, then the reachability relation  0 N   is guaranteed by Lemma 14. Assume that the reaction network contains at least one directed cycle   =  ](1) . . . ](ℎ) 0  N   ⇒  0 + Γ =   It follows from the definition of reachability.(2)0+ Γ =   ⇒  0  N   Since in the initial state  0 the number of each catalyzer molecule is higher than or equal to 1 and there is no reaction in N consuming a catalyzer species, it follows that for each state reachable from  0 the number of each catalyzer molecule is higher than or equal to 1. Let us remove all the catalyzer species of N from the reactions where they act as a catalyzer; i.e., for each  ∈ R of the form  =  +  1 →  +  2 we erase the catalyzer  to obtain  =  1 →  2 .In such a way a d-CRN N  is obtained so that for each  ⪰ ,  0  N   iff  0  N .N  satisfies the conditions of Proposition 15; hence the reachability relation  0  N    holds implying that  0  N   also holds.
. Then the reachability relation  0  N   holds if and only if there exists a vector  ∈ Z  ≥0 satisfying the state equation  0 + Γ =   .Proof.By changing the sign of the entries in the stoichiometric matrix Γ, we obtain a subconservative d-CRN N  of stoichiometric matrix −Γ satisfying the conditions of Proposition 17.We can consider the reachability problem   ? N   0 .