Source-only realizations, weakly reversible deficiency one networks and dynamical equivalence

It is in general difficult to analyze the dynamics exhibited by nonlinear systems. They can exhibit a wide range of behaviour ranging from periodicity, chaos, limit cycles etc. One way of analyzing dynamical systems is to represent them as reaction networks. Reaction networks are ubiquitous in biological systems. The dynamics generated by reaction networks often assume polynomial or power-law form. A key feature here is how the graphical structure of the reaction network regulates its dynamics. Such graphical properties include but are not limited to weak reversibility, complex balancing, endotacticity etc. Reaction networks possessing these properties are conjectured to to exhibit dynamical properties like persistence (which implies that no species can go extinct), permanence (which implies convergence to a compact set) and the existence of a globally attracting steady state.

It turns out that there exists networks with different graphical structures that exhibit the same dynamics. This property is called dynamical equivalence and makes the problem of network identification difficult. Dynamical equivalence has also been called “macro-equivalence” by Horn and Jackson [22] and “confoundability” by Craciun and Pantea [10]. A natural application is to infer the dynamics of networks that do not possess properties like weak reversibility or endotacticity, but are dynamically equivalent to weakly reversible or endotactic networks. A large body of literature exists on designing efficient algorithms for detecting weakly reversible or endotactic realizations, given a dynamical system and a set of source vertices. One of our contributions in this paper is to give an upper bound on the number of source vertices required for detecting certain families of reaction networks.

Another quantity that is useful in the analysis of reaction networks is the deficiency of the reaction network. Informally, deficiency is a measure of how far the reaction vectors in the network are from their general position. In particular, weakly reversible deficiency zero reaction networks exhibit robust dynamics. It is known that for weakly reversible deficiency zero reaction networks, there exists a Lyapunov function. In addition, there exists a unique steady state that is locally asymptotically stable [22]. Recently, it was shown [8] that weakly reversible deficiency zero realizations are unique, i.e., there cannot exist two different weakly reversible deficiency zero networks that generate the same dynamical system. Further, an algorithm has been proposed to check whether a weakly reversible deficiency zero realization exists given a dynamical system and a set of source vertices [6]. A natural extension is to ask whether weakly reversible deficiency one realizations are unique. We show that the answer to this question is no. Another question worth exploring is to analyze the property of dynamical equivalence between different types of weakly reversible deficiency one networks, and between weakly reversible deficiency one and weakly reversible deficiency zero networks. We attempt to answer these questions in this paper.

Our paper is organized as follows: In Section 2 we give a short primer on reaction networks. More specifically we introduce reaction networks as graphs embedded in Euclidean space. We define weakly reversible, strongly endotactic, endotactic, consistent and conservative networks. In Section 3 we introduce net vectors and ghost vertices. The idea of network corresponds to the net flux out a vertex. A ghost vertex is a source vertex whose net vector is zero. We will use both these notions later in the paper to prove results about source-only realizations and dynamical equivalence between different families of reaction networks. In Section 4 we define the notion of dynamical equivalence. We show in Lemma 4.3 that two dynamically equivalent weakly reversible or strongly endotactic realizations have the stoichiometric subspace. In Section 5 in Theorem 5.4 we show that certain families of reaction networks possess the following property: A mass-action system $\mathcal{G}_{k}$ is dynamically equivalent to a mass-action system generated by networks from this family iff it is dynamically equivalent to a mass-action system $\tilde{\mathcal{G}}_{\tilde{k}}$ generated by networks from this family that uses only the source vertices of $\mathcal{G}$. The family of such networks include weakly reversible, consistent, endotactic and strongly endotactic networks. Section 6 deals with the uniqueness and non-uniqueness of realizations. In particular, we show that weakly reversible deficiency one realizations are not unique, i.e., there exists different weakly reversible deficiency one networks that generate the same dynamics. In Theorem 6.3 we show that two weakly reversible deficiency one realizations with the property that one realization has all linkage classes of deficiency zero and the other realization has all linkage classes of deficiency zero except one which has deficiency one, cannot be dynamically equivalent.

Reaction networks can be characterized as directed graphs in Euclidean space called Euclidean embedded graphs [4, 5, 11, 3]. In particular, a reaction network is a tuple $\mathcal{G}=(V,E)$, where $V\subset\mathbb{R}^{n}$ is the set of vertices (that correspond to complexes) and $E$ is the set of edges (that correspond to reactions). An edge $(\boldsymbol{y},\boldsymbol{y}^{\prime})\in E$ corresponds to a reaction $\boldsymbol{y}\to\boldsymbol{y}^{\prime}$. Here $\boldsymbol{y}$ is the source vertex and $\boldsymbol{y}^{\prime}$ is the target vertex.

The stoichiometric subspace corresponding to a reaction network $\mathcal{G}$ is given by the set $S:=\rm{span}\{\boldsymbol{y}^{\prime}-\boldsymbol{y}\,|\,\boldsymbol{y}\to\boldsymbol{y}^{\prime}\in E\}$. Given $\boldsymbol{x}_{0}\in\mathbb{R}^{n}_{>0}$, the stoichiometric compatibility class of $\boldsymbol{x}_{0}$ is the polyhedron $(\boldsymbol{x}_{0}+S)\cap\mathbb{R}^{n}_{>0}$.

A linkage class $\mathcal{L}$ is a maximal connected component of $\mathcal{G}$. The deficiency of a reaction is the non-negative integer $\delta$, given by the formula $\delta=|V|-\ell-\rm{dim}(S)$, where $|V|$ is the number of vertices, $\ell$ is the number of linkage classes and $\rm{dim}(S)$ is the dimension of the stoichiometric subspace. Let $V_{\mathcal{S}}$ denote the set of source vertices and let $\rm{conv}(V_{\mathcal{S}})$ denote the convex hull of the source vertices of $\mathcal{G}$. Then

• •

  $\mathcal{G}$ is weakly reversible if every reaction is part of a strongly connected component.

• •

  $\mathcal{G}$ is strongly endotactic [15] if for every vector $\boldsymbol{v}\in\mathbb{R}^{n}$ and $\boldsymbol{y}\to\boldsymbol{y}^{\prime}\in E$ with $\boldsymbol{v}\cdot(\boldsymbol{y}^{\prime}-\boldsymbol{y})<0$, there exists $\tilde{\boldsymbol{y}}\to\tilde{\boldsymbol{y}}^{\prime}\in E$ such that $\boldsymbol{v}\cdot\tilde{\boldsymbol{y}}<\boldsymbol{v}\cdot\boldsymbol{y}$, $\boldsymbol{v}\cdot(\tilde{\boldsymbol{y}}\to\tilde{\boldsymbol{y}}^{\prime})>0$ and $\boldsymbol{v}\cdot\tilde{\boldsymbol{y}}<\boldsymbol{v}\cdot\hat{\boldsymbol{y}}$ for all $\hat{\boldsymbol{y}}\in V_{\mathcal{S}}$.

• •

  $\mathcal{G}$ is endotactic [9] if for every vector $\boldsymbol{v}\in\mathbb{R}^{n}$ and $\boldsymbol{y}\to\boldsymbol{y}^{\prime}\in E$ with $\boldsymbol{v}\cdot(\boldsymbol{y}^{\prime}-\boldsymbol{y})<0$, there exists $\tilde{\boldsymbol{y}}\to\tilde{\boldsymbol{y}}^{\prime}\in E$ such that $\boldsymbol{v}\cdot\tilde{\boldsymbol{y}}<\boldsymbol{v}\cdot\boldsymbol{y}$, $\boldsymbol{v}\cdot(\tilde{\boldsymbol{y}}\to\tilde{\boldsymbol{y}}^{\prime})>0$.

• •

  $\mathcal{G}$ is consistent if there exists positive constants $\lambda_{\boldsymbol{y}\to\boldsymbol{y}^{\prime}}$ such that $\displaystyle\sum_{\boldsymbol{y}\to\boldsymbol{y}^{\prime}\in E}\lambda_{\boldsymbol{y}\to\boldsymbol{y}^{\prime}}(\boldsymbol{y}^{\prime}-\boldsymbol{y})=0$.

• •

  $\mathcal{G}$ is conservative if there exists a vector $\boldsymbol{v}\in\mathbb{R}^{n}_{>0}$ such that $\boldsymbol{v}\in S^{\perp}$.

Figure 2 shows a few examples of such networks.

The chain of inclusions among the above mentioned networks is as follows: Weakly reversible networks consisting of a single linkage class are strongly endotactic. Weakly reversible networks networks are endotactic. Strongly endotactic networks are also endotactic. Further endotactic networks are consistent. It is possible for a network to be

1. (i)

   consistent, but not conservative

2. (ii)

   conservative, but not consistent.

3. (iii)

   consistent and conservative

4. (iv)

   neither consistent nor conservative.

Figure 3 gives examples of networks that are of the four types mentioned above.

Conservative networks are also related to the notion of critical sets, which are relevant in the context of autocatalysis [12, 18, 19, 20]. It can be shown that if a network is conservative, then the set of species is not a critical set [12]. Critical sets also show up in the context of persistence; it is known that if there are no critical siphons in the network, then the dynamics exhibited by it is persistent [2].

It turns out that there is a geometrically intuitive way to deduce whether the network is endotactic or strongly endotactic. This is given by the parallel sweep test [9, 15]. The parallel sweep test for strongly endotactic networks can be summarized as follows [15]

The parallel sweep test for endotactic networks can be summarized as follows [9]

A reaction network can generate a wide array of dynamical systems. If we assume mass-action kinetics [13, 31, 16, 32, 17, 1], the dynamical system generated by $\mathcal{G}$ is given by

$\displaystyle\frac{d\boldsymbol{x}}{dt}=\displaystyle\sum_{\boldsymbol{y}\rightarrow\boldsymbol{y}^{\prime}\in E}k_{\boldsymbol{y}\rightarrow\boldsymbol{y}^{\prime}}{\boldsymbol{x}}^{\boldsymbol{y}}(\boldsymbol{y}^{\prime}-\boldsymbol{y})$

(1)

Let us denote vector of rate constants in this dynamical system by $\boldsymbol{k}$. Then the associated mass-action system will be represented as $(\mathcal{G},\boldsymbol{k})$. A mass-action system is complex balanced if there exists a point $\boldsymbol{x}_{0}\in\mathbb{R}^{n}_{>0}$ such that for any vertex $\boldsymbol{y}\in V$, we have

$\displaystyle\displaystyle\sum_{\boldsymbol{y}\rightarrow\boldsymbol{y}^{\prime}\in E}k_{\boldsymbol{y}\rightarrow\boldsymbol{y}^{\prime}}{\boldsymbol{x}_{0}}^{\boldsymbol{y}}=\displaystyle\sum_{\tilde{\boldsymbol{y}}\rightarrow\boldsymbol{y}\in E}k_{\tilde{\boldsymbol{y}}\rightarrow\boldsymbol{y}}{\boldsymbol{x}_{0}}^{\tilde{\boldsymbol{y}}}$

(2)

As remarked before, graphical properties like weak reversbility, endotacticity and complex balancing are closely related to dynamical properties like persistence, permanence and global stability. Anderson has proved the Global Attractor Conjecture, which says that complex balanced systems have a globally attracting steady state within every stoichiometric compatibility class when the associated reaction graph consists of a single linkage class. It is know from the work of Craciun, Nazarov and Pantea [9] that two-dimensional endotactic networks are permanent. This has been extended by Pantea [24] to endotactic networks with two-dimensional subspace. Gopalkrishnan, Miller and Shiu [15] have shown that strongly endotactic networks are permament. Recently, Craciun [4] has proposed a proof of the Global Attractor Conjecture.

The idea of this section is to introduce the notions of net vectors and ghost vertices, which will be used later in the paper to analyze the dynamics of several families of networks. The dynamical system generated by $(\mathcal{G},\boldsymbol{k})$ (given by Equation 1) can be rewritten as

$\displaystyle\frac{d\boldsymbol{x}}{dt}=\displaystyle\sum_{\boldsymbol{y}\in V_{S}}{\boldsymbol{x}}^{\boldsymbol{y}}\displaystyle\sum_{\boldsymbol{y}^{\prime}\in V}k_{\boldsymbol{y}\rightarrow\boldsymbol{y}^{\prime}}{\boldsymbol{x}}^{\boldsymbol{y}}(\boldsymbol{y}^{\prime}-\boldsymbol{y})$

(3)

Then $\displaystyle\sum_{\boldsymbol{y}^{\prime}\in V}k_{\boldsymbol{y}\rightarrow\boldsymbol{y}^{\prime}}{\boldsymbol{x}}^{\boldsymbol{y}}(\boldsymbol{y}^{\prime}-\boldsymbol{y})$ will denote the net vector corresponding to the source vertex $\boldsymbol{y}$.

If the net reaction vector corresponding to a source vertex $\boldsymbol{y}_{0}$ is zero, then $\boldsymbol{y}_{0}$ will be called a ghost vertex. Previous works [7] have referred to such vertices as virtual sources. Figure 4 shows an example of a ghost vertex.

Associated with every reaction network is a dynamical system. However, it is possible for two different reaction networks to generate the same dynamical system. This phenomenon is called dynamical equivalence. It is also referred to as macro-equivalence [22] or confoundability [10]. We formally define the notion of dynamical equivalence below.

Figure 5 shows a a few examples of dynamically equivalent networks.

Dynamical systems generated by reaction networks are also sometimes referred to as realizations. The next lemma collects a useful observation about the stoichiometric subspaces of dynamically equivalent weakly reversible reaction networks.

Network identification is a challenging problem in general due to the phenomenon of dynamical equivalence. However, certain properties of the network like complex balancing, weak reversibility are often required to design efficient systems [23, 30]. There exists a large body of work that design algorithms for detecting properties like weak reversibility, complex balancing or for finding dense realization, i.e., realizations that have the maximum number of reactions. Many of these approaches use techniques like linear and mixed linear integer programming to perform a constrained optimization [30, 26, 27, 25]. In such approaches, the set of source vertices is part of the input search space over which the optimization is performed. Previous works [28, 29, 7] have shown that to deduce whether a detailed balanced, complex balanced or weakly reversible realizations exists, it suffices to check over the space of input source vertices. In our paper, we show that this property extends to endotactic, strongly endotactic, consistent and conservative networks. Theorem 5.4 illustrates this point.

Before we move to Theorem 5.4, we recall the following lemma from linear algebra which will be used in the context of consistent networks.

As remarked before, realizations of reaction networks may not be unique due to the property of dynamical equivalence. However, under certain restrictions on the reaction network structure, one can show that the realization is unique. For example in [8], Craciun et. al. showed that weakly reversible deficiency one realizations are unique. Weakly reversible deficiency zero networks form an important class of reaction networks given the rich dynamics that they exhibit. Weakly reversible deficiency zero reaction networks are known to be complex balanced for all values of the rate constants [21]. The property of being complex balanced implies that there exists a Lyapunov function that decreases along trajectories; further there exists a unique equilibria in each stoichiometric compatibility class which is locally asymptotically stable [22]. In the language of Feinberg [14], the dynamical system is quasi-thermodynamic.

Note that it follows from the deficiency zero theorem that for a deficiency zero network that is not weakly reversible, there exists no positive equilibrium for all choices of $\boldsymbol{k}$. In other words, a consistent deficiency zero reaction network is weakly reversible. Since weakly reversible deficiency zero realizations are unique [8], this implies that consistent deficiency zero realizations are unique.

A natural extension is to ask the following question: “Are weakly reversible deficiency one realizations unique?” We show that the answer to this question is no in general, i.e., weakly reversible deficiency one realizations are not unique.

A network with deficiency one can be classified into the following three types:

1. (i)

   All linkages classes of deficiency zero.

2. (ii)

   One linkage of deficiency one.

3. (iii)

   One linkage of deficiency one and the remaining linkage classes of deficiency zero.

It is interesting to see the interplay of dynamical equivalence between the above mentioned cases in the context of weakly reversible deficiency one networks; and also the relationship between the dynamics of weakly reversible deficiency zero and deficiency one networks.

Figure 7 shows an example of two dynamically equivalent weakly reversible networks that have deficiency one and each of them consists of two linkage classes of deficiency zero.

Figure 8 shows an example of two dynamically equivalent weakly reversible networks that have deficiency one and each of them consists of a single linkage class.

One can also generate examples of reaction networks consisting of different number of linkage classes. Figure 9.(a) shows a weakly reversible network of deficiency zero having two linkage classes. This is dynamically equivalent to the network in Figure 9.(b) which is a weakly reversible network of deficiency one having one linkage class.

In Theorem 6.3, we show that two reaction networks having a deficiency of one (one reaction network whose linkage classes have deficiency zero and another reaction network whose linkage classes all of deficiency zero except one which has deficiency one) and same number of linkages cannot be dynamically equivalent. Towards this, we need a result from [8] as the next remark states.

Network identification is an important problem in the analysis of biological systems. The property of dynamical equivalence makes the problem of network identification quite difficult. Identifying a particular type of realization like weakly reversible or endotactic given a dynamical system and a set of source vertices places no premature bound on the number of vertices in the network. In this paper, we show that to identify a strongly endotactic, endotactic, consistent or conservative realizations, it suffices to consider only the given source vertices. This builds upon earlier work by Craciun, Jin and Yu [7], who showed that this property is true for detailed balanced, complex balanced and weakly reversible realizations.

In addition, we have studied the uniqueness/non-uniqueness of realizations for weakly reversible networks having low deficiency. In particular, we show that weakly reversible deficiency one realizations are not unique. This builds upon earlier work by Craciun, Jin and Yu [8] who showed that weakly reversible deficiency zero realizations are unique. Further, we have tried to characterize the dynamical relationships between several types of weakly reversible deficiency one networks. In particular, Theorem 6.3 says that two weakly reversible deficiency one realizations such that one realization has all linkage classes of deficiency zero and the other realization has all linkage classes of deficiency zero except one which has deficiency one, cannot be dynamically equivalent.

A future line of research is to extend the set of inclusions obtained among several types of deficiency one networks in the context of dynamical equivalence to networks of higher deficiency. Another direction worthy of exploration is to design efficient algorithms to detect weakly reversible networks of high deficiency.

I thank Gheorghe Craciun for useful discussions and suggestions for improving the presentation of the paper.