Complex Balancing Reconstructed to the Asymptotic Stability of Mass-action Chemical Reaction Networks with Conservation Laws

TL;DR

This paper introduces a generalized Lyapunov function for analyzing the asymptotic stability of mass-action systems with conservation laws by reconstructing networks and examining their stability properties.

Contribution

It develops a novel reconstruction method and a generalized pseudo-Helmholtz function to analyze stability in complex balanced and more general mass-action systems.

Findings

Reverse reconstruction has a unique equilibrium in each positive class.

Complex balanced reconstruction ensures asymptotic stability of the original network.

The method provides a systematic way to compute reconstructions with conservation laws.

Abstract

Motivated by the fact that the pseudo-Helmholtz function is a valid Lyapunov function for characterizing asymptotic stability of complex balanced mass action systems (MASs), this paper develops the generalized pseudo-Helmholtz function for stability analysis for more general MASs assisted with conservation laws. The key technique is to transform the original network into two different MASs, defined by reconstruction and reverse reconstruction, with an important aspect that the dynamics of the original network for free species is equivalent to that of the reverse reconstruction. Stability analysis of the original network is then conducted based on an analysis of how stability properties are retained from the original network to the reverse reconstruction. We prove that the reverse reconstruction possesses only an equilibrium in each positive stoichiometric compatibility class if the corresponding reconstruction is complex balanced. Under this complex balanced reconstruction strategy, the asymptotic stability of the reverse reconstruction, which also applies to the original network, is thus reached by taking the generalized pseudo-Helmholtz function as the Lyapunov function. To facilitate applications, we further provide a systematic method for computing complex balanced reconstructions assisted with conservation laws. Some representative examples are presented to exhibit the validity of the complex balanced reconstruction strategy.

Tables

Table 1. Table 1: Some mass action systems and their complex balanced reconstructions

No.	MAS $(𝒳,𝒞,ℛ,𝒦)$	Equilibrium $x^{\ast }$	Complex Balanced Reconstruction $(𝒳,{\widehat{𝒞}}_C,{\widehat{ℛ}}_C,{\widehat{𝒦}}_C)$	Reconstructing Matrix $D$
1		$\left[\begin{array}{c}1\\ 2\end{array}\right]$
2		$\left[\begin{array}{c}1\\ 1\end{array}\right]$
3		$\left[\begin{array}{c}\frac{1}{3}\\ \frac{1}{3}\\ \frac{1}{3}\end{array}\right]$
4		$\left[\begin{array}{c}0.5\\ 0.5\\ 0.5\end{array}\right]$
5		$\left[\begin{array}{c}0.3077\\ 0.8077\\ 1.0128\\ 0.6\end{array}\right]$
6		$\left[\begin{array}{c}1\\ 1\\ 2\\ 2\end{array}\right]$
§ $\varnothing$ represents zero complex.

Equations

$$3X_2\to 3X_1\to 2X_1+X_2$$

$$S=ZB.$$

$$v_j(x)=k_j{x}^{Z_{\cdot j}}:=k_j\prod _{i=1}^n{x}_i^{Z_{ij}},$$

$$\dot{x}=Sv(x),x\in {\mathbb{R}}_{\ge 0}^n,$$

$$x=x_0+\sum _{j=1}^r{S}_{\cdot j}{\int }_0^t{v}_j(\tau )d\tau .$$

$$G(x)=\sum _{i=1}^n(x_i^{\ast }-x_i-x_i\text{ln}\frac{x_i^{\ast }}{x_i}),x\in {\mathbb{R}}_{>0}^n$$

$$G_e(x)=\sum _{i=1}^n{d}_i(x_i^{\ast }-x_i-x_i\text{ln}\frac{x_i^{\ast }}{x_i}),x\in {\mathbb{R}}_{>0}^n,d_i\in {\mathbb{R}}_{\text{textgreater}0}.$$

$$\tilde{G}(x)=\sum _{i={\tau }_1}^{{\tau }_p}{d}_i(x_i^{\ast }-x_i-x_i\ln \frac{x_i^{\ast }}{x_i})$$

$$\tilde{G}(x)=\sum _{i=1}^p{d}_i(x_i^{\ast }-x_i-x_i\ln \frac{x_i^{\ast }}{x_i}).$$

$$x^{\ast }=\left[\begin{array}{c}{x}_{\perp }^{\ast }\\ x_{\top }^{\ast }\end{array}\right]$$

$$C^{\top }x=\left[\begin{array}{cc}{C}_l^{\top }& C_r^{\top }\end{array}\right]\left[\begin{array}{c}{x}_{\perp }\\ x_{\top }\end{array}\right]=C_l^{\top }{x}_{\perp }+C_r^{\top }{x}_{\top }=C^{\top }{x}^{\ast },$$

$$x_{\top }=C_r^{-\top }{C}_l^{\top }(x_{\perp }^{\ast }-x_{\perp })+x_{\top }^{\ast },$$

$$x_{\top }={C^{\prime }}_r^{-\top }{C^{\prime }}_l^{\top }(x_{\perp }^{\ast }-x_{\perp })+x_{\top }^{\ast }.$$

$$\tilde{G}(x)=\sum _{i=1}^{n-q}{d}_i(x_i^{\ast }-x_i-x_i\ln \frac{x_i^{\ast }}{x_i}),\forall x\in {\mathbb{R}}_{\text{textgreater}0}^{n}and\forall d_i\text{textgreater}0,$$

$$\nabla \tilde{G}(x)=\mathrm{diag}(d_1,\cdots ,d_{n-q},0,\cdots ,0)\nabla G(x).$$

$$\dot{\tilde{G}}(x)$$

$$={\nabla }^{\top }G(x)\mathrm{diag}(d_1,\cdots ,d_{n-q},0,\cdots ,0)Sv(x)$$

$$\le 0.$$

$$DSv(x)=\left[\begin{array}{c}\hat{S}\hat{v}(\hat{x})\\ \text{mathbbold}{0}_q\end{array}\right],$$

$$D=\left[\begin{array}{cc}{D}_1& \text{mathbbold}{0}_{(n-q)\times q}\\ C_l^{\top }& C_r^{\top }\end{array}\right],\hat{x}=x_{\perp }\in {\mathbb{R}}_{\ge 0}^{n-q},\hat{S}\in {\mathbb{Z}}^{(n-q)\times \hat{r}}\mathrm{and}\hat{v}:{\mathbb{R}}_{\ge 0}^{n-q}\to {\mathbb{R}}_{\ge 0}^{\hat{r}},$$

$$D_1{S}_{1}v(x)=\hat{S}\hat{v}(\hat{x}),$$

$$D^{-1}=\left[\begin{array}{cc}{D}_1^{-1}& \text{mathbbold}{0}_{(n-q)\times q}\\ -C_r^{-\top }{C}_l^{\top }{D}_1^{-1}& C_r^{-\top }\end{array}\right].$$

$$\frac{d\hat{V}(\hat{x})}{dt}$$

$$V(x)={\int }_{x^{\ast }}^x{\nabla }^{\top }\hat{V}(\hat{y})D_{2}dy,$$

$$V(x)={\nabla }^{\top }\hat{V}(\hat{\bar{x}})D_2(x-x^{\ast })=\sum _{i=1}^{n-q}{d}_i\frac{\partial \hat{V}(\hat{\bar{x}})}{\partial {\hat{x}}_i}(x_i-x_i^{\ast }),$$

$$H(V(x))=\left[\begin{array}{cc}\mathrm{diag}\left(d_i\frac{{\partial }^2\hat{V}(\hat{x})}{\partial {\hat{x}}_i^2}\right)& \text{mathbbold}{0}_{(n-q)\times q}\\ \text{mathbbold}{0}_{q\times (n-q)}& \text{mathbbold}{0}_{q\times q}\end{array}\right],i=1,\cdots ,n-q.$$

$$\frac{dV(x)}{dt}$$

$$={\nabla }^{\top }\hat{V}(\hat{x})D_{2}Sv(x)$$

$$={\nabla }^{\top }\hat{V}(\hat{x})\left[\begin{array}{cc}{D}_1& \text{mathbbold}{0}_{(n-q)\times q}\end{array}\right]\left[\begin{array}{c}{S}_1\\ S_2\end{array}\right]v(x)$$

$$={\nabla }^{\top }\hat{V}(\hat{x})D_1{S}_{1}v(x)$$

$$\le 0.$$

$$\dot{\tilde{x}}=\tilde{S}\tilde{v}(\tilde{x})$$

$$Sv(x)=\left[\begin{array}{c}{E}_{n-q}\\ -C_r^{-\top }{C}_l^{\top }\end{array}\right]\tilde{S}\tilde{v}(\tilde{x}),$$

Full text

{AMS}

34D20, 80A30, 93C15, 93D20

1 Introduction

Chemical reaction networks (CRNs) taken with mass action kinetics arise often in chemistry, biology and engineering. The dynamics of a mass action system (MAS) can be modeled by polynomial ordinary differential equations (ODEs) with respect to the concentration of each chemical species that is a nonnegative real number. Although the concentration evolution takes place in a certain invariant set termed the stoichiometric compatibility class, these ODEs will typically lead to a high-dimensional and complex nonlinear system, often referred to as a kinetic system, with unknown model parameters. The model parameters and the network structure both have great effect on the kinetic system, so it is usually extremely difficult to characterize the dynamical properties. However, the past decades witnessed the rapid progress in related directions, including characterizing stability [1, 3, 8, 9, 13, 21, 23, 27], multi-stationarity [6, 10, 17] and persistence [2, 4, 5, 7, 20]. Now, the dynamical properties have gradually been coming to light, especially for those network systems with special structures, such as reversible and weakly reversible structures. In line with these studies, this paper continues to be concerned with characterizing (local asymptotic) stability of equilibria in general MASs.

The pioneering work on stability analysis of equilibria in MASs originated from [13], in which complex balanced MASs were proved to have no self-sustaining oscillation or bistability. This paper [13] reported the famous Deficiency Zero Theorem, saying that for a deficiency zero weakly reversible CRN with any choice of rate constants, the dynamical equation admits a unique positive equilibrium in each positive stoichiometric compatibility class, and these positive equilibria are locally asymptotically stable. This theorem was further extended to the Deficiency One Theorem [10], which gives alternative conditions for the existence and uniqueness of equilibrium in each positive stoichiometric compatibility class, i.e., weakly reversible but not necessary deficiency zero. Recently, the stability issue of complex balanced MASs [21], and of detailed balanced MASs (a special case of complex balanced systems) [27] were revisited, and a compact mathematical formulation about the network were derived that allows to characterize the set of positive equilibria and their stability properties, easily. For complex balanced systems, if they are assumed to admit persistence [9], then the stability of equilibria is global [2, 20].

In this paper, we consider a class of CRNs, like

[TABLE]

which are not weakly reversible, for characterising stability of the corresponding MASs. This class of CRNs are often encountered in enzymic catalytic reactions and epidemic models [12]. With respect to stability analysis, the key issue is how to construct the Lyapunov function. Anderson et. al. [3] treated scaling limits of non-equilibrium potential as Lyapunov functions for some special species balanced systems, for example, birth-death systems, etc. Al-Radhawi and Angeli [1] developed Lyapunov functions that are piecewise linear in rates. The existence of such Lyapunov functions could guarantee stability of equilibria, and further serve to establish asymptotic stability if the Lyapunov functions also satisfy the LaSalle’s condition. Clearly, it poses a great challenge to find a Lyapunov function for this class of MASs. Another alternative strategy for capturing their stability properties is to apply the (asymptotic) stability results of complex balanced systems, after all, for which there is a ready-made Lyapunov function, the pseudo-Helmholtz function. Linear conjugacy [14, 15, 16] is a related notion for this purpose that defines a non-degenerative linear transformation to map the mass action system under study to another one. If the transformed system is complex balanced, then the original system has local asymptotic stability [14]. To obtain expected linear conjugacy, this strategy usually calls for solving linear programming and mixed-integer linear programming problems, analogous to computing complex balanced or detailed balanced realizations [25].

Motivated by the facts that the pseudo-Helmholtz function is a valid Lyapunov function for asymptotic stability analysis for complex balanced MASs and general MASs admitting a complex balanced MAS as the linear conjugacy, we develop the generalized pseudo-Helmholtz function aiming to characterize asymptotic stability for general MASs. Based on the mass conservation law, the reconstruction and reverse reconstruction notions are adaptively developed, which are highly related to the original network in properties of equilibria and stability. More importantly, the dynamics of the reverse reconstruction is equivalent to that of the original network at the free species. The stability analysis for the latter is thus transformed into that for the former. Taking the generalized pseudo-Helmholtz function as the Lyapunov function, the reverse reconstruction is proved asymptotically stable when the corresponding reconstruction is complex balanced. The key point lies that under this condition the reverse reconstruction possesses only an equilibrium in each positive stoichiometric compatibility class. Utilizing dynamical equivalence and mass conservation law, the state of the original network is proved to converges asymptotically to the equilibrium of interest. To facilitate applications of the proposed strategy of stability analysis, we also present a systematic method to compute complex balanced reconstructions based on the computation algorithms for realizations [24, 25, 26].

The remainder of this paper is organized as follows. Section 2 gives a brief introduction on CRNs and further formulates the research motivation. In Section 3, the generalized pseudo-Helmholtz function is developed with the potential to behave as a Lyapunov function. Section 4 defines the reconstruction concept to match the use of the generalized pseudo-Helmholtz function for characterizing stability of general MASs. This is followed by further defining reverse reconstruction and complex balanced reconstruction to characterize asymptotic stability of MASs in Section 5. Section 6 provides a systematic method to compute complex balanced reconstructions, and then illustrates it through some representative examples. Finally, Section 7 concludes the paper.

Notations: Throughout the paper, $\mathbb{R}^{n}$, $\mathbb{R}^{n}_{>0}$, $\mathbb{R}^{n}_{\geq 0}$, $\mathbb{Z}^{n}$, $\mathbb{Z}^{n}_{\geq 0}$ denote the space of $n$ dimensional real, positive real, nonnegative real, integer and nonnegative integer vectors, respectively; for any $x\in\mathbb{R}^{n}_{>0}$, the mapping $\mbox{Ln}:\mathbb{R}^{n}_{\textgreater 0}\rightarrow\mathbb{R}^{n},x\mapsto\mbox{Ln}(x)$ is defined by $(\mbox{Ln}(x))_{i}:=\ln(x_{i})$, and the mapping $\mbox{Exp}(x)$ is similarly defined by $(\mbox{Exp}(x))_{i}=\mbox{exp}(x_{i})$; for matrix $A$, $A_{\cdot i}$ refers to the $i$th column of $A$; $\mathbbold{0}_{n}$ and $\mathbbold{1}_{n}$ denote the vector with all the entries equal to [math] and $1$, respectively; matrix $D=\text{diag}(d_{i})\in\mathbb{R}^{n\times n}$ is a diagonal matrix with the $i$th diagonal element to be $d_{i}$; $E_{n}$ denotes the $n\times n$ identity matrix; $\mathrm{Card}(\cdot)$ represents the number of elements in set.

2 Chemical reaction networks and Motivation statement

In this section, we will review some basic concepts related to CRNs [10, 13] and set forth the motivation of the current study.

2.1 Chemical reaction networks

For a CRN involving $n$ species $X_{i}~{}(i=1,\cdots,n)$, $c$ complexes $\mathcal{C}_{l}~{}(l=1,\cdots,c)$ and $r$ reactions $\mathcal{R}_{j}~{}(j=1,\cdots,r)$, it is defined as follows [10].

Definition 2.1** (Chemical Reaction Network).**

A CRN is composed of three finite sets:

• (1)

A set $\mathcal{X}=\bigcup_{i=1}^{n}\{X_{i}\}$ of chemical species;

• (2)

*A set $\mathcal{C}=\bigcup_{j=1}^{r}\{Z_{\cdot j},Z^{\prime}_{\cdot j}\}$ of complexes with $\mathrm{Card}(\mathcal{C})=c$, $Z_{\cdot j}$, $Z^{\prime}_{\cdot j}\in\mathbb{Z}^{n}_{\geq 0}$ and the *$i$th entry of $Z_{\cdot j}$ representing the coefficient of $X_{i}$ in this complex;

• (3)

A set $\mathcal{R}=\bigcup_{j=1}^{r}\{Z_{\cdot j}\to Z^{\prime}_{\cdot j}\}$ of reactions such that $\forall~{}Z_{\cdot j}\in\mathcal{C}$, $Z_{\cdot j}\to Z_{\cdot j}\notin\mathcal{R}$ but $\exists~{}Z^{\prime}_{\cdot j}\in\mathcal{C}$ to support either $Z_{\cdot j}\to Z^{\prime}_{\cdot j}\in\mathcal{R}$ or $Z^{\prime}_{\cdot j}\to Z_{\cdot j}\in\mathcal{R}$.

*The triple $(\mathcal{X,C,R})$ is usually used to express a CRN. *

The complexes $Z_{\cdot j},Z^{\prime}_{\cdot j}$ actually express the existence of all of $n$ species in the reactant complex and the product complex, respectively. If $Z_{ij}\neq 0$ then there includes the species $X_{i}$ in the complex $Z_{.j}$, and vice versa. Organize distinct complexes according to $Z=[Z_{\cdot 1},\cdots,Z_{\cdot c}]\in\mathbb{Z}^{n\times c}_{\geq 0}$ to form the complex stoichiometric matrix, whose $l$th column expresses the $l$th complex $Z_{\cdot l}$ in the CRN.

For any CRN, its structure can be described by a digraph, termed as reaction graph [13], with the vertices indicating the complexes and the directed edges representing the reactions. Every tail vertex in a directed edge denotes a reactant complex while the head vertex denotes a product complex. Mathematically, this digraph can be defined by an incident matrix $B\in\mathbb{Z}^{c\times r}$ with $B_{lj}=1$ if vertex $l$ is the head vertex of edge $j$ and $B_{lj}=-1$ if vertex $l$ is the tail vertex of edge $j$, while $B_{lj}=0$ otherwise. The incident matrix and the complex stoichiometric matrix further define the stoichiometric matrix $S\in\mathbb{Z}^{n\times r}$ satisfying

[TABLE]

It is not difficult to find that the $j$th column of $S$ is given by $S_{\cdot j}=Z^{\prime}_{\cdot j}-Z_{\cdot j}$, $j=1,\cdots,r$, which is also called the reaction vector of the $j$th reaction. The rank of $S$ usually represents the rank the CRN $(\mathcal{X,C,R})$, denoted by $s$, satisfying $s\leq\min\{n,r\}$.

Definition 2.2** (Weakly Reversible CRN).**

*A CRN $(\mathcal{X,C,R})$ is weakly reversible if for each reaction $Z_{\cdot j}\to Z^{\prime}_{\cdot j}\in\mathcal{R}$, there exists a sequence of reactions in $\mathcal{R}$, which starts with $Z^{\prime}_{\cdot j}$ and ends with $Z_{\cdot j}$, i.e., $Z^{\prime}_{\cdot j}\to Z_{\cdot j_{1}}\in\mathcal{R}$, $Z_{\cdot j_{1}}\to Z_{\cdot j_{2}}\in\mathcal{R}$, $\cdots$, $Z_{\cdot j_{(m-1)}}\to Z_{\cdot j_{m}}\in\mathcal{R}$, $Z_{\cdot j_{m}}\to Z_{\cdot j}\in\mathcal{R}$. *

Definition 2.3** (Stoichiometric Compatibility Class).**

*For a CRN $(\mathcal{X,C,R})$, the linear subspace $\mathrm{Im}(S)=\mathrm{span}\{Z^{\prime}_{\cdot 1}-Z_{\cdot 1},\cdots,Z^{\prime}_{\cdot r}-Z_{\cdot r}\}$ is called the stoichiometric subspace. Let $x_{0}\in\mathbb{R}^{n}_{\geq 0}$, then the sets $\mathscr{S}(x_{0})=\{x_{0}+\xi|\xi\in\mathrm{Im}(S)\}$, $\mathscr{S}(x_{0})\cap\mathbb{R}^{n}_{\geq 0}$ and $\mathscr{S}(x_{0})\cap\mathbb{R}^{n}_{>0}$ are named the stoichiometric compatibility class, the nonnegative stoichiometric compatibility class and the positive stoichiometric compatibility class of $x_{0}$, respectively. *

When a CRN is assigned a mass action kinetics, the rate $v_{j}$ for reaction $Z_{\cdot j}\to Z^{\prime}_{\cdot j}$ is evaluated by

[TABLE]

where $k_{j}\in\mathbb{R}_{\textgreater 0}$ is the rate constant for this reaction, and $x\in\mathbb{R}^{n}_{\geq 0}$ is the vector of concentrations $x_{i}$ of the chemical species $X_{i}$.

Definition 2.4** (Mass Action System).**

*A mass action system is a CRN $(\mathcal{X,C,R})$ taken together with a mass action kinetics. Denote the vector of reaction rate constants by $\mathcal{K}=(k_{1},\cdots,k_{r})$ with $k_{j}$ representing the rate constant for reaction $Z_{\cdot j}\to Z^{\prime}_{\cdot j},~{}j=1,\cdots,r$, then a MAS is often referred to as a quadruple $(\mathcal{X},\mathcal{C},\mathcal{R},\mathcal{K})$. *

The dynamics of a MAS $(\mathcal{X},\mathcal{C},\mathcal{R},\mathcal{K})$ that captures the changes of concentrations of every species is thus expressed as

[TABLE]

where $v(x)\in\mathbb{R}_{\geq 0}^{r}$ is the vector of reaction rate $v_{j}(x)$ of $\mathcal{R}_{j}$, $j=1,\cdots,r$.

Remark 2.5**.**

Whatever the initial state $x_{0}$ is, the state of a MAS $(\mathcal{X},\mathcal{C},\mathcal{R},\mathcal{K})$ will evolve in the nonnegative stoichiometric compatibility class of $x_{0}$, i.e., in $\mathscr{S}(x_{0})\cap\mathbb{R}^{n}_{\geq 0}$, which can be easily verified by integrating Eq. 3 from [math] to $t$, that is,

[TABLE]

Remark 2.6**.**

*The dynamical equation Eq. 3 is essentially a polynomial system of ODEs, often referred to as a kinetic system. Note that there is not a one-to-one relation between a kinetic system and a MAS. Theoretically, the former can be realized by structurally and/or parametrically different MASs that are dynamically equivalent [13]. *

Here, the dynamical equivalence is defined as follows.

Definition 2.7** (Dynamical Equivalence).**

*Two MASs are said to be dynamically equivalent if they have identical differential equations produced under the assumption of mass action kinetics. *

Definition 2.8** (Equilibrium).**

*A concentration vector $x^{\ast}\in\mathbb{R}^{n}_{>0}$ is an equilibrium of a MAS $(\mathcal{X},\mathcal{C},\mathcal{R},\mathcal{K})$ if the dynamics of Eq. 3 supports $Sv(x^{\ast})=\mathbbold{0}_{n}$, and a complex balanced equilibrium if $Bv(x^{\ast})=\mathbbold{0}_{c}$. A MAS $(\mathcal{X},\mathcal{C},\mathcal{R},\mathcal{K})$ is called complex balanced if it admits a complex balanced equilibrium. *

Remark 2.9**.**

*Any complex balanced equilibrium is an equilibrium, but not necessarily the other way around, since $Z$ need not be injective. *

A nice property of complex balanced MASs is that each positive stoichiometric compatibility class contains only one equilibrium, and each such equilibrium is locally asymptotically stable[13].

Theorem 2.10** (Locally asymptotic stability of complex balanced equilibrium [13]).**

*Let $x^{*}\in\mathbb{R}^{n}_{>0}$ be any complex balanced equilibrium of the MAS $(\mathcal{X},\mathcal{C},\mathcal{R},\mathcal{K})$, then $x^{*}$ is locally asymptotically stable with respect to any initial condition in $\mathscr{S}(x^{\ast})\cap\mathbb{R}^{n}_{>0}$ near $x^{*}$. *

The proof of the above theorem is carried through with the following well-known Lyapunov function

[TABLE]

which is referred to as the pseudo-Helmholtz function [13] in the context.

Remark 2.11**.**

*Note that the domain of the pseudo-Helmholtz function $G(x)$ is $\mathbb{R}^{n}_{>0}$, which seems not completely consistent with the nonnegative stoichiometric compatibility class of the equilibrium, i.e., $\mathscr{S}(x^{\ast})\cap\mathbb{R}^{n}_{\geq 0}$, the domain of the complex balanced MAS. Actually, the domain of the complex balanced MAS can shrink to the positive stoichiometric compatibility class of the equilibrium, $\mathscr{S}(x^{\ast})\cap\mathbb{R}^{n}_{\textgreater 0}$, as long as the initial point is selected carefully[10]. The reason is that $\lim_{x\to\mathbbold{0}_{n}}G(x)=\sum_{i=1}^{n}x^{*}_{i}$ and $\frac{dG(x)}{dt}\leq 0$, as long as the energy value at an initial point, evaluated by $G(x)$, is lower than $\lim_{x\to\mathbbold{0}_{n}}G(x)=\sum_{i=1}^{n}x^{*}_{i}$, then the trajectory will not run towards zero point. A feasible initial points set is $\{x_{0}|G(x_{0})\textless\lim_{x\to\mathbbold{0}_{n}}G(x)\}\cap\mathscr{S}(x^{\ast})\cap\mathbb{R}^{n}_{\textgreater 0}$. *

2.2 Motivation statement

The pseudo-Helmholtz function $G(x)$ exhibits high efficacy in characterizing the asymptotic stability of a complex balanced equilibrium, but it is not in general Lyapunov function for systems that are not complex balanced. In order to reach the latter purpose, it is necessary to develop new Lyapunov functions. Due to high relevance between complex balanced equilibria and general ones (the former is a subset of the latter), a possible solution to find the new Lyapunov function is to extend the pseudo-Helmholtz function to a more general form. A naive extension is to multiply every term in the function $G(x)$ by a positive number $d_{i},~{}i=1,\cdots,n$, then we get a new function

[TABLE]

Clearly, like $G(x)$, this extension also has a strict minimum 0 at $x^{*}$, which indicates that $G_{e}(x)$ satisfies the necessary condition to behave as a Lyapunov function, and is thus potential for stability analysis for some MASs. Following this superficial motivation, we will discuss a more general extension of the pseudo-Helmholtz function so that more MASs may be served, and analyze its properties from viewpoint of becoming the Lyapunov function in the following.

3 Generalized pseudo-Helmholtz function

A more general extension of the pseudo-Helmholtz function is to further set the number of terms in $G(x)$ not fixed. This leads to the definition of the generalized pseudo-Helmholtz function.

Definition 3.1** (Generalized pseudo-Helmholtz function).**

Let $x^{\ast}\in\mathbb{R}^{n}_{>0}$ be an equilibrium of a MAS $(\mathcal{X,C,R,K})$ governed by Eq. 3. The function $\tilde{G}:\mathbb{R}^{n}_{\textgreater 0}\mapsto\mathbb{R}_{\geq 0}$ given by

[TABLE]

*is called the generalized pseudo-Helmholtz function, where $d_{i}\textgreater 0$, $\tau_{1}\leq\cdots\leq\tau_{p}$ and $\{\tau_{1},\cdots,\tau_{p}\}\subseteq\mathcal{I}=\{1,\cdots,n\}$. *

Remark 3.2**.**

*For any $x\in\mathbb{R}^{n}_{\textgreater 0}$, the generalized pseudo-Helmholtz function satisfies $\tilde{G}(x)\geq 0$ with equality hold if and only if $(x_{\tau_{1}},\cdots,x_{\tau_{p}})^{\top}=(x_{\tau_{1}}^{*},\cdots,x_{\tau_{p}}^{*})^{\top}$. *

For simplicity but without loss of generality, we set $\{\tau_{1},\cdots,\tau_{p}\}=\{1,\cdots,p\}$ in (8) for the subsequent analysis, i.e., keeping the first $p$ terms of $G_{e}(x)$. As a result,

[TABLE]

Denote $(x_{1},\cdots,x_{p})^{\top}=x_{\bot}$ and $(x_{p+1},\cdots,x_{n})^{\top}=x_{\top}$, then $\forall~{}x_{\top}\in\mathbb{R}^{n-p}_{\textgreater 0}$ the function $\tilde{G}(x)$ attains the minimum [math] at $x=({x_{\bot}^{\ast}}^{\top},x_{\top}^{\top})^{\top}$. This implies that the equilibrium

[TABLE]

is not the strict minimum point for the generalized pseudo-Helmholtz function. At this point, $\tilde{G}(x)$ seems not able to behave as a Lyapunov function. However, if the degree of freedom of the concerning MAS is $p$, and the remaining $n-p$ state variables may be determined uniquely by those free $p$ variables, then $\tilde{G}(x)$ may be a Lyapunov function. We analyze the number of non-free variables in a MAS $(\mathcal{X,C,R,K})$ through considering the law of mass conservation.

Definition 3.3** (Mass conservation law [13]).**

*The law of conservation of mass states that mass in an isolated system is neither created nor destroyed by chemical reactions or physical transformations. For a MAS $(\mathcal{X,C,R,K})$ governed by Eq. 3, this means that there exists a $n$-dimensional positive vector $\rho\in\mathbb{R}^{n}_{>0}$ such that $\rho^{\top}S=\mathbbold{0}_{r}^{\top}$, and $\rho$ is called a conserved vector. *

Remark 3.4**.**

*Clearly, $\rho\in\mathrm{Ker}(S^{\top})$. Moreover, for a mass-conserved MAS, there must exist a conserved vector, the most simple candidate for which is the vector $M$ of molecular weight $M_{i}$ of the species $X_{i}$, $i=1,\cdots,n$. *

We then define an important concept, conserved matrix, that will play a key role on the subsequent analysis.

Definition 3.5** (Conserved Matrix).**

*Consider a mass-conserved MAS $(\mathcal{X,C,R,K})$ governed by Eq. 3 with $\mathrm{rank}(S)=s$. Assume that there exist $q$ linearly independent vectors $\xi_{m}\in\mathbb{R}^{n}_{>0}~{}(m=1,\cdots,q)$ in $\mathrm{Ker}(S^{\top})$, where $1\leq q\leq\mathrm{dim}\left(\mathrm{Ker}(S^{\top})\right)=n-s$, then the matrix $C\in\mathbb{R}^{n\times q}_{\textgreater 0}$ equipped with $\xi_{m}$ as column vectors is called a conserved matrix for $(\mathcal{X,C,R,K})$. *

Remark 3.6**.**

*For a mass-conserved MAS, the conserved matrix is not unique, but all of them are full column rank. Under the influence of the same conserved matrix $C$, the state of the MAS in $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\geq 0}$ is constrained by this matrix and the equilibrium point $x^{*}$ according to $C^{\top}x=C^{\top}x^{*}$. For a MAS not admitting the mass conservation law, the conserved matrix does not exist. In this case, $q=0$. *

Based on the conserved matrix, we can be aware of the number of non-free variables in a mass-conserved MAS.

Theorem 3.7**.**

*For a mass-conserved MAS $(\mathcal{X,C,R,K})$ described by Eq. 3, $C\in\mathbb{R}^{n\times q}_{\textgreater 0}$ $(q\textgreater 0)$ is a conserved matrix and $x^{*}$ is an equilibrium point, then in $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\geq 0}$ there are at least $q$ non-free state variables, and moreover, these variables are uniquely determined by other $n-q$ state variables and $x^{*}$ but independent of the selection of $C$. *

Proof 3.8**.**

For this mass-conserved MAS, $\forall x\in\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\geq 0}$, we have $C^{\top}x=C^{\top}x^{*}$. Note that rank$(C)=q$ ($1\leq q\leq\mathrm{dim}\left(\mathrm{Ker}(S^{\top})\right)$), so there exist at least $q$ non-free state variables among all $n$ variables contained in the MAS. For simplicity but without loss of generality, the last $q$ state variables from $x_{n-q+1}$ to $x_{n}$ are thought as non-free variables, then we get

[TABLE]

where $C_{l}^{\top}\in\mathbb{R}^{q\times(n-q)}_{>0}$, $C_{r}^{\top}\in\mathbb{R}^{q\times q}_{>0}$ are two sub-blocks of $C^{\top}$, $x_{\bot}=(x_{1},\cdots,x_{n-q})^{\top}$ and $x_{\top}=(x_{n-q+1},\cdots,x_{n})^{\top}$. Since $C$ is full column rank, $C_{r}^{\top}$ is invertible. Also, utilizing $C^{\top}x^{*}=C_{l}^{\top}x_{\bot}^{*}+C_{r}^{\top}x_{\top}^{*}$, the last equality of Eq. 9 may be rewritten as

[TABLE]

which indicates that $x_{\top}$ is uniquely determined by $C$, $x^{*}$ and $x_{\bot}$.

Since the conserved matrix $C$ is not unique for the MAS, we need to further prove the above equality independent of $C$. Assume that $C^{\prime}\in\mathbb{R}^{n\times q}_{\textgreater 0}$ is another conserved matrix for the MAS under consideration, which has the same dimension as $C$. Likewise, we utilize the conservation relation ${C^{\prime}}^{\top}x={C^{\prime}}^{\top}x^{*}$, and further partition $C^{\prime}$ according to the same pattern as $C$, then we obtain

[TABLE]

*By combining it with Eq. 10, we get $C_{r}^{-\top}C^{\top}_{l}(x_{\bot}^{*}-x_{\bot})={C^{\prime}}_{r}^{-\top}{C^{\prime}}_{l}^{\top}(x_{\bot}^{*}-x_{\bot}).$ This means that whatever the conserved matrix is, $C_{r}^{-\top}C^{\top}_{l}(x_{\bot}^{*}-x_{\bot})$ is always equal to ${C^{\prime}}_{r}^{-\top}{C^{\prime}}_{l}^{\top}(x_{\bot}^{*}-x_{\bot})$. Therefore, $x_{\top}$ is uniquely determined by $x^{*}$ and $x_{\bot}$ but independent of the selection of $C$, which completes the proof. *

Based on this theorem, we can thus address the issue that the generalized pseudo-Helmholtz function has a strict minimum [math] at $x^{*}$ through setting $p=n-q$ in Eq. 8.

Theorem 3.9**.**

For a mass-conserved MAS $(\mathcal{X,C,R,K})$ captured by Eq. 3, let $x^{*}$ be an equilibrium, and $C\in\mathbb{R}^{n\times q}_{\textgreater 0}$ be a conserved matrix. Then the generalized pseudo-Helmholtz function given by

[TABLE]

has a strict minimum at $x^{*}$ within $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\textgreater 0}$, while

[TABLE]

Proof 3.10**.**

*From Remark 3.2, we have $\tilde{G}(x)\geq 0$ with equality hold if $x_{\bot}=x_{\bot}^{*}\in\mathbb{R}^{n-q}_{\textgreater 0}$ and $x_{\top}\in\mathbb{R}^{q}_{\textgreater 0}$ may be arbitrary in its domain. Also from Theorem 3.7, we get that $x_{\top}$ is uniquely determined by $x_{\bot}^{*}$ in $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\textgreater 0}$, i.e., $x_{\top}=x_{\top}^{*}$. This means that $\tilde{G}(x)$ has a strict minimum at $x^{*}$ in $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\textgreater 0}$. It is straightforward to obtain Eq. 12 if $\tilde{G}(x)$ and $G(x)$ in Eq. 5 are both taken partial derivative with respect to $x$. *

Remark 3.11**.**

*The more the linearly independent conserved vectors are for a mass-conserved MAS, the more sparse the diagonal matrix $\mathrm{diag}(d_{1},\cdots,d_{n-q},0,\cdots,0)$ is, which implies that the dynamical equation Eq. 3 of this MAS is more reduced. However, it is extremely difficult for a MAS to find a maximum linear independence group of all conserved vectors, i.e., maximizing $q$. In fact, this is not necessary, since a portion of linear independence group, such as a conserved vector composed of the molecular weight of all species, is enough to work as a conserved matrix to derive nice properties of the generalized pseudo-Helmholtz function, like strict convexity at equilibrium. *

Remark 3.12**.**

*If the MAS does not support the mass conservation law, the generalized pseudo-Helmholtz function reduce to the naive extension of the pseudo-Helmholtz function, i.e., $G_{e}(x)$ in Eq. 6 which, obviously, can behave as a Lyapunov function. *

The generalized pseudo-Helmholtz function, as discussed in Remark 3.2, is nonnegative and has lower bound only at the equilibrium, so it may be a Lyapunov function. However, it is still needed the time derivative of $\tilde{G}(x)$ to follow $\dot{\tilde{G}}(x)\leq 0$, i.e., we should have

[TABLE]

where the second equality is based on Eq. 12. As one might know, $G(x)$ is an available Lyapunov function for complex balanced MASs, and its time derivative follows $\dot{G}(x)=\nabla^{\top}G(x)Sv(x)$. However, in the result of Eq. 13 there exists a diagonal matrix between $\nabla^{\top}G(x)$ and $Sv(x)$. In order to fully utilize the performance of $G(x)$ as a Lyapunov function, we manage to eliminate the diagonal matrix $\mathrm{diag}(d_{1},\cdots,d_{n-q},0,\cdots,0)$ through defining reconstruction notation in the next section.

4 Reconstruction

In this section, the reconstruction notation will be defined for characterizing stability of equilibria in some MASs.

4.1 Definition

The reconstruction concept has some relevance to the realization notation proposed in [25]. We firstly revisit the realization definition.

Definition 4.1** (Realization).**

*Given a nonnegative autonomous polynomial system $\dot{x}=f(x)$ with $f:\mathbb{R}^{n}_{>0}\rightarrow\mathbb{R}^{n}$ to be locally Lipschitz, if there exists a MAS $(\mathcal{X,C,R,K})$ modeled by $\dot{x}=Sv(x)$ such that $f(x)=Sv(x)$, then the autonomous polynomial system is kinetically realizable, and the MAS $(\mathcal{X,C,R,K})$ is termed as a kinetic realization, or realization for short, of this polynomial system. *

We then give the definition of reconstruction following that of realization.

Definition 4.2**.**

For a mass-conserved MAS $(\mathcal{X,C,R,K})$ described by $\dot{x}=Sv(x)$, let $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$ be an equilibrium and $C\in\mathbb{R}^{n\times q}_{\textgreater 0}$ be a conserved matrix. If there exists a positive diagonal matrix $D_{1}=\mathrm{diag}(d_{1},\cdots,d_{n-q})$ and a MAS $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ given by $\dot{\hat{x}}=\hat{S}\hat{v}(\hat{x})$ (the parameters follow the same symbols as for $(\mathcal{X,C,R,K})$ but with a hat above) such that

[TABLE]

where

[TABLE]

*then $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ is called a reconstruction of $(\mathcal{X,C,R,K})$, and $D$ is termed as the reconstructing matrix from $(\mathcal{X,C,R,K})$ to $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$. *

Note that the state variables $\hat{x}$ in the reconstructed system is actually the projection of $x$ onto the first $n-q$ coordinations, and Eq. 14 holds for all $x$ in $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\geq 0}$. Some examples of reconstructions are given in the subsequent Table 1.

Remark 4.3**.**

The relation Eq. 14 can be equivalently expressed as

[TABLE]

where $S_{1}\in\mathbb{Z}^{(n-q)\times r}$ comes from $S=(S_{1}^{\top},S_{2}^{\top})^{\top}$. In addition, since $D_{1}$ and $C_{r}^{\top}$ are both invertible, the reconstructing matrix $D$ is invertible with the inverse matrix to be

[TABLE]

*The nonsingular property of $D$ is quite important that will be used to eliminate the diagonal matrix $\mathrm{diag}(d_{1},\cdots,d_{n-q},0,\cdots,0)$ in the result of $\dot{\tilde{G}}(x)$ in Eq. 13. *

Remark 4.4**.**

*If the MAS $(\mathcal{X,C,R,K})$ dose not support the mass conservation law, i.e., $q=0$, then $D=D_{1}=\mathrm{diag}(d_{1},\cdots,d_{n})$. In this case, the reconstruction means to make a positive diagonal transformation from the right hand side of the differential equation $\dot{x}=Sv(x)$ to that of $\dot{\hat{x}}=\hat{S}\hat{v}(\hat{x})$. At this time, the reconstruction is the same with the well-known linear conjugacy concept [14, 15, 16] which defines a linear transformation $\hat{x}=Dx$, and then by substituting $x=D^{-1}\hat{x}$ into $\dot{x}=Sv(x)$ and dynamical equivalence to $\dot{\hat{x}}=\hat{S}\hat{v}(\hat{x})$ produces the relation $DSv(D^{-1}\hat{x})=\hat{S}\hat{v}(\hat{x})$. Further, if $\forall i,~{}d_{i}=1$, then the reconstruction will have the same effect as the realization. *

From Definition 4.2, apparently, for a MAS even though under the same reconstructing matrix, there may exist a few parametrically and/or structurally different reconstructions that are dynamically equivalent. This property is the same as in the realization concept. However, unlike realization the notion of reconstruction only characterizes a binary relation between MASs, but not between the MAS and a general polynomial system. Therefore, there are at least two MASs involved when we mention reconstruction. The whole process seems to reconstruct a new MAS from a known one. The reconstruction notation together with the generalized pseudo-Helmholtz function will play an important role on capturing the stability of some MASs.

4.2 Reconstruction for characterizing stability

In this subsetion, we will characterize stability of equilibria in MASs utilizing the reconstruction strategy.

Lemma 4.5**.**

*Assume that $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ is a reconstruction of a mass-conserved MAS $(\mathcal{X,C,R,K})$ under the reconstructing matrix $D$ given by Eq. 15, and their dynamical equations are $\dot{\hat{x}}=\hat{S}\hat{v}(\hat{x})$ and Eq. 3, respectively. If $x^{*}\in\mathbb{R}^{n}_{>0}$ is an equilibrium in $(\mathcal{X,C,R,K})$, then $\hat{x}^{*}=x_{\bot}^{*}\in\mathbb{R}^{n-q}_{>0}$ is an equilibrium in $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$. *

Proof 4.6**.**

*Since $x^{*}\in\mathbb{R}^{n}_{>0}$ is an equilibrium in $(\mathcal{X,C,R,K})$, we get $Sv(x^{*})=\mathbbold{0}_{n}$, i.e., $D_{1}S_{1}v(x^{*})=\mathbbold{0}_{n-q}$. From Eq. 16, we have $D_{1}S_{1}v(x^{*})=\hat{S}\hat{v}(x_{\bot}^{*})$. Therefore, $\hat{x}^{*}=x_{\bot}^{*}$ is an equilibrium in $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$. *

Theorem 4.7**.**

*For a mass-conserved MAS $(\mathcal{X,C,R,K})$ modeled by Eq. 3 with an equilibrium $x^{*}\in\mathbb{R}^{n}_{>0}$, assume that the MAS $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ in the form of $\dot{\hat{x}}=\hat{S}\hat{v}(\hat{x})$ is its reconstruction, bridged by the reconstructing matrix $D$ defined in Eq. 15, then $\hat{x}^{*}\in\mathbb{R}^{n-q}_{>0}$ is an equilibrium in $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$. If $\hat{x}^{*}$ is stable rendered by a Lyapunov function $\hat{V}:\mathbb{R}_{\geq 0}^{n-q}\to\mathbb{R}_{\geq 0}$ with $\hat{V}\in\mathscr{C}^{2}$ and its Hessian matrix to be a $(n-q)$-dimensional positive diagonal matrix, then $x^{*}$ is stable in $(\mathcal{X,C,R,K})$. *

Proof 4.8**.**

According to Lemma 4.5, $\hat{x}^{*}$ is naturally an equilibrium in $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$. As $\hat{V}\in\mathscr{C}^{2}$ is a Lyapunov function to render the system $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ stable at $\hat{x}^{*}$, we have $\hat{V}(\hat{x})\geq 0$ with equality holding if and only if $\hat{x}=\hat{x}^{*}$, and moreover, $\nabla\hat{V}(\hat{x}^{*})=\mathbbold{0}_{n-q}$ and $\forall~{}\hat{y}\in\mathbb{R}^{n-q}_{\textgreater 0}$, if $\hat{y}\neq\hat{x}^{*}$, then $\nabla\hat{V}(\hat{y})\neq\mathbbold{0}_{n-q}$ and $\partial\hat{V}(\hat{y})/\partial\hat{x}_{i}$ is strictly increasing. In addition,

[TABLE]

where $D_{1}$ and $S_{1}$ are given according to Eq. 16. Construct the following $\mathscr{C}^{2}$ function $V:\mathbb{R}^{n}_{\geq 0}\to\mathbb{R}_{\geq 0}$ as

[TABLE]

where $D_{2}=[D_{1},\mathbbold{0}_{(n-q)\times q}]$ and $y\in\mathbb{R}^{n}_{\geq 0}$ with $\hat{y}$ to be the projection of $y$ onto the first $n-q$ coordinates. Note that $x^{*}$ and $x$ are vectors, so the integral is over the line connecting these two points.

It is clear that $V(x^{*})=0$ and $\nabla V(x^{*})=D_{2}^{\top}\nabla\hat{V}(\hat{x}^{*})=\mathbbold{0}_{n}$. Utilizing the mean value theorem of integrals, we can rewrite Eq. 19 as

[TABLE]

where $\bar{x}\in\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\textgreater 0}$ lies between $x^{*}$ and $x$, and $\hat{\bar{x}}$ is the projection of $\bar{x}$ onto the first $n-q$ coordinates. Since $\frac{\partial\hat{V}(\hat{\bar{x}})}{\partial\hat{x}_{i}}$ is strictly increasing and $\nabla\hat{V}(\hat{x}^{*})=\mathbbold{0}_{n-q}$, we have $\frac{\partial\hat{V}(\hat{\bar{x}})}{\partial\hat{x}_{i}}(x_{i}-x_{i}^{*})\geq 0$. Therefore, $V(x)=0$ if and only if $x_{\bot}=x_{\bot}^{*}$, but $x_{\top}$ may be arbitrary, where $x_{\bot}=(x_{1},\cdots,x_{n-q})^{\top}$ and $x_{\top}=(x_{n-q+1},\cdots,x_{n})^{\top}$. Also since $x_{\top}$ is uniquely determined by $x_{\bot}^{*}$ in $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\textgreater 0}$, we have $V(x)=0$ iff $x=x^{*}$. Again from Eq. 19 we can write the Hessian matrix of $V(x)$ to be

[TABLE]

Since the Hessian matrix of $\hat{V}(\hat{x})$ is a positive diagonal matrix, $\forall~{}x\in\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\textgreater 0}$ the Hessian matrix $H(V(x))$ is semi-positive definite. This means that $V(x)$ is a convex function in $\mathscr{S}(x^{*})\cap\mathbb{R}^{n}_{\textgreater 0}$, and $x^{*}$ is a minimum point for $V(x)$, i.e., $V(x)\geq V(x^{*})$. Therefore, we have $V(x)\geq 0$ with equality hold iff $x=x^{*}$. Further, we have

[TABLE]

*where the last inequality holds from Eq. 18. Moreover, when $x=x^{*}$, we get $\frac{dV(x)}{dt}=0$. Hence, $x^{*}$ is a stable equilibrium in the original MAS $(\mathcal{X,C,R,K})$ with $V(x)$ as the Lyapunov function. *

Corollary 4.9**.**

*Consider a mass-conserved MAS $(\mathcal{X,C,R,K})$ described by Eq. 3 with $x^{*}\in\mathbb{R}^{n}_{>0}$ as an equilibrium. Let $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ be a reconstruction. If $\hat{x}^{*}\in\mathbb{R}^{n-q}_{>0}$ is a complex balanced equilibrium in the reconstructed system , then $x^{*}$ is stable in the original MAS. *

Proof 4.10**.**

*From Theorem 2.10, the complex balanced system $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ is locally asymptotically stable at $\hat{x}^{*}$ with $G(\hat{x})$ defined like Eq. 5 as the Lyapunov function. Note that $G\in\mathscr{C}^{2}$ and $\frac{\partial G^{2}(\hat{x})}{\partial\hat{x}^{2}}=\mathrm{diag}(1/\hat{x}_{1},\cdots,1/\hat{x}_{n-q})$, which satisfy the conditions requested by Theorem 4.7, therefore, $(\mathcal{X,C,R,K})$ is stable at $x^{*}$. *

In the following special case, the stability properties can be transferred mutually between the original network and the reconstruction.

Proposition 4.11**.**

*Consider a MAS $(\mathcal{X,C,R,K})$ with the linear kinetics $\dot{x}=Px$, where $P\in\mathbb{R}^{n\times n}$ is a Metzler matrix. Assume that $x^{*}\in\mathbb{R}^{n}_{>0}$ is an equilibrium in $(\mathcal{X,C,R,K})$ and $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ is a reconstruction of $(\mathcal{X,C,R,K})$ under a positive diagonal matrix $D$, then $(\mathcal{X,C,R,K})$ is stable at $x^{*}$ if and only if $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ is stable at $\hat{x}^{*}\in\mathbb{R}^{n}_{>0}$. *

Proof 4.12**.**

*Clearly, $\hat{x}^{*}$ is an equilibrium in $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$. Note that the reconstruction $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ induces a linear kinetic dynamics $\dot{\hat{x}}=DP\hat{x}$, in which $DP$ is also a Metzler matrix. Based on the $D$-stability of matrices [19], the system $(\mathcal{X,C,R,K})$ is stable at $x^{*}$ if and only if $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ is stable at $\hat{x}^{*}$. *

As can be seen from Corollary 4.9 and Proposition 4.11, the stability properties are retained from the reconstruction to the original network usually requesting the former to have special network structure. As far as the current technique is concerned, it is difficult to say that a stable reconstruction with a general structure must lead to a stable original network. We place this issue as a possible point of future research. The high validity of complex balancing in the reconstruction for stability retainment motivates us to further use this strategy to characterize asymptotic stability of the original network.

5 Asymptotic stability characterized through complex balanced reconstruction

In Theorem 4.7, we are able to characterize stability of equilibrium points in MASs, but fail to capture local asymptotic stability. Here, the stability of an equilibrium means that the system trajectory runs in a bounded ball centered around this equilibrium while the local asymptotic stability means that the system trajectory will asymptotically converge to this equilibrium if the trajectory starts from a neighborhood of the equilibrium. The main reason is that there may exist multiple equilibria in the MAS so that the time derivative of the Lyapunov function is negative semidefinite with respect to the concerning equilibrium. In this section, we will address this issue to characterize the local asymptotic stability of some MASs through setting the reconstruction to have a complex balancing structure. We firstly introduce a special reconstruction.

5.1 Reverse reconstruction

Reverse reconstruction is a generalized MAS [13], obtained from the reconstruction for a given MAS, which is defined as follows.

Definition 5.1**.**

Given a mass-conserved MAS $(\mathcal{X,C,R,K})$, let $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ denote its reconstruction under the reconstructing matrix $D$ defined in Eq. 15. A MAS $(\mathcal{\tilde{X},\tilde{C},\tilde{R},\tilde{K}})$ is the reverse reconstruction of $(\mathcal{X,C,R,K})$ with respect to $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ if its complexes set $\tilde{\mathcal{C}}=\bigcup_{j=1}^{\tilde{r}}\{\tilde{Z}_{\cdot j},\tilde{Z}^{\prime}_{\cdot j}\}$ and reactions set \tilde{\mathcal{R}}=\bigcup_{j=1}^{\tilde{r}}\{\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern 7.8689pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern 0.0pt\offinterlineskip\halign{\entry@#!@&&\entry@@#!@\cr&\crcr}}}\ignorespaces{\hbox{\kern-7.8689pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{\textstyle{\tilde{Z}{\cdot j}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 11.32698pt\raise 7.65953pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.56271pt\hbox{\scriptstyle{\tilde{k}{j}~{}}}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 31.8689pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 31.8689pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{\textstyle{\tilde{Z}^{\prime}_{\cdot j}}}}}}}}}\ignorespaces}}}}\ignorespaces\} satisfy

• (i)

$\tilde{r}=\hat{r}$;

• (ii)

$\tilde{Z}_{\cdot j}=\hat{Z}_{\cdot j},~{}\tilde{Z}^{\prime}_{\cdot j}=\hat{Z}_{\cdot j}+D_{1}^{-1}(\hat{Z}^{\prime}_{\cdot j}-\hat{Z}_{\cdot j}),~{}\tilde{k}_{j}=\hat{k}_{j},~{}\forall~{}j=1,\cdots,\tilde{r}$.

Remark 5.2**.**

The dynamics of the reverse reconstruction $(\mathcal{\tilde{X},\tilde{C},\tilde{R},\tilde{K}})$ follows

[TABLE]

*where $\tilde{S}\in\mathbb{R}^{(n-q)\times\tilde{r}}$, $\tilde{v}(\cdot)$ is a $\tilde{r}$-dimensional vector valued function, and $\tilde{x}$ is the projection of $\hat{x}$ onto itself. It is clear from the reverse reconstruction definition that $\tilde{S}=D_{1}^{-1}\hat{S}$ and $\tilde{v}(\tilde{x})=\hat{v}(\hat{x})$. *

Remark 5.3**.**

*The underlying CRN in the reverse reconstruction does not fit into the CRN class defined in Definition 2.1, since the elements in $\tilde{Z}^{\prime}_{\cdot j}$ are not necessarily nonnegative integer. However, this change does not invalidate the subsequent results. The main reason is that the reverse reconstruction belongs to generalized MASs, but some results of CRN theory (used for the subsequent analysis)[9], like power-law rate functions, definitions of stoichiometric subspace and stoichiometric compatibility classes, carry over to the case of generalized mass action kinetics[13, 18]. *

Proposition 5.4**.**

*For a mass-conserved MAS $(\mathcal{X,C,R,K})$ with dynamics Eq. 3 and an equilibrium $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$, let $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$ be its reconstruction bridged by the reconstructing matrix $D$ defined in Eq. 15, and $(\mathcal{\tilde{X},\tilde{C},\tilde{R},\tilde{K}})$ be its reverse reconstruction with respect to $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$. Then $\tilde{x}^{*}=x^{*}_{\bot}\in\mathbb{R}^{n-q}_{\textgreater 0}$ is an equilibrium in $(\mathcal{\tilde{X},\tilde{C},\tilde{R},\tilde{K}})$, and the dynamics of $(\mathcal{\tilde{X},\tilde{C},\tilde{R},\tilde{K}})$ is equivalent to that of $(\mathcal{X,C,R,K})$ at the first $n-q$ species. *

Proof 5.5**.**

Remark 12 reveals that $\tilde{S}=D_{1}^{-1}\hat{S}$ and $\tilde{v}(\tilde{x})=\hat{v}(\hat{x})$. Hence, $D_{1}^{-1}\hat{S}\hat{v}(\hat{x})=\tilde{S}\tilde{v}(\tilde{x})$. Further, utilizing Eq. 14 we get

[TABLE]

*where $E_{n-q}$ is a $(n-q)$-dimensional unit matrix. Therefore, $S_{1}v(x)=\tilde{S}\tilde{v}(\tilde{x})$. Note that the left term measures the dynamics of $(\mathcal{X,C,R,K})$ at the first $n-q$ species, i.e., $\dot{x}_{\bot}(t)=S_{1}v(x)$. We thus obtain $\dot{x}_{\bot}(t)=\dot{\tilde{x}}(t)$. $\tilde{S}\tilde{v}(\tilde{x}^{*})=\mathbbold{0}_{n-q}$ is straightforward since $Sv(x^{*})=\mathbbold{0}_{n}$. Therefore, the results are true. *

Corollary 5.6**.**

*Assume that $(\mathcal{X,C,R,K})$, $(\mathcal{\hat{X},\hat{C},\hat{R},\hat{K}})$, $(\mathcal{\hat{X},\tilde{C},\tilde{R},\tilde{K}})$ and $D$ have the same meanings with those in Proposition 5.4. If $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$ is an equilibrium in $(\mathcal{X,C,R,K})$, then the reconstruction and reverse reconstruction systems have the common equilibrium point $\hat{x}^{*}=\tilde{x}^{*}=x_{\bot}^{*}$, and moreover, $\mathrm{Im}(\tilde{S})=\mathrm{Im}(D_{1}^{-1}\hat{S})$. *

Reverse reconstruction is a very special MAS, which has the same dynamics with that of the original MAS at the first $n-q$ species. Therefore, if an equilibrium point in the original network system is locally asymptotically stable whose first $n-q$ entries are naturally locally asymptotically stable, then this stability behavior can be potentially characterized based on an similar analysis in the reverse reconstruction and on how stability properties are retained under the adopted nonsingular transformation. We tackle these tasks in the following by using a particular network class for reconstruction, i.e., so-called complex balanced MASs, for which the local asymptotic stability has been proved in Theorem 2.10.

5.2 Complex balanced reconstruction

As said in Corollary 4.9, when a mass-conserved MAS has a complex balanced MAS as a reconstruction, then its equilibrium must be stable. We thus set the reconstruction to be a complex balanced system for further analysis. This kind of reconstructions is referred to as the complex balanced reconstruction in the context. For simplicity of symbols, we denote the set of MASs admitting complex balanced reconstructions by $\mathscr{N}_{\mathrm{Com}}$, i.e., for any mass-conserved MAS $(\mathcal{X,C,R,K})\in\mathscr{N}_{\mathrm{Com}}$, there exist a complex balanced MAS, denoted by $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}}_{C})$, and a reconstructing matrix $D$ defined in Eq. 15 satisfying

[TABLE]

where $\hat{S}_{C}\in\mathbb{Z}^{(n-q)\times\hat{r}_{C}}$.

In the case of existence of multiple equilibria in a MAS, the key to extend stability of equilibrium to asymptotic stability is to identify that there is only/at most an equilibrium in each positive stoichiometric compatibility class. The main reason is that the state of the MAS will evolve only in a stoichiometric compatibility class, as given in Remark 2.5. We try to reach this conclusion for the reverse reconstruction firstly, since its dynamics links that of the original MAS in $\mathscr{N}_{\mathrm{Com}}$ by a constant matrix.

As one may know, there exists a close relation between equilibria for a complex balanced MAS, shown as follows.

Lemma 5.7** ([10, 13]).**

For any complex balanced MAS $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}}_{C})$ admitting an equilibrium $\hat{x}^{*}\in\mathbb{R}^{n-q}_{>0}$, the concentration vector $\hat{x}^{{\dagger}}\in\mathbb{R}^{n-q}_{>0}$ is another equilibrium if and only if

[TABLE]

Remark 5.8**.**

*This lemma suggests that if a complex balanced MAS is conservative then it must admit multiple equilibrium points. However, mass conservation law is not necessary to say the existence of multiple equilibria in a complex balanced system. It is still possible to possess multiple equilibria for the complex balanced MAS even if mass conservation does not hold. Of course, whatever the case is, if there are multiple equilibria in a complex balanced MAS, then each equilibrium is in its individual positive stoichiometric compatibility class. *

We reach the following proposition to approach the result that each positive stoichiometric compatibility class in the reverse reconstruction of a MAS $(\mathcal{X,C,R,K})\in\mathscr{N}_{\mathrm{Com}}$ contains only an equilibrium.

Proposition 5.9**.**

For a mass-conserved MAS $(\mathcal{X,C,R,K})$ admitting an equilibrium $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$, let $(\mathcal{\hat{X}},\hat{\mathcal{C}},\hat{\mathcal{R}},\hat{\mathcal{K}})$ be its reconstruction bridged by the reconstructing matrix $D$ given in Eq. 15, i.e., satisfying Eq. 14. Assume that $x_{0}\in\mathbb{R}^{n}_{\textgreater 0}$ represents any initial state of $(\mathcal{X,C,R,K})$, then for $\forall~{}\hat{x}_{0}\in\mathbb{R}^{n-q}_{\textgreater 0}$ there exists a unique $\mu\in\mathrm{Ker}(\hat{S}^{\top})$ such that

[TABLE]

*where $\hat{x}^{\ast}$ is the vector of the first $n-q$ entries in $x^{*}$. *

Proof 5.10**.**

Clearly, $\hat{x}^{*}$ is an equilibrium in $(\mathcal{\hat{X}},\hat{\mathcal{C}},\hat{\mathcal{R}},\hat{\mathcal{K}})$, we thus define a function $\varphi_{D_{1}}:\mathbb{R}^{n-q}\mapsto\mathbb{R}$ as

[TABLE]

Following the almost completely same proof as in Appendix B of [10] where the function $\sum^{n}_{i=1}[\hat{x}_{i}^{\ast}\mathrm{exp}(u_{i})-\hat{x}_{0_{i}}u_{i}]$ was proved strictly convex, $\varphi_{D_{1}}(u)$ in Eq. 23 can be proved strictly convex and $\lim_{a\rightarrow\infty}\varphi_{D_{1}}(au)=\infty$ for any $u\neq\mathbbold{0}_{n-q}$. Note that the gradient of $\varphi_{D_{1}}(u)$ takes

[TABLE]

and $\mathrm{Ker}(\hat{S}^{\top})\subset\mathbb{R}^{n-q}$. In the same way as in [10], we can prove that there exists a unique $\mu\in\mathrm{Ker}(\hat{S}^{\top})$ such that $\nabla\varphi_{D_{1}}^{\top}(\mu)u=0$ for any $u\in\mathrm{Ker}(\hat{S}^{\top})$, so $\nabla\varphi_{D_{1}}(\mu)\in\mathrm{Im}(\hat{S})$. Namely, for any $\hat{x}_{0}\in\mathbb{R}^{n-q}_{\textgreater 0}$ there exists a unique $\mu\in\text{Ker}(\hat{S}^{\top})$ such that

[TABLE]

*which completes the proof. *

We then give an important result about the reverse reconstruction for a MAS in $\mathscr{N}_{\mathrm{Com}}$.

Theorem 5.11**.**

*Consider a mass-conserved MAS $(\mathcal{X,C,R,K})\in\mathscr{N}_{\mathrm{Com}}$ governed by Eq. 3 with an equilibrium $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$. Let $x_{0}\in\mathbb{R}^{n}_{>0}$ represent any of its initial state, $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}}_{C})$ be a complex balanced reconstruction with $D$ given by Eq. 15 as the reconstructing matrix, and $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$ be the corresponding reverse reconstruction. Then $\forall~{}\tilde{x}_{0}\in\mathbb{R}^{n-q}_{>0}$ there exists a unique equilibrium $\tilde{x}^{\dagger}$ in $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$ such that $\tilde{x}^{\dagger}\in\tilde{\mathscr{S}}_{C}(\tilde{x}_{0})\cap\mathbb{R}^{n-q}_{>0}$. *

Proof 5.12**.**

*According to Proposition 5.9, for any $\hat{x}_{0}=x_{0_{\bot}}\in\mathbb{R}^{n-q}_{>0}$ there exists a unique $\mu\in\mathrm{Ker}(\hat{S}_{C}^{\top})$ such that $D_{1}[\hat{x}^{\ast}\cdot\mathrm{Exp}(\mu)-\hat{x}_{0}]\in\text{Im}(\hat{S}_{C})$, which together with Corollary 5.6 yields $[\tilde{x}^{\ast}\cdot\mathrm{Exp}(\mu)-\tilde{x}_{0}]\in\text{Im}(\tilde{S}_{C})$. Note that $\hat{x}^{\ast}$ and $\tilde{x}^{\ast}$ are an equilibrium in $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}}_{C})$ and $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$, respectively. Set $\hat{x}^{\dagger}=\hat{x}^{\ast}\cdot\mathrm{Exp}(\mu)$, which naturally satisfies $\hat{x}^{\dagger}\in\mathbb{R}^{n-q}_{>0}$, then we get $\mu=\text{Ln}(\frac{\hat{x}^{\dagger}}{\hat{x}^{\ast}})\in\text{Ker}(\hat{S}_{C}^{\top})$. Further based on Lemma 5.7, we obtain that $\hat{x}^{\dagger}$ is an equilibrium in $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}}_{C})$, so $\tilde{x}^{\dagger}$ is an equilibrium in $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$. Therefore, $\tilde{x}^{\dagger}\in\tilde{\mathscr{S}}_{C}(\tilde{x}_{0})\cap\mathbb{R}^{n-q}_{>0}$. The uniqueness of $\mu$ guarantees that $\tilde{x}^{\dagger}$ is a unique equilibrium in $\tilde{\mathscr{S}}_{C}(\tilde{x}_{0})\cap\mathbb{R}^{n-q}_{>0}$. *

Remark 5.13** (Robustness of the first $n-q$ species at equilibrium).**

*This theorem means that the reverse reconstruction of a mass-conserved MAS $(\mathcal{X,C,R,K})\in\mathscr{N}_{\mathrm{Com}}$ admitting an equilibrium $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$ and starting from any initial sate $x_{0}\in\mathbb{R}^{n}_{\textgreater 0}$ must contain an equilibrium, and moreover, contains precisely an equilibrium $\tilde{x}^{\dagger}\in\mathbb{R}^{n-q}_{\textgreater 0}$ in each positive stoichiometric compatibility class. However, it should be noted that $\tilde{x}^{\dagger}$ does not guarantee to map an equilibrium in the original system, since from the fact $S_{1}v(x^{\dagger})=\tilde{S}\tilde{v}(\tilde{x}^{\dagger})=\mathbbold{0}_{n-q}$ the last $q$ species in $x^{\dagger}$ can not guarantee to be positive despite being arbitrary. If there exist multiple equilibria in $(\mathcal{X,C,R,K})$ mapped from $\tilde{x}^{\dagger}$, the concentrations of the first $n-q$ species at these different equilibria will be robust with respect to the last $q$ species in the initial state $x_{0}$, i.e., robust with respect to $x_{0_{\top}}$. Specially, if $\mathrm{dim}\left(\mathrm{Im}(\hat{S}_{C})\right)=n-q$, then $\mu=\mathbbold{0}_{n-q}$. At this time, $x^{\dagger}_{\bot}=x^{*}_{\bot}$, which indicates the concentrations of the first $n-q$ species will keep unchanged at all of different equilibria, i.e., absolute concentration robustness[22]. *

5.3 Asymptotic stability characterized for MASs in $\mathscr{N}_{\mathrm{Com}}$

Now we can characterize the local asymptotic stability of any MAS in $\mathscr{N}_{\mathrm{Com}}$ utilizing the complex balanced reconstruction strategy.

Theorem 5.14**.**

*For any mass-conserved MAS $(\mathcal{X,C,R,K})\in\mathscr{N}_{\mathrm{Com}}$ described by Eq. 3, assume that $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}_{C}})$ is its complex balanced reconstruction with $D$ given by Eq. 15 as the reconstructing matrix, and $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$ is the corresponding reverse reconstruction. Let $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$ be an equilibrium in $(\mathcal{X,C,R,K})$, then $\tilde{x}^{*}\in\mathbb{R}^{n-q}_{\textgreater 0}$ is locally asymptotically stable with respect to all initial conditions in $\tilde{\mathscr{S}}_{C}(\tilde{x}^{*})\cap\mathbb{R}^{n-q}_{>0}$ near $\tilde{x}^{*}$. *

Proof 5.15**.**

From $S_{1}v(x^{*})=\tilde{S}\tilde{v}(\tilde{x}^{*})$, $\tilde{x}^{*}$ is an equilibrium in $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$. Further based on Theorem 5.11, $\tilde{x}^{*}$ is the unique equilibrium in $\tilde{\mathscr{S}}_{C}(\tilde{x}^{*})\cap\mathbb{R}^{n-q}_{>0}$. For the reverse reconstruction $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$ , take the generalized pseudo-Helmholtz function $\tilde{G}(\tilde{x})$ as the Lyapunov function, which is naturally positive definite with respect to $\tilde{x}-\tilde{x}^{*}$, as said in Theorem 3.9. Further, we have its time derivative as

[TABLE]

*Note that the pseudo-Helmholtz function $G(\hat{x})$ is a valid Lyapunov function for a complex balanced MAS to suggest the local asymptotic stability at the equilibrium [10, 13], so $\forall~{}\hat{x}$ we have $\nabla^{\top}G(\hat{x})\hat{S}_{C}\hat{v}_{C}(\hat{x})\leq 0$ with equality hold at equilibria of the complex balanced reconstruction $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}_{C}})$, i.e., at equilibria of the reverse reconstruction $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$. Since it is impossible for $\tilde{\mathscr{S}}_{C}(\tilde{x}^{*})\cap\mathbb{R}^{n-q}_{>0}$ to include other equilibrium of the reverse reconstruction, $\dot{\tilde{G}}(\tilde{x})\leq 0$ with equality hold if and only if $\tilde{x}=\tilde{x}^{*}$ in this positive stoichiometric compatibility class. Therefore, $\tilde{x}^{*}$ is locally asymptotically stable with all initial conditions in $\tilde{\mathscr{S}}_{C}(\tilde{x}^{*})\cap\mathbb{R}^{n-q}_{>0}$ near $\tilde{x}^{*}$. *

Remark 5.16**.**

Theorem 5.14* indicates that for the reverse reconstruction $(\mathcal{\tilde{X}},\tilde{\mathcal{C}}_{C},\tilde{\mathcal{R}}_{C},\tilde{\mathcal{K}}_{C})$, $\forall~{}\epsilon,\epsilon^{\prime}\textgreater 0$, $\exists~{}\delta(\epsilon)\textgreater 0$ such that every solution of $\dot{\tilde{x}}=\tilde{S}_{C}\tilde{v}_{C}(\tilde{x})$ having initial condition within distance $\delta(\epsilon)$, i.e., $\|\tilde{x}_{0}-\tilde{x}^{*}\|\textless\delta(\epsilon)$, and $\tilde{x}_{0}\in\tilde{\mathscr{S}}_{C}(\tilde{x}^{*})\cap\mathbb{R}^{n-q}_{>0}$, the equilibrium $\tilde{x}^{*}$ remains within distance $\epsilon$, i.e., $\|\tilde{x}-\tilde{x}^{*}\|\textless\epsilon$ for all $t\geq 0$ and $\lim_{t\rightarrow{\infty}}\|\tilde{x}(t)-\tilde{x}^{*}\|\leq\epsilon^{\prime}$. *

Theorem 5.17**.**

*For a mass-conserved MAS $(\mathcal{X,C,R,K})\in\mathscr{N}_{\mathrm{Com}}$, let $x^{*}\in\mathbb{R}^{n}_{\textgreater 0}$ be any equilibrium, then it is locally asymptotically stable. *

Proof 5.18**.**

Assume $x_{0}\in\mathbb{R}^{n}_{\textgreater 0}$ to be an initial point for the MAS $(\mathcal{X,C,R,K})$ such that $\tilde{x}_{0}=x_{0_{\bot}}\in\tilde{\mathscr{S}}_{C}(\tilde{x}^{*})\cap\mathbb{R}^{n-q}_{>0}$. From Theorem 5.14 there must exist a neighbourhood of $\tilde{x}^{*}$ so that if $\tilde{x}_{0}$ is in this neighbourhood then $\tilde{x}^{*}$, i.e., $x^{*}_{\bot}$, is locally asymptotically stable. This means that, based on Remark 5.16 and Proposition 5.4, $\forall~{}\epsilon,\epsilon^{\prime}\textgreater 0$, $\exists~{}\delta(\epsilon)\textgreater 0$, when $\|x_{0_{\bot}}-x^{*}_{\bot}\|\textless\delta(\epsilon)$ we have $\|x_{\bot}-x^{*}_{\bot}\|\textless\epsilon$ for all $t\geq 0$ and $\lim_{t\rightarrow{\infty}}\|x_{\bot}(t)-x^{*}_{\bot}\|\leq\epsilon^{\prime}$. From Eq. 10, we get

[TABLE]

Hence,

[TABLE]

*Note that $\sqrt{1+\|{C}_{r}^{-\top}{C}_{l}^{\top}\|^{2}}$ is a constant, so we have $\forall~{}\epsilon_{1}=\sqrt{1+\|{C}_{r}^{-\top}{C}_{l}^{\top}\|^{2}}\epsilon\textgreater 0$ and $\epsilon_{1}^{\prime}=\sqrt{1+\|{C}_{r}^{-\top}{C}_{l}^{\top}\|^{2}}\epsilon^{\prime}\textgreater 0$, $\exists~{}\delta(\epsilon_{1})\textgreater 0$, when $\|x_{0}-x^{*}\|\textless\delta(\epsilon_{1})$ we have $\|x-x^{*}\|\textless\epsilon_{1}$ for all $t\geq 0$ and $\lim_{t\rightarrow{\infty}}\|x(t)-x^{*}\|\leq\epsilon_{1}^{\prime}$. This implies that the arbitrary equilibrium $x^{*}$ in $(\mathcal{X,C,R,K})$ is locally asymptotically stable. *

Although all of the above results, including the reverse reconstruction possessing an equilibrium in each positive stoichiometric compatibility class, the local asymptotic stability of any MAS in $\mathscr{N}_{\mathrm{Com}}$, etc., are achieved under the condition of mass conservation in the network, they also apply to the case that mass conservation is not admitted by the network. For the latter, the case will become simpler since the reconstructing matrix is a positive diagonal matrix. There will be no dimension reduction in the reconstruction and the corresponding reverse reconstruction compared to the original network. In that case, the reconstruction concept has the same signification as the linear conjugacy proposed in [14, 15, 16], as stated in Remark 4.4. However, there is a large difference in using Lyapunov function to characterize asymptotic stability of the MAS between these two strategies. The current reconstruction strategy uses the extension of the pseudo-Helmholtz function as the Lyapunov function while the latter uses the traditional pseudo-Helmholtz function. Intuitively, the extension of the pseudo-Helmholtz function and the generalized one lay a basis for deeper studies. In the case that the mass conservation law is admitted, these two concepts are completely different. The utilization of the mass conservation law will reduce the complexity of the dynamical model greatly. In this case there will be no definition for the linear conjugacy.

6 Systematic method to find complex balanced reconstructions and cases studies

This section contributes to designing algorithm to compute the complex balanced reconstructions for a MAS in $\mathscr{N}_{\mathrm{Com}}$, and then illustrating the proposed algorithm through some representative examples.

6.1 Algorithm for finding complex balanced reconstructions

Finding complex balanced reconstructions for a MAS $(\mathcal{X,C,R,K})\in\mathscr{N}_{\mathrm{Com}}$ means to design $D,\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C}$ and $\hat{\mathcal{K}}_{C}$ such that

[TABLE]

where $\hat{B}_{C}$ and $\hat{x}^{\ast}$ are the incident matrix and complex balanced equilibrium in the complex balanced MAS $(\mathcal{\hat{X}},\hat{\mathcal{C}}_{C},\hat{\mathcal{R}}_{C},\hat{\mathcal{K}}_{C})$, respectively. Every MAS with the dynamics $\dot{\hat{x}}=\hat{S}_{C}\hat{v}_{C}(\hat{x})$ is a reconstruction needed to be found. Essentially, it is to find realizations with a complex balancing structure for the nonnegative autonomous polynomial system with the vector field of $D_{1}S_{1}v(x)$. Therefore, the known algorithms for finding realizations [24, 25, 26, 16] naturally constitute the base to achieve complex balanced reconstructions.

Note that the dynamics $\dot{x}=Sv(x)$ is expressed in a reaction centered formulation, which is structure-hidden. To reflect the underlying CRN structure, we use complex centered formulation [25] instead of the current one through substituting Eq. 1 into the dynamical equation, and then replacing $Bv(x)$ by $L\varPsi(x)$, where $L\in\mathbb{R}^{c\times c}$ is the Kirchhoff matrix that stores the reaction rate coefficients according to

[TABLE]

with $k_{\pi\rho}(\rho,\pi=1,\cdots,c)$ indicating the reaction rate coefficient of reaction $Z_{\cdot\pi}\to Z_{\cdot\rho}$, and $\varPsi(x)\in\mathbb{R}^{c}$ is a polynomial vector with the $\rho$th entry to be

[TABLE]

Clearly, $L$ is column conservative, constrained by $\mathbbold{1}_{c}^{\top}L=\mathbbold{0}_{c}^{\top}$, and $\varPsi(x)$ can be solely identified by $Z$. In Eq. 25, $k_{\pi\rho}=0$ when $Z_{\cdot\pi}\to Z_{\cdot\rho}$ is not existing. Note that these reaction rate coefficients are identified by two complexes, which are essentially the same with those in Eq. 2 identified by reactions. Utilizing the complex centered formulation, we can rewrite the complex balanced reconstruction condition Eq. 24 as

[TABLE]

where the third constraint serves for eliminating those non-free species in the complex balanced reconstruction. Therefore, the essence of finding complex balanced reconstructions is to find complex balanced realizations, uniquely characterized by the pair $(\hat{Z}_{C},\hat{L}_{C})$, for the nonnegative autonomous polynomial system with the vector field of $DZL\varPsi(x)$ under given $D$.

There exist a few optimization based algorithms to find realizations, like minimizing/maximizing the sum of reaction rate coefficients [25], minimizing/maximizing the number of reactions [24, 26], minimizing the deficiency of weakly reversible realizations [16], etc. Some linear programming and mixed-integer linear programming problems are designed towards these objects. Here, for simplicity we use the strategy of minimizing the sum of reaction rate coefficients to design algorithm for finding complex balanced reconstructions. In addition, we calculate the conserved matrix $C$ for the original MAS $(\mathcal{X,C,R,K})$ beforehand, which may be easily obtained through solving the linear system of equations $S^{\top}\xi=\mathbbold{0}_{r}$, and then following Definition 3.5.

With the above information, we can design the linear programming problem towards solving complex balanced reconstructions for the kinetic system $\dot{x}=ZL\varPsi(x)$ to be

[TABLE]

where $D_{1}=\mathrm{diag}(d_{i})$ and $\hat{L}_{C}$ are the decision variables, and $\epsilon>0$ is a given very small positive number. Once $Z$ and $L$ are given, then $x^{*}$, i.e., $\hat{x}^{*}=x^{*}_{\bot}$, may be calculated based on $ZL{\varPsi}(x^{\ast})=\mathbbold{0}_{n}$. These parameters together with a given $\epsilon$ are known data for Eq. P1. The pair $(\hat{Z}_{C},\hat{L}_{C})$ solved from Eq. P1, namely, represents the desired complex balanced reconstruction.

6.2 Cases studies

In this subsection, some examples of MASs and their complex balanced reconstructions computed through solving the linear programming Eq. P1 are shown. Table 1 reports the computation results of complex balanced reconstructions for $6$ MASs. The first two examples are two-species CRNs, and the third and fourth examples are three-species but having the conserved matrices of rank $1$ and $2$, respectively. The fifth and sixth CRNs contain four species. Clearly, every MAS has a complex balanced reconstruction, and for the first, second, fourth and sixth examples the complex balanced reconstructions are even detailed balanced. The reconstructing matrix in each case is also reported in Table 1.

7 Conclusions

This paper has developed the generalized pseudo-Helmholtz function for stability analysis for general MASs, which fully utilizes the fact that the pseudo-Helmholtz function is a valid Lyapunov function for characterizing local asymptotic stability for complex balanced MASs. To match the use of the generalized pseudo-Helmholtz function as the Lyapunov function, the notions of reconstruction and reverse reconstruction for the original network have been defined in succession through an invertible matrix storing a positive diagonal matrix and the conserved matrix that captures mass conservation of the original network. These two kinds of defined networks will retain some properties of the original one, such as their equilibrium obtained from the concentrations of the first $n-q$ species of the equilibrium in the original network, the dynamics of the reverse reconstruction bridged to that of the original network by a constant matrix, etc. We have subsequently proved that if the original network has a complex balanced reconstruction, then it is stable. The asymptotic stability of the original network has been further reached based on the complex balanced reconstruction strategy in which the corresponding reverse reconstruction is proved to possess only an equilibrium in each positive stoichiometric compatibility class, and thus to be asymptotically stable. To facilitate applications of the proposed reconstruction strategy, a constructive approach is proposed to compute complex balanced reconstructions for general MASs.

Acknowledgments

We would like to acknowledge anonymous reviewers for their valuable comments, especially on the definition of reconstruction. Thanks are also given to Dr. Xiaoyu Zhang for her computation on examples in Table 1, and to the National Natural Science Foundation of China under Grant Nos. 11671418 and 61611130124 for financial support.