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

Authors: Gergely Szlobodnyik; Gábor Szederkényi
Journal: Complexity
 (2019)

---

## Abstract

In this paper we study the reachability problem of sub- and 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 sub- and 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.


---

## Body

## 1. 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–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–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–15].In the case of discrete state reaction networks the so-called 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–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 Propositions15 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 Section2 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.


## 2. Notations and Mathematical Background

In Table1 we summarize the notations and concepts of discrete chemical reaction networks which will be extensively used later.Table 1
Notations.
R the set of real numbers Z the set of integer numbers Z≥0 the set of non-negative integer numbers Tn×m the set of(n×m)-dimensional vectors over the set T 0n×m a zero matrix of dimensionn×m 1n×m a matrix of dimensionn×m for which all the entries are equal to 1 {0,1}n×m the set of(n×m)-dimensional binary vectors (all the entries are equal to 0 or 1) {-1,0,1}n×m the set of(n×m)-dimensional vectors composed of the entries -1,0,1 [A]i,: theith row of the matrix A a≺b fora,b∈Rn, ai<bi for i=1,…,n a⪯b fora,b∈Rn, ai≤bi for i=1,…,n σX an ordered sequence of states σr an ordered sequence of reaction vectors σS an ordered sequence of species σC an ordered sequence of complexes


### 2.1. Discrete State Chemical Reaction Networks

A discrete state Chemical Reaction Network (d-CRN) withnspecies, mcomplexes, and lreactions is a triple N=(S,C,R) so that (1)S=si∣i=1,…,nC=yj=∑i=1nαjisi∣si∈S,αji∈Z≥0,i=1,…,n,j=1,…,mR=rv=ysourcerv→yproductrv∣ysourcerv,yproductrv∈C,v=1,…,lwhere si is the i’thspecies, yj is the j’thcomplex, and rv is the v’threaction of the network. Moreover, αji is thestoichiometric coefficient of the i’th species in the j’th complex. For areaction  rv=ysource(rv)→yproduct(rv) of R,  ysource(rv) and yproduct(rv) are the source complex and the product complex, respectively.For each complexyj∈C,j∈{1,…,m}, the stoichiometric coefficients of the species can be represented as a vector of the following form:(2)y¯j=αj1αj2…αjn⊤For each r∈R, a reaction vector rij∈Zn can be associated with the track of the net molecular count changes of the species upon firing the reaction:(3)rij=y¯j-y¯iso that yj and yi are the corresponding source and product complexes of r. In the sequel the notation ri will be used for denoting both the i’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 r1,r2,…,rl is fixed and l=|R|.A d-CRN can also be represented by a directed graphG=G(V,E) such that the vertices and edges correspond to the complexes and the reactions, respectively, i.e.,(4)V=C(5)E=R.The direction of the edges is determined by the reactions of R, so that if yi→yj∈R, then there exists an edge e∈E from the vertex representing yi to the vertex of yj. 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-CRNN=(S,C,R). Thestoichiometric matrix  Γ∈Zn×l of N is defined as follows:(6)Γ=r1…rlNote that[Γ]ij encodes the net molecule count change on species si upon the occurrence of reaction rj. Besides Γ we also define the following matrices:(7)Γ+=y¯r1+…y¯rl+⊤(8)Γ-=y¯r1-…y¯rl-⊤where y¯ri+ denotes the vector form of the product complex belonging to reaction ri while y¯rj- represents the vector of the source complex associated with reaction rj. The relation among the above defined matrices is as follows: (9)Γ=Γ+-Γ-Example 2.
Let us consider the d-CRNN=(S,C,R) depicted in Figure 1. N characterizes a simple network of a bifunctional enzyme E having both phosphorylation and dephosphorylation activities on species I and Ip, respectively. This network is characterized by the following sets. (10)S=I,Ip,E,EI,EIpC=I+E,EI,Ip+E,EIpR=E+I→EI,IE→E+I,EI→Ip+E,E+Ip→EIp,EIp→E+Ip,EIp→E+IWe 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.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. [Γ]ij is the net change in the number of the i’th species upon occurring the j’th reaction.The molecular count of each species of a d-CRN at any timet≥0 is given by its state vector X(t)∈Z≥0n and the time evolution of the system is characterized by the followingdiscrete state equation:(11)Xt=X0+ΓNtwhere X(0) is the state vector belonging to the initial time instant and N(t)=[N1(t),N2(t),…,Nm(t)]⊤∈Z≥0l such that Nk(t)∈Z≥0 stores the number of occurrences of the k’th reaction up to time t. We note that N(t) 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 timet in the formulas.Definition 3.
Let us consider a d-CRNN=(S,C,R). It is said that:(1)
a speciess∈S is acatalyzer of a reaction r∈R if it has the form of r=s+s1→s+s2 with s1,s2∈S,(2)
a complexy∈C ischarged at state X∈Z≥0n if X⪰y,(3)
a reactionr∈R is charged if its respective source complex is charged,(4)
a stateX∈Z≥0nreacts to a state X′∈Z≥0n (denoted by X→X′) if there exists a reaction r∈R such that r is charged at state X and X+r=X′,(5)
areaction (vector) sequence  σr is an ordered set of reaction vectors σr=r1…rv where ri∈R,  i=1,…,v,(6)
astate transition sequence  σX is an ordered set states X0,X1,…,Xp so that X1→X2→…→Xp-1→Xp,(7)
a stateX′∈Z≥0n isreachable from a state X∈Z≥0n  (denoted by X⇝NX′ ) if there exists a directed path in the state space so that X=Xν(1)→Xν(2)→⋯→Xν(v)=X′.Considering a state transition sequenceσX=X0X1…Xp-1Xp, we call X0 and Xp the initial and target states, respectively, while Xi for i∈{1,…,p-1} are called transition states of σX.The condition that a reactionr∈R is charged at state X∈Z≥0n can be expressed by the inequality X⪰y¯r-. For a reaction sequence σr a state transition sequence σX=X0X1…Xv can be uniquely associated so that(12)Xj=Xj-1+rj,j∈1,…,vwhere the initial state X0 is assumed to be given. A state transition sequence σX is said to beadmissible if Xi⪰ri+1 for Xi∈σX,  i∈{0,…,v-1}; moreover, we say that a reaction sequence σr is admissible if the corresponding state transition sequence is admissible.From the reachability of a stateX′∈Z≥0n from an initial state X0∈Z≥0n, it follows that the following equation has a nonnegative integer solution c∈Z≥0l:(13)X′=X0+Γcwhere [c]i encodes the number of occurrences for reaction ri∈R for i∈{1,…,l}. However, it is important to note that from the existence of a nonnegative integer solution c of (13), the reachability relation X0⇝NX′ does not necessary follow.We note thatc of (13) corresponds to N(t) of (11). Since a solution c∈Z≥0l of (13) encodes the number of occurrences for each reaction in a fixed order, the following equality is fulfilled:(14)Γc=∑i=1hriwhere h=∑i=1l[c]i and ri∈R for i∈{1,…,h}. When we want to emphasize that a reaction vector sequence is encoded by a particular c∈Z≥0l, we will use the notation σrc=r1,…,rh and the state transition sequence determined by σrc will be denoted by σXc.Definition 4.
Let us consider a d-CRNN with stoichiometric matrix Γ∈Zn×l and an initial state X0∈Z≥0n. Thereachable state space  Reach(N,X0) of N with initial state X0 is the set of nonnegative discrete states reachable from X0.(15)ReachN,X0=X∣X∈Z≥0n,X0⇝NX


## 2.1. Discrete State Chemical Reaction Networks

A discrete state Chemical Reaction Network (d-CRN) withnspecies, mcomplexes, and lreactions is a triple N=(S,C,R) so that (1)S=si∣i=1,…,nC=yj=∑i=1nαjisi∣si∈S,αji∈Z≥0,i=1,…,n,j=1,…,mR=rv=ysourcerv→yproductrv∣ysourcerv,yproductrv∈C,v=1,…,lwhere si is the i’thspecies, yj is the j’thcomplex, and rv is the v’threaction of the network. Moreover, αji is thestoichiometric coefficient of the i’th species in the j’th complex. For areaction  rv=ysource(rv)→yproduct(rv) of R,  ysource(rv) and yproduct(rv) are the source complex and the product complex, respectively.For each complexyj∈C,j∈{1,…,m}, the stoichiometric coefficients of the species can be represented as a vector of the following form:(2)y¯j=αj1αj2…αjn⊤For each r∈R, a reaction vector rij∈Zn can be associated with the track of the net molecular count changes of the species upon firing the reaction:(3)rij=y¯j-y¯iso that yj and yi are the corresponding source and product complexes of r. In the sequel the notation ri will be used for denoting both the i’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 r1,r2,…,rl is fixed and l=|R|.A d-CRN can also be represented by a directed graphG=G(V,E) such that the vertices and edges correspond to the complexes and the reactions, respectively, i.e.,(4)V=C(5)E=R.The direction of the edges is determined by the reactions of R, so that if yi→yj∈R, then there exists an edge e∈E from the vertex representing yi to the vertex of yj. 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-CRNN=(S,C,R). Thestoichiometric matrix  Γ∈Zn×l of N is defined as follows:(6)Γ=r1…rlNote that[Γ]ij encodes the net molecule count change on species si upon the occurrence of reaction rj. Besides Γ we also define the following matrices:(7)Γ+=y¯r1+…y¯rl+⊤(8)Γ-=y¯r1-…y¯rl-⊤where y¯ri+ denotes the vector form of the product complex belonging to reaction ri while y¯rj- represents the vector of the source complex associated with reaction rj. The relation among the above defined matrices is as follows: (9)Γ=Γ+-Γ-Example 2.
Let us consider the d-CRNN=(S,C,R) depicted in Figure 1. N characterizes a simple network of a bifunctional enzyme E having both phosphorylation and dephosphorylation activities on species I and Ip, respectively. This network is characterized by the following sets. (10)S=I,Ip,E,EI,EIpC=I+E,EI,Ip+E,EIpR=E+I→EI,IE→E+I,EI→Ip+E,E+Ip→EIp,EIp→E+Ip,EIp→E+IWe 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.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. [Γ]ij is the net change in the number of the i’th species upon occurring the j’th reaction.The molecular count of each species of a d-CRN at any timet≥0 is given by its state vector X(t)∈Z≥0n and the time evolution of the system is characterized by the followingdiscrete state equation:(11)Xt=X0+ΓNtwhere X(0) is the state vector belonging to the initial time instant and N(t)=[N1(t),N2(t),…,Nm(t)]⊤∈Z≥0l such that Nk(t)∈Z≥0 stores the number of occurrences of the k’th reaction up to time t. We note that N(t) 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 timet in the formulas.Definition 3.
Let us consider a d-CRNN=(S,C,R). It is said that:(1)
a speciess∈S is acatalyzer of a reaction r∈R if it has the form of r=s+s1→s+s2 with s1,s2∈S,(2)
a complexy∈C ischarged at state X∈Z≥0n if X⪰y,(3)
a reactionr∈R is charged if its respective source complex is charged,(4)
a stateX∈Z≥0nreacts to a state X′∈Z≥0n (denoted by X→X′) if there exists a reaction r∈R such that r is charged at state X and X+r=X′,(5)
areaction (vector) sequence  σr is an ordered set of reaction vectors σr=r1…rv where ri∈R,  i=1,…,v,(6)
astate transition sequence  σX is an ordered set states X0,X1,…,Xp so that X1→X2→…→Xp-1→Xp,(7)
a stateX′∈Z≥0n isreachable from a state X∈Z≥0n  (denoted by X⇝NX′ ) if there exists a directed path in the state space so that X=Xν(1)→Xν(2)→⋯→Xν(v)=X′.Considering a state transition sequenceσX=X0X1…Xp-1Xp, we call X0 and Xp the initial and target states, respectively, while Xi for i∈{1,…,p-1} are called transition states of σX.The condition that a reactionr∈R is charged at state X∈Z≥0n can be expressed by the inequality X⪰y¯r-. For a reaction sequence σr a state transition sequence σX=X0X1…Xv can be uniquely associated so that(12)Xj=Xj-1+rj,j∈1,…,vwhere the initial state X0 is assumed to be given. A state transition sequence σX is said to beadmissible if Xi⪰ri+1 for Xi∈σX,  i∈{0,…,v-1}; moreover, we say that a reaction sequence σr is admissible if the corresponding state transition sequence is admissible.From the reachability of a stateX′∈Z≥0n from an initial state X0∈Z≥0n, it follows that the following equation has a nonnegative integer solution c∈Z≥0l:(13)X′=X0+Γcwhere [c]i encodes the number of occurrences for reaction ri∈R for i∈{1,…,l}. However, it is important to note that from the existence of a nonnegative integer solution c of (13), the reachability relation X0⇝NX′ does not necessary follow.We note thatc of (13) corresponds to N(t) of (11). Since a solution c∈Z≥0l of (13) encodes the number of occurrences for each reaction in a fixed order, the following equality is fulfilled:(14)Γc=∑i=1hriwhere h=∑i=1l[c]i and ri∈R for i∈{1,…,h}. When we want to emphasize that a reaction vector sequence is encoded by a particular c∈Z≥0l, we will use the notation σrc=r1,…,rh and the state transition sequence determined by σrc will be denoted by σXc.Definition 4.
Let us consider a d-CRNN with stoichiometric matrix Γ∈Zn×l and an initial state X0∈Z≥0n. Thereachable state space  Reach(N,X0) of N with initial state X0 is the set of nonnegative discrete states reachable from X0.(15)ReachN,X0=X∣X∈Z≥0n,X0⇝NX


## 3. 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:(16)ILPminx⁡a⊤xsubjecttoAx≤bx∈Znwhere x is the n-dimensional vector of decision variables while a∈Zn,  A∈Zm×n, and b∈Zm 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 ax∈Zn vector satisfying the set of specified constraints, then the problem is called ILP feasibility problem.(17)FPP=x∣Ax≤b,A∈Zm×n,b∈Zm,x∈RnP∩Zn=?∅An ILP feasibility problem, as a decision problem, addresses the question of whether the polytope P contains an integer lattice point, formally P∩Zn=?∅. 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 A and b. Therefore, a feasibility problem of the form (17), assuming fixed dimension n, can be decided in polynomial time in the number of constraints m and the maximum of the absolute values of the coefficients A and b by means of the Lenstra algorithm [28, 29]. Moreover, the number of integer lattice points in P can also be numerated in polynomial time in m and the maximum of the absolute value of the coefficients using Barvinok’s integer lattice point counting algorithm [30–33]. We note that for the Barvinok algorithm there exists an effective implementation called LattE [34].


## 4. Sub- and Superconservative d-CRNs

We define conservativity and subconservativity in the same way as they were introduced, e.g. in [4, 16].Definition 5.
A d-CRNN=(S,C,R) having stoichiometric matrix Γ∈Zn×l is calledsubconservative (superconservative) if there exists a strictly positive vector z∈R>0n for which z⊤Γ≤01×l (z⊤Γ≥01×l) holds. The vector z is called aconservation vector.An important property related to subconservativity is the strong boundedness which is defined as follows.Definition 6.
A d-CRNN is said to bestrongly bounded if, for any X0∈Z≥0n initial state, the reachable state space Reach(N,X0) is bounded.The subconservative property of the reaction network structure is a necessary and sufficient condition of strong boundedness [16, 35].Proposition 7 (see [35]).
Let us consider a d-CRNN. The following propositions are equivalent: (1)
N is subconservative,(2)
N is strongly bounded.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-CRNN=(S,C,R) with stoichiometric matrix Γ∈Zn×l. The d-CRN N is said to beconservative if there exists a vector z∈R>0n satisfying the matrix equation z⊤Γ=01×l.We note that the above structural properties can be easily decided in polynomial time by means of an LP of the following form:(18)min∑j=1nzjs.t.z⊤Γ≤01×lorz⊤Γ≥01×lz⪰0n×1+εn×1,ε≻0n×1The relationship between sub- and superconservativity can be expressed by the following proposition.Proposition 9.
A d-CRNN with stoichiometric matrix ΓN∈Zn×l is subconservative if and only if the d-CRN N′ with stoichiometric matrix ΓN′=-ΓN is superconservative.Proof.
(19) z ⊤ Γ ≤ 0 1 × m ⇔ z ⊤ - Γ ≥ 0 1 × lWe note that-ΓN means the change of the direction of each reaction in the d-CRN N of stoichiometric matrix ΓN.Example 10.
Figure2 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.Figure 2
A pair of sub- and superconservative reaction network structures denoted byN 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) (b)From Proposition9 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.Proposition 11.
Let us consider a subconservative d-CRNN characterized by the matrices ΓN=Γ, ΓN+=Γ+ and a superconservative d-CRN N′ with matrices ΓN′=-Γ, ΓN′-=Γ+. Let us take an initial state X0∈Z≥0n and a target state X′∈Z≥0n. Then the reachability X0⇝NX′ holds if and only if X′⇝N′X0 also holds.Proof.
(1)
X 0 ⇝ N X ′ ⇒ ∃ c ∈ Z ≥ 0 l such that X0+Γc=X′ which is equivalent to X′+(-Γ)c=X0. FromX0⇝NX′ it follows that the solution c∈Z≥0l can be decomposed to an admissible reaction vector sequence σrc=r1c…rhc, h=∑i=1l[c]i; i.e., all the states of σXc determined by σrc are composed of nonnegative entries. Then, by reversing σXc, we obtain a nonnegative state transition sequence σ^Xc from X′ to X0 which is uniquely determined by means of the reaction vector sequence σ^rc=-rhc…-r1c. It is also needed to show thatσ^rc is an admissible reaction sequence. This can be done as follows: for each state X∈σXc\X0 there exists a reaction r∈σrc so that upon firing r the resulting state is X, from which it follows that X⪰y¯r+; moreover, considering the reversed reaction sequence σ^rc, the reaction vector that will occur at state X is -r∈σ^rc which is charged at X even if X⪰y¯r+. Then the admissibility ofσ^rc follows. (2)
The proof for the other directionX′⇝N′X0 works analogously as above.The importance of Proposition11 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.


## 5. Reachability Analysis

### 5.1. Low-Dimensional Case

In this section the case of low-dimensional (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 stateM=M(Γ-) as follows:(20)MΓ-i=max⁡Γ-ij:j=1,…,li=1,…,n.Here Γ- is defined by (8). Note that the set X∣X∈Z≥0n,X⪰M contains all the states where each reaction is charged.Proposition 12.
Let us consider a subconservative d-CRNN with stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l. Assume that rank(Γ)≤2. We consider an initial state X0∈Z≥0n and a target state X′∈Z≥0n such that X0⪰M and X′⪰M hold where M=M(Γ-) is defined by (20). Then the state X′ is reachable from X0 through a state transition sequence σX=X0X1…X′ for which, ∀X∈σX,X⪰M if and only if the equation(21)Γc=X′-X0has a nonnegative integer solution c.Proof.
(1)
IfX′ is reachable from X0 through an admissible state transition sequence σX, then it follows that a solution c∈Z≥0l exists. (2)
Assume that there existsc∈Z≥0n such that X0+Γc=X′ holds. Let us consider any reaction vector decomposition σr=rν(1)…rν(h) of c where ∑j=1hrν(j)=c and ∑j=1l[c]j=h. We show that Algorithm 1 returns a permutation of σr so that for all the transition states X the inequality X⪰M holds. Let us assume that there exists a transition stateXi,Xi⪰M, so that the forthcoming state Xi+1 satisfies the inequality [Xi+1]d<[M]d for some d∈{1,2}. For the target state X′ to be reached the inequality X′⪰M holds; hence there exists a reaction increasing the state variable along the coordinate d. Let us assume that all the reactions increasing the state variable along Xi decrease the other coordinate d′ so that the resulting forthcoming state Xi+1 satisfies the inequality [Xi+1]d′<[M]d′. Then Xi=M holds. Now there are two different cases: ( P 1 )  If X′=M, then Algorithm 1 terminates, and the correctness follows. ( P 2 )  If X′≠M, then the subconservativity of N implies that it is not possible to reach a state X, X⪰M, X≠M; i.e., X′ is not reachable from Xi. This is contradiction, since arbitrary permutation of the initial ordering σr results in the same target state X′, given the initial state X0. Then the correctness of Algorithm 1 follows.Algorithm 1
1:procedure  Reorder(X0rν(1)rν(2)…rν(h), M) 2:Xcurrent←X0 3:for  i=1 to h  do 4:if  Xcurrent=X′  then 5:return   rν1rν2…rνh 6:end if 7:if[Xcurrent+rν(i)]l<[M]l  for some  l∈{1,…,n}  then 8:Choose a transition vectorrν(j),  i<j≤h  so that X c u r r e n t + r ν ( j ) ⪰ M 9:r′←rν(i) 10:rν(i)←rν(j) 11:rν(j)←r′ 12:end if 13:Xcurrent←Xcurrent+rν(i) 14:end for 15:return   rν(1)rν(2)…rν(h) 16:end procedureAlgorithm1 can be easily extended to the class of superconservative reaction networks.Corollary 13.
Let us consider a superconservative d-CRNN with stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈0,1n×l. Assume that rank(Γ)≤2 holds and consider an initial state X0∈Z≥0n and a target state X′∈Z≥0n for which X0⪰M and X′⪰M hold where M is defined by (20). Then the state X′∈Z≥0n is reachable from X0 if and only if the equation(22)Γc=X′-X0has a nonnegative integer solution c.Proof.
According to Proposition11 we can consider a subconservative d-CRN N′ of stoichiometric matrix -Γ and take the reachability problem X′⇝?N′X0. Then Proposition 7 can be applied.


### 5.2. 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 relationX0⇝NX′ is equivalent to the existence of a c∈Z≥0l solution of the d-CRN state equation X0+Γc=X′. 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.Lemma 14.
Let us consider a d-CRNN with stoichiometric matrix Γ∈{-1,0,1}n×l such that for all r∈R reactions ∑i=1n[y¯+]i≤1 and ∑i=1n[y¯-]i=1 holds. Assume that the reaction network of N does not contain directed cycle (i.e., N has an acyclic network structure). Consider two states X0,X′∈Z≥0n. Then the reachability relation X0⇝NX′ holds if and only if there exists c∈Z≥0l vector satisfying the state equation X0+Γc=X′.Now we can state the result on the reachability of subconservative d-CRNs composed of reaction having at most one input and one output species.Proposition 15.
Let us consider a subconservative d-CRNN=(S,C,R) of stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l for which C=S∪{0}. Assume that for all r∈R reactions ∑i=1n[y¯+]i≤1 and ∑i=1n[y¯-]i=1 hold. Let us consider two states X0,X′∈Z≥0n so that X0⪰M and X′⪰M hold where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.Proof.
(1)
X 0 ⇝ N X ′ ⇒ ∃ c ∈ Z ≥ 0 l : X 0 + Γ c = X ′ By the definition of reachability it is guaranteed that the state equation is satisfied with somec∈Z≥0l. (2)
X 0 ⇝ N X ′ ⇐ ∃ c ∈ Z ≥ 0 l : X 0 + Γ c = X ′ For this side an inductive proof is employed.
(a)
k = 2 If a d-CRN is 2-dimensional, according to Proposition12, the existence of a solution c∈Z≥0l of the state equation implies that the reachability relation holds. (b)
Inductive assumption Fork=n-1 we assume that the reachability relation X0⇝NX′ holds. (c)
k = n We have two different cases with respect to the existence of directed cycles. If the reaction network has no directed cycle, then the reachability relationX0⇝NX′ is guaranteed by Lemma 14. Assume that the reaction network contains at least one directed cycle(23)σS=sν1…sνh whereh≤n,sν(1)=sν(h) and sν(i)≠sν(j) for i,j∈{1,…,h},i≠j. Note again that C=S∪{0}, and hence σS can be considered as a directed cycle of complexes in the reaction network (i.e., σS=σC=sν(1)…sν(h)). Let us consider an arbitrary rk1k2∈R reaction defined between some sk1,sk2∈σS, i.e. rk1k2=sk1→sk2. Now we construct a d-CRNN′=(S′,C′,R′) from the stoichiometric matrix Γ∈{-1,0,1}(n-1)×l and Γ-∈{0,1}(n-1)×l as follows:(24)Γ′i,:=Γi,:,i<kmax,i≠kmin,Γkmin,:+Γkmax,:,i=kmin,Γi+1,:,kmax≤i≤n-1, and(25)Γ-′i,:=Γ-i,:,i<kmax,i≠kmin,Γ-kmin,:+Γ-kmax,:,i=kmin,Γ-i+1,:,kmax≤i≤n-1. Herekmin=min⁡{k1,k2} and kmax=max⁡{k1,k2}. 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′∈R′ the ordered pair of source complex and product complex of r∈R from which it is obtained. In such a way every reaction of N′ is uniquely described by an ordered pair (r′,r)∈R′×R. Then by the mapping P((r′,r))=r one can uniquely determine the reaction r∈R from which r′∈R′ is derived. Let us construct the statesX0m∈Z≥0n-1 and Xm′∈Z≥0n-1 as follows:(26)X0mi=X0i,i<kmax,i≠kmin,X0kmin+X0kmax,i=kmin,X0i+1,kmax≤i≤n-1.(27)X′mi=X′i,i<kmax,i≠kmin,X′kmin+X′kmax,i=kmin,X′i+1,kmax≤i≤n-1. Then we have thatX0m⪰M(Γ-′) and X′m⪰M(Γ-′), and hence the (n-1)-dimensional d-CRN N′ with the initial and final states X0m and X′m satisfies the assumptions of the proposition. From X0+Γc=X′ we have that X0m+Γ′c=X′m holds; hence, according to the (n-1)-dimensional inductive assumption, the reachability relation(28)X0m⇝N′X′m follows. Let us consider an admissible reaction vector sequenceσr′ associated with relation (28). Since for each r′∈R′ we associated the reaction r∈R from which r′ is obtained, making use of the mapping P:R′×R→R, we can consider the reaction vector sequence σr  (r∈R∀r∈σr) uniquely determined by σr′. We start from X0 and modify the state variable X∈Z≥0n according to the reaction vector sequence σr. We may get to two invalid cases: ( C 1 )  [X]k2=0, but the source complex of the forthcoming reaction rcurrent∈σr is sk2. Then, according to the (n-1)-dimensional reachability, it is guaranteed that sk1 is charged at the current state X. Let us insert rk1k2 into σr before the current reaction rcurrent. ( C 2 )  [X]k1=0, but the source complex of the forthcoming reaction rcurrent∈σr is sk1. Then, according to the (n-1)-dimensional reachability, it is guaranteed that sk2 is charged at the current state X. It is known that sk1 can be reached from sk2 along a reaction vector sequence σr∗ in the reaction network of N. Let us insert σr∗ into σr before the current reaction rcurrent. By modifyingσr according to the above discussed cases (C1) and (C2), we obtain an admissible reaction vector sequence σrmod with respect to the reachability relation(29)X0⇝N′X∗ whereX∗⪰0n, [X∗]i=[X′]i for i∈{1,…n}, i≠k1 and i≠k2; moreover, [X∗]k1+[X∗]k2=X′. According to the assumptions N contains directed paths both from sk1 to sk2 and from sk2 to sk1; hence the reachability relation X∗⇝NX′ follows. Then, due to the transitivity of the relation ⇝N, we have that X0⇝NX′ also holds.Figure 3
Graphical explanation of how the reaction network structure ofN′ in the proof of Proposition 15 is constructed. (a) Reaction network structure of an n-dimensional d-CRN N. (b) Reaction network structure of N′ resulting from merging the species sk1 and sk2 of N along their shared reaction rk1k2 (and reverse counterpart reaction rk2k1). Note that by merging sk1 and sk2 we obtain a stoichiometric matrix Γ′ having redundant reactions (e.g., (s1,sk1), (s1,sk2) resulting in (s1,sk1′), (s1,sk1′)) and zero reaction vectors (i.e., self-loops on sk1′), but they are omitted in (b). A directed cycle on which the chosen reaction rk1k2 lies is depicted in gray.
(a) (b)Proposition15 can be extended to the case of superconservative d-CRNs.Corollary 16.
Let us consider a superconservative d-CRNN=(S,C,R) with stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l for which C=S. Assume that for all r∈R reactions ∑i=1n[y¯+]i=1 and ∑i=1n[y¯-]i≤1 hold. Let us consider two states X0,X′∈Z≥0n so that X0⪰M and X′⪰M hold where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.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 X′⇝?N′X0.We can extend Proposition15 by allowing the restricted application of catalyzer species as follows.Proposition 17.
Let us consider a subconservative d-CRNN=(S,C,R) of stoichiometric matrix Γ∈{-1,0,1}n×l and  Γ-∈{0,1}n×l. Assume that for each reaction r: (1)
r=s1→s2 for some s1,s2∈S, s1≠s2, s1≠0, OR(2)
r=s+s1→s+s2 where s,s1,s2∈S, s≠s1≠s2, s≠0, s1≠0 and ∀r′∈Rr′ does not consume s.
Let us consider two statesX0,X′∈Z≥0n for which X0⪰M and X′⪰M where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.Proof.
(1)
X 0 ⇝ N X ′ ⇒ X 0 + Γ c = X ′ It follows from the definition of reachability. (2)
X 0 + Γ c = X ′ ⇒ X 0 ⇝ N X ′ Since in the initial stateX0 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 X0 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∈R of the form r=s+s1→s+s2 we erase the catalyzer s to obtain r′=s1→s2. In such a way a d-CRN N′ is obtained so that for each X⪰M,  X0⇝N′X iff X0⇝NX. N′ satisfies the conditions of Proposition 15; hence the reachability relation X0⇝N′X′ holds implying that X0⇝NX′ also holds.According to the duality of the sub- and superconservativity properties, we can extend Proposition17 to the case of superconservative d-CRNs.Corollary 18.
Let us consider a superconservative d-CRNN=(S,C,R) of stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l. Assume that for each reaction r(1)
r=s1→s2 for some s1,s2∈S, s1≠s2, s2≠0, OR(2)
r=s+s1→s+s2 where s,s1,s2∈S, s≠s1≠s2, s≠0, s2≠0 and ∀r′∈Rr′ does not produce s.
Let us consider two statesX0,X′∈Z≥0n for which X0⪰M and X′⪰M where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.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 X′⇝?N′X0.By the above corollary, any reachability problem on a superconservative d-CRN satisfying the conditions of Corollary18 can be easily traced back to that of a subconservative network; hence the problem is equivalent to finding a c∈Z≥0l solution for the respective d-CRN state equation.The reaction network class covered by the above statements might be beneficial in modeling first- and second-order (bio)chemical reaction networks. For a representative example, see Example19 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 Figure4. 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 z=156 and Proposition 15 can be applied.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. aj,Aj and ij,Ij for j=0,…13 denote the active and inactive proteins, respectively. Lower indices denote the number of phosphorylated residues.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.Remark 20.
Let us consider a subconservative (superconservative) d-CRNN=(S,C,R) of n species, m complexes, and l reactions. Assume that N satisfies the conditions of Proposition 17 (Corollary 18). Then for any X0,X′∈Z≥0n initial and target states for which X0⪰M(Γ-),X′⪰M(Γ-) hold we have that the number of distinct trajectories σX satisfying the reachability relation X0⇝?NX′ can be determined in polynomial time in the distance of X0 and X′, given the fixed number of reactions l in the d-CRN. The explanation of this is the following. According to Proposition 17 (Corollary 18) the reachability problem X0⇝?NX′ is equivalent to the existence of a nonnegative integer solution c∈Z≥0l of the state equation X0+Γc=X′. In this way the reachability problem can be reformulated as an ILP feasibility problem in terms of c, 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 reactionsl,(2)
polynomial in the number of constraints, that is, in the number of speciesn,(3)
polynomial in the maximum of the absolute values of the coefficientsΓ, X′-X0.The particular importance of Remark20 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 X0 initial state.


## 5.1. Low-Dimensional Case

In this section the case of low-dimensional (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 stateM=M(Γ-) as follows:(20)MΓ-i=max⁡Γ-ij:j=1,…,li=1,…,n.Here Γ- is defined by (8). Note that the set X∣X∈Z≥0n,X⪰M contains all the states where each reaction is charged.Proposition 12.
Let us consider a subconservative d-CRNN with stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l. Assume that rank(Γ)≤2. We consider an initial state X0∈Z≥0n and a target state X′∈Z≥0n such that X0⪰M and X′⪰M hold where M=M(Γ-) is defined by (20). Then the state X′ is reachable from X0 through a state transition sequence σX=X0X1…X′ for which, ∀X∈σX,X⪰M if and only if the equation(21)Γc=X′-X0has a nonnegative integer solution c.Proof.
(1)
IfX′ is reachable from X0 through an admissible state transition sequence σX, then it follows that a solution c∈Z≥0l exists. (2)
Assume that there existsc∈Z≥0n such that X0+Γc=X′ holds. Let us consider any reaction vector decomposition σr=rν(1)…rν(h) of c where ∑j=1hrν(j)=c and ∑j=1l[c]j=h. We show that Algorithm 1 returns a permutation of σr so that for all the transition states X the inequality X⪰M holds. Let us assume that there exists a transition stateXi,Xi⪰M, so that the forthcoming state Xi+1 satisfies the inequality [Xi+1]d<[M]d for some d∈{1,2}. For the target state X′ to be reached the inequality X′⪰M holds; hence there exists a reaction increasing the state variable along the coordinate d. Let us assume that all the reactions increasing the state variable along Xi decrease the other coordinate d′ so that the resulting forthcoming state Xi+1 satisfies the inequality [Xi+1]d′<[M]d′. Then Xi=M holds. Now there are two different cases: ( P 1 )  If X′=M, then Algorithm 1 terminates, and the correctness follows. ( P 2 )  If X′≠M, then the subconservativity of N implies that it is not possible to reach a state X, X⪰M, X≠M; i.e., X′ is not reachable from Xi. This is contradiction, since arbitrary permutation of the initial ordering σr results in the same target state X′, given the initial state X0. Then the correctness of Algorithm 1 follows.Algorithm 1
1:procedure  Reorder(X0rν(1)rν(2)…rν(h), M) 2:Xcurrent←X0 3:for  i=1 to h  do 4:if  Xcurrent=X′  then 5:return   rν1rν2…rνh 6:end if 7:if[Xcurrent+rν(i)]l<[M]l  for some  l∈{1,…,n}  then 8:Choose a transition vectorrν(j),  i<j≤h  so that X c u r r e n t + r ν ( j ) ⪰ M 9:r′←rν(i) 10:rν(i)←rν(j) 11:rν(j)←r′ 12:end if 13:Xcurrent←Xcurrent+rν(i) 14:end for 15:return   rν(1)rν(2)…rν(h) 16:end procedureAlgorithm1 can be easily extended to the class of superconservative reaction networks.Corollary 13.
Let us consider a superconservative d-CRNN with stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈0,1n×l. Assume that rank(Γ)≤2 holds and consider an initial state X0∈Z≥0n and a target state X′∈Z≥0n for which X0⪰M and X′⪰M hold where M is defined by (20). Then the state X′∈Z≥0n is reachable from X0 if and only if the equation(22)Γc=X′-X0has a nonnegative integer solution c.Proof.
According to Proposition11 we can consider a subconservative d-CRN N′ of stoichiometric matrix -Γ and take the reachability problem X′⇝?N′X0. Then Proposition 7 can be applied.


## 5.2. 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 relationX0⇝NX′ is equivalent to the existence of a c∈Z≥0l solution of the d-CRN state equation X0+Γc=X′. 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.Lemma 14.
Let us consider a d-CRNN with stoichiometric matrix Γ∈{-1,0,1}n×l such that for all r∈R reactions ∑i=1n[y¯+]i≤1 and ∑i=1n[y¯-]i=1 holds. Assume that the reaction network of N does not contain directed cycle (i.e., N has an acyclic network structure). Consider two states X0,X′∈Z≥0n. Then the reachability relation X0⇝NX′ holds if and only if there exists c∈Z≥0l vector satisfying the state equation X0+Γc=X′.Now we can state the result on the reachability of subconservative d-CRNs composed of reaction having at most one input and one output species.Proposition 15.
Let us consider a subconservative d-CRNN=(S,C,R) of stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l for which C=S∪{0}. Assume that for all r∈R reactions ∑i=1n[y¯+]i≤1 and ∑i=1n[y¯-]i=1 hold. Let us consider two states X0,X′∈Z≥0n so that X0⪰M and X′⪰M hold where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.Proof.
(1)
X 0 ⇝ N X ′ ⇒ ∃ c ∈ Z ≥ 0 l : X 0 + Γ c = X ′ By the definition of reachability it is guaranteed that the state equation is satisfied with somec∈Z≥0l. (2)
X 0 ⇝ N X ′ ⇐ ∃ c ∈ Z ≥ 0 l : X 0 + Γ c = X ′ For this side an inductive proof is employed.
(a)
k = 2 If a d-CRN is 2-dimensional, according to Proposition12, the existence of a solution c∈Z≥0l of the state equation implies that the reachability relation holds. (b)
Inductive assumption Fork=n-1 we assume that the reachability relation X0⇝NX′ holds. (c)
k = n We have two different cases with respect to the existence of directed cycles. If the reaction network has no directed cycle, then the reachability relationX0⇝NX′ is guaranteed by Lemma 14. Assume that the reaction network contains at least one directed cycle(23)σS=sν1…sνh whereh≤n,sν(1)=sν(h) and sν(i)≠sν(j) for i,j∈{1,…,h},i≠j. Note again that C=S∪{0}, and hence σS can be considered as a directed cycle of complexes in the reaction network (i.e., σS=σC=sν(1)…sν(h)). Let us consider an arbitrary rk1k2∈R reaction defined between some sk1,sk2∈σS, i.e. rk1k2=sk1→sk2. Now we construct a d-CRNN′=(S′,C′,R′) from the stoichiometric matrix Γ∈{-1,0,1}(n-1)×l and Γ-∈{0,1}(n-1)×l as follows:(24)Γ′i,:=Γi,:,i<kmax,i≠kmin,Γkmin,:+Γkmax,:,i=kmin,Γi+1,:,kmax≤i≤n-1, and(25)Γ-′i,:=Γ-i,:,i<kmax,i≠kmin,Γ-kmin,:+Γ-kmax,:,i=kmin,Γ-i+1,:,kmax≤i≤n-1. Herekmin=min⁡{k1,k2} and kmax=max⁡{k1,k2}. 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′∈R′ the ordered pair of source complex and product complex of r∈R from which it is obtained. In such a way every reaction of N′ is uniquely described by an ordered pair (r′,r)∈R′×R. Then by the mapping P((r′,r))=r one can uniquely determine the reaction r∈R from which r′∈R′ is derived. Let us construct the statesX0m∈Z≥0n-1 and Xm′∈Z≥0n-1 as follows:(26)X0mi=X0i,i<kmax,i≠kmin,X0kmin+X0kmax,i=kmin,X0i+1,kmax≤i≤n-1.(27)X′mi=X′i,i<kmax,i≠kmin,X′kmin+X′kmax,i=kmin,X′i+1,kmax≤i≤n-1. Then we have thatX0m⪰M(Γ-′) and X′m⪰M(Γ-′), and hence the (n-1)-dimensional d-CRN N′ with the initial and final states X0m and X′m satisfies the assumptions of the proposition. From X0+Γc=X′ we have that X0m+Γ′c=X′m holds; hence, according to the (n-1)-dimensional inductive assumption, the reachability relation(28)X0m⇝N′X′m follows. Let us consider an admissible reaction vector sequenceσr′ associated with relation (28). Since for each r′∈R′ we associated the reaction r∈R from which r′ is obtained, making use of the mapping P:R′×R→R, we can consider the reaction vector sequence σr  (r∈R∀r∈σr) uniquely determined by σr′. We start from X0 and modify the state variable X∈Z≥0n according to the reaction vector sequence σr. We may get to two invalid cases: ( C 1 )  [X]k2=0, but the source complex of the forthcoming reaction rcurrent∈σr is sk2. Then, according to the (n-1)-dimensional reachability, it is guaranteed that sk1 is charged at the current state X. Let us insert rk1k2 into σr before the current reaction rcurrent. ( C 2 )  [X]k1=0, but the source complex of the forthcoming reaction rcurrent∈σr is sk1. Then, according to the (n-1)-dimensional reachability, it is guaranteed that sk2 is charged at the current state X. It is known that sk1 can be reached from sk2 along a reaction vector sequence σr∗ in the reaction network of N. Let us insert σr∗ into σr before the current reaction rcurrent. By modifyingσr according to the above discussed cases (C1) and (C2), we obtain an admissible reaction vector sequence σrmod with respect to the reachability relation(29)X0⇝N′X∗ whereX∗⪰0n, [X∗]i=[X′]i for i∈{1,…n}, i≠k1 and i≠k2; moreover, [X∗]k1+[X∗]k2=X′. According to the assumptions N contains directed paths both from sk1 to sk2 and from sk2 to sk1; hence the reachability relation X∗⇝NX′ follows. Then, due to the transitivity of the relation ⇝N, we have that X0⇝NX′ also holds.Figure 3
Graphical explanation of how the reaction network structure ofN′ in the proof of Proposition 15 is constructed. (a) Reaction network structure of an n-dimensional d-CRN N. (b) Reaction network structure of N′ resulting from merging the species sk1 and sk2 of N along their shared reaction rk1k2 (and reverse counterpart reaction rk2k1). Note that by merging sk1 and sk2 we obtain a stoichiometric matrix Γ′ having redundant reactions (e.g., (s1,sk1), (s1,sk2) resulting in (s1,sk1′), (s1,sk1′)) and zero reaction vectors (i.e., self-loops on sk1′), but they are omitted in (b). A directed cycle on which the chosen reaction rk1k2 lies is depicted in gray.
(a) (b)Proposition15 can be extended to the case of superconservative d-CRNs.Corollary 16.
Let us consider a superconservative d-CRNN=(S,C,R) with stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l for which C=S. Assume that for all r∈R reactions ∑i=1n[y¯+]i=1 and ∑i=1n[y¯-]i≤1 hold. Let us consider two states X0,X′∈Z≥0n so that X0⪰M and X′⪰M hold where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.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 X′⇝?N′X0.We can extend Proposition15 by allowing the restricted application of catalyzer species as follows.Proposition 17.
Let us consider a subconservative d-CRNN=(S,C,R) of stoichiometric matrix Γ∈{-1,0,1}n×l and  Γ-∈{0,1}n×l. Assume that for each reaction r: (1)
r=s1→s2 for some s1,s2∈S, s1≠s2, s1≠0, OR(2)
r=s+s1→s+s2 where s,s1,s2∈S, s≠s1≠s2, s≠0, s1≠0 and ∀r′∈Rr′ does not consume s.
Let us consider two statesX0,X′∈Z≥0n for which X0⪰M and X′⪰M where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.Proof.
(1)
X 0 ⇝ N X ′ ⇒ X 0 + Γ c = X ′ It follows from the definition of reachability. (2)
X 0 + Γ c = X ′ ⇒ X 0 ⇝ N X ′ Since in the initial stateX0 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 X0 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∈R of the form r=s+s1→s+s2 we erase the catalyzer s to obtain r′=s1→s2. In such a way a d-CRN N′ is obtained so that for each X⪰M,  X0⇝N′X iff X0⇝NX. N′ satisfies the conditions of Proposition 15; hence the reachability relation X0⇝N′X′ holds implying that X0⇝NX′ also holds.According to the duality of the sub- and superconservativity properties, we can extend Proposition17 to the case of superconservative d-CRNs.Corollary 18.
Let us consider a superconservative d-CRNN=(S,C,R) of stoichiometric matrix Γ∈{-1,0,1}n×l and Γ-∈{0,1}n×l. Assume that for each reaction r(1)
r=s1→s2 for some s1,s2∈S, s1≠s2, s2≠0, OR(2)
r=s+s1→s+s2 where s,s1,s2∈S, s≠s1≠s2, s≠0, s2≠0 and ∀r′∈Rr′ does not produce s.
Let us consider two statesX0,X′∈Z≥0n for which X0⪰M and X′⪰M where M=M(Γ-) is defined by (20). Then the reachability relation X0⇝NX′ holds if and only if there exists a vector c∈Z≥0l satisfying the state equation X0+Γc=X′.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 X′⇝?N′X0.By the above corollary, any reachability problem on a superconservative d-CRN satisfying the conditions of Corollary18 can be easily traced back to that of a subconservative network; hence the problem is equivalent to finding a c∈Z≥0l solution for the respective d-CRN state equation.The reaction network class covered by the above statements might be beneficial in modeling first- and second-order (bio)chemical reaction networks. For a representative example, see Example19 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 Figure4. 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 z=156 and Proposition 15 can be applied.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. aj,Aj and ij,Ij for j=0,…13 denote the active and inactive proteins, respectively. Lower indices denote the number of phosphorylated residues.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.Remark 20.
Let us consider a subconservative (superconservative) d-CRNN=(S,C,R) of n species, m complexes, and l reactions. Assume that N satisfies the conditions of Proposition 17 (Corollary 18). Then for any X0,X′∈Z≥0n initial and target states for which X0⪰M(Γ-),X′⪰M(Γ-) hold we have that the number of distinct trajectories σX satisfying the reachability relation X0⇝?NX′ can be determined in polynomial time in the distance of X0 and X′, given the fixed number of reactions l in the d-CRN. The explanation of this is the following. According to Proposition 17 (Corollary 18) the reachability problem X0⇝?NX′ is equivalent to the existence of a nonnegative integer solution c∈Z≥0l of the state equation X0+Γc=X′. In this way the reachability problem can be reformulated as an ILP feasibility problem in terms of c, 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 reactionsl,(2)
polynomial in the number of constraints, that is, in the number of speciesn,(3)
polynomial in the maximum of the absolute values of the coefficientsΓ, X′-X0.The particular importance of Remark20 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 X0 initial state.


## 6. Computational Example: A Superconservative d-CRN of First-Order Reactions

Let us consider the d-CRN depicted in Figure5. This system is superconservative with a particular conservation vector z=121×1 implying the unboundedness of its reachable state space regardless of the initial state X0. Making use of the above results, the reachability problem of X0⇝?X′ for any X0,X′∈Z≥021 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.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 X0∈Z≥0n.As initial state we considerX0 given by (31) that was randomly generated from [10,100]21. In order to find a target state X′ satisfying the reachability relation X0⇝X′ 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 X0 and X′, we need to solve the following decision problem.(30)Γc=X′-X0c∈Z≥028Clearly, Corollary 18 guarantees that Γc=X′-X0 is satisfied with some c∈Z≥028 if and only if the reachability relation X0⇝X′ holds. Let us consider the following initial and final states.(31)X0=561035876675876060555089587252714871574768⊤(32)X′=504665778895715659544376557840625171536491⊤We found that for the target state X′ 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 fromX0 to X′ with additional constraints. One can observe that X′ results in a significant increase in the number of molecules in the species s5, s6, s20, and s21 and any trajectory from X0 to X′ 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 s5, s6, s20, and s21 can take place through different paths. We assume that the directed paths(33)σS1=s3s12s13s14s15s4,(34)σS2=s3s16s17s4are slow compared to the other ones; hence we wish to minimize the flow through them in order to lower their effect in c. This can be easily expressed by posing addition linear constraints on c as is done in the decision problem (35).(35)Γc=X′-X0c∈Z≥028c15≤10c27≤10We also determined a particular solution c by equipping (35) with the objective function ∑i=128[c]i to be minimized.(36)c=11210611211111848829121816111178131716158171316131916951⊤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.


## 7. 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 high-dimensional superconservative d-CRN.



---

*Source: 1035974-2019-03-26.xml*