Input-to-state stability of chemical reaction networks with application to molecular computationThanks: Submitted to the editors DATE.

Renlei Jiang Email: jiangrl@zju.edu.cn Thanks: School of Mathematical Sciences, Zhejiang University, Hangzhou, China (). Xiaoyu Zhang Email: Xiaoyu_Z@seu.edu.cn Thanks: School of Mathematics, Southeast University, Nanjing, China (). Chuanhou Gao Email: gaochou@zju.edu.cn Thanks: Corresponding author. School of Mathematical Sciences, Zhejiang University, Hangzhou, China and Center for Interdisciplinary Applied Mathematics, Zhejiang University, Hangzhou, China (). Denis Dochain Email: denis.dochain@uclouvain.be Thanks: ICTEAM, UCLouvain, Bâtiment Euler, avenue Georges Lemaître 4-6, 1348 Louvain-la-Neuve, Belgium ().

Abstract

In biological reaction systems, reaction rates may vary over time due to environmental fluctuations, regulation, or coupling with other reaction modules. Input-to-state stability (ISS) provides a useful tool for analyzing the robustness of time-varying chemical reaction networks (CRNs). Existing ISS results for CRNs typically rely on restrictive structural assumptions, such as zero deficiency, a single linkage class, or weak reversibility. This paper makes two main contributions. First, we establish ISS for a broader class of weakly reversible CRNs, allowing nonzero deficiency and multiple linkage classes. Second, we extend the analysis to certain CRNs that are not weakly reversible by using network transformation techniques (linear conjugacy and reconstruction). Together, these results enlarge the class of CRNs for which robustness under time-varying reaction-rate inputs can be certified. Since CRNs are a standard framework for biomolecular computation, our results further enable the stability analysis of parallel molecular computing systems, an important problem in biomolecular computation where multiple CRN-based computing modules operate simultaneously and perturb one another through time-varying effective reaction rates.

keywords

chemical reaction network, input-to-state stability, mass-action kinetics, network transformation, parallel biomolecular computation

Funding.

This work was funded by the National Natural Science Foundation of China under Grant No. 12320101001, and the National Foreign Expert Project of China under Grant No. S20250211.

^†^†runningheads: Input-to-state stability of chemical reaction networks / R. Jiang, X. Zhang, C. Gao, D. Dochain

MSC

92C45, 93D09, 93D25

1 Introduction

Chemical reaction network (CRN) is a fundamental mathematical framework for describing chemical reactions and studying biological processes. It has been extensively studied in both deterministic and stochastic formulations [1, 21, 23]. When a CRN is endowed with mass-action kinetics, its dynamics are represented by a system of polynomial ordinary differential equations, commonly referred to as a mass-action system (MAS). The dynamical analysis of MASs has been extensively studied over the past several decades, with particular attention devoted to the roles of network structure [12, 20], equilibria [5], stability [19, 25], and persistence [4]. Among these studies, the most influential results are the Deficiency Zero Theorem [20] and complex balanced theory [25], which establish strong relationships between weak reversibility, deficiency, and the existence and stability of positive equilibria.

However, most of the above studies focused on the case in which the reaction rate constants are fixed. Relatively little attention has been paid to time-varying systems in which the reaction rate constants change over time [4, 8], despite the fact that, in practice, such constants may be affected by external factors such as temperature [4]. Input-to-state stability (ISS) provides a natural framework for analyzing the robustness of such systems with time-varying inputs. Roughly speaking, the ISS property states that, for a bounded input $u$ , the trajectory can be bounded by a function of $u$ . Furthermore, as $t$ increases, the trajectory approaches (in a Lyapunov stability manner) a ball whose radius depends on $u$ [34]. In the context of rate-controlled biochemical networks, Chaves and Sontag established semiglobal ISS results for a class of weakly reversible CRNs that have zero deficiency and a single linkage class, and applied these results to observer design [8, 7].

As the development of silicon-based computation is increasingly constrained by physical limitations, many researchers have turned their attention to biomolecular computation because of its advantages, including low energy consumption, parallel computing, and biocompatibility [9, 10]. Many of these theoretical studies are based on the mathematical framework of CRNs [2, 18, 37]. Recently, several studies have explored the possibility of operating multiple modules of biomolecular computations in parallel [3, 6], some of which make use of the ISS concept.

Despite these important results, the existing ISS theory for CRNs remains limited in several respects. First, many CRNs of interest have nonzero deficiency and multiple linkage classes. Second, a large class of biochemical and synthetic reaction networks is not weakly reversible. For these systems, classical complex balanced theory cannot be directly applied. Moreover, in applications to molecular computations, a bounded trajectory is not sufficient. It also needs to establish that a downstream system converges to the equilibrium associated with the limiting value of an upstream input signal. These issues motivate the development of ISS conditions for broader classes of CRNs and the investigation of their implications for interconnected molecular reaction modules.

In this paper, we extend the results of Chaves and Sontag to CRNs with more general structures, including networks with nonzero deficiency, multiple linkage classes, and networks that are not weakly reversible. The contributions of this paper are summarized as follows:

•

To establish the semiglobal ISS property for weakly reversible networks, we use the toric locus introduced in [12] to characterize the admissible input-value set.
•

To establish the semiglobal ISS property for some non-weakly reversible networks, we exploit the concepts of linear conjugacy [30] and reconstruction [31] to relate their dynamical properties to those of weakly reversible networks, thereby deriving the corresponding admissible input-value sets.
•

Based on the above results, we show that when the reaction-rate input converges to a constant value, the state variable also converges to an equilibrium associated with the limiting input, thereby providing a theoretical foundation for the implementation of parallel molecular computations.

The rest of this paper is organized as follows. Section 2 introduces some preliminaries on CRNs and MASs, and presents the motivation of this paper. Section 3 reviews several existing concepts and results, and establishes the semiglobal ISS property for weakly reversible networks. Section 4 establishes the semiglobal ISS property for non-weakly reversible networks. In Section 5, we discuss the application of the established results to the implementation of parallel molecular computations. Finally, Section 6 concludes this paper.

Mathematical Notation:

$\mathbb{R}^{n},\mathbb{R}_{\geq 0}^{n},\mathbb{R}_{>0}^{n},\mathbb{Z}_{\geq 0}^{n},\mathbb{Z}_{>0}^{n}$: : $n$ -dimensional real space, nonnegative real space, positive real space, nonnegative integer space, and positive integer space, respectively.
$x^{v_{\cdot j}}$: : $\prod_{i=1}^{n}x_{i}^{v_{ij}}$ , where $x\in\mathbb{R}^{n}$ , $v_{\cdot j}\in\mathbb{Z}^{n}$ and $0^{0}=1$ .
$|\cdot|$: : Euclidean norm.
ess.sup.: : essential supremum.
$\|u(t)-u^{*}\|$: : $\mathrm{ess.sup.}\{|u(t)-u^{*}|:t\geq 0\}$
$\mathcal{K},\mathcal{K}_{\infty},\mathcal{KL}$: : $\gamma:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ is a class $\mathcal{K}$ function if it is continuous, strictly increasing and satisfies $\gamma(0)=0$ ; $\gamma$ is a class $\mathcal{K}_{\infty}$ function if it is a class $\mathcal{K}$ function and satisfies $\lim_{s\to\infty}\gamma(s)=\infty$ ; $\beta:\mathbb{R}_{\geq 0}\times\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ is a class $\mathcal{KL}$ function if for each fixed $t$ the mapping $\beta(\cdot,t)$ is a class $\mathcal{K}$ function and for each fixed $s$ the function $\beta(s,t)$ decreases to zero on $t$ as $t\to\infty$ .

2 Preliminaries and motivation

In this section, we introduce some basic knowledge about CRN and MAS [21], and then present the motivation for the current study.

2.1 CRN

Consider a CRN with $n$ species, denoted by $X_{1},...,X_{n}$ , and $r$ reactions with the $j$ th reaction written as

\sum_{i=1}^{n}v_{ij}X_{i}\to\sum_{i=1}^{n}v^{\prime}_{ij}X_{i},

where $v_{ij}$ represents the stoichiometric coefficient of species $X_{i}$ in the $j$ th reaction. Furthermore, define $v_{.j}=(v_{1j},...,v_{nj})^{\top},v^{\prime}_{.j}=(v_{1j}^{\prime},...,v_{nj}^{\prime})^{\top}\in\mathbb{Z}_{\geq 0}^{n}$ as the reactant and product complexes of the $j$ th reaction, respectively. The reaction can then be written as $v_{.j}\to v^{\prime}_{.j}$ . Mathematically, a CRN is defined as follows.

Definition 1 (CRN).

A CRN consists of three finite sets:

(i)

a finite species set $\mathcal{S}=\{X_{1},...,X_{n}\}$ ;
(ii)

a finite complex set $\mathcal{C}=\bigcup_{j=1}^{r}{\left\{v_{\cdot j},v_{\cdot j}^{\prime}\right\}}$ ;
(iii)
a finite reaction set $\mathcal{R}=\bigcup_{j=1}^{r}{\left\{v_{\cdot j}\rightarrow v_{\cdot j}^{\prime}\right\}}$ satisfying
1. (a)
  
  $\forall~v_{\cdot j}\rightarrow v_{\cdot j}^{\prime}\in\mathcal{R},~v_{\cdot j}\neq v_{\cdot j}^{\prime}$ ,
2. (b)
  
  $\forall~v_{\cdot j}\in\mathcal{C},\exists~v_{\cdot j}^{\prime}\in\mathcal{C}$ such that $v_{\cdot j}\rightarrow v_{\cdot j}^{\prime}\in\mathcal{R}$ or $v_{\cdot j}^{\prime}\rightarrow v_{\cdot j}\in\mathcal{R}$ .

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

A CRN can be represented as a directed graph whose vertices correspond to complexes and whose directed edges correspond to reactions. This graph-theoretic representation allows several important structural classes of CRNs to be defined.

Definition 2 (Reversible/Weakly reversible CRN).

A CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ is

(i)

reversible if $\forall~v_{\cdot j}\to v_{\cdot j}^{\prime}\in\mathcal{R}$ , it holds $v_{\cdot j}^{\prime}\to v_{\cdot j}\in\mathcal{R}$ ;
(ii)

weakly reversible if $\forall~v_{\cdot j}\to v_{\cdot j}^{\prime}\in\mathcal{R}$ , there exists a sequence of complexes $v_{\cdot j_{1}},\cdots,v_{\cdot j_{p}}\in\mathcal{C}$ such that $v_{\cdot j}^{\prime}\to v_{\cdot j_{1}}\in\mathcal{R},v_{\cdot j_{1}}\to v_{\cdot j_{2}}\in\mathcal{R}$ , $\cdots$ , $v_{\cdot j_{p-1}}\to v_{\cdot j_{p}}\in\mathcal{R}$ , $v_{\cdot j_{p}}\to v_{\cdot j}\in\mathcal{R}$ .

Clearly, every reversible CRN is weakly reversible, but not vice versa. A CRN is weakly reversible if and only if each of its linkage class, defined as a connected component of the underlying undirected reaction graph, is strongly connected, or equivalently, if every reaction lies on a directed cycle.

For the reaction $v_{\cdot j}\rightarrow v_{\cdot j}^{\prime}$ in the CRN, we call $v_{\cdot j}^{\prime}-v_{\cdot j}$ its reaction vector, all of which induce two important concepts for dynamic analysis.

Definition 3 (Stoichiometric subspace).

For a CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ , the linear subspace $\mathscr{S}\triangleq\mathrm{span}\{v_{\cdot 1}^{\prime}-v_{\cdot 1},...,v_{\cdot r}^{\prime}-v_{\cdot r}\}$ is called the stoichiometric subspace of the network. The dimension of $\mathscr{S}$ , denoted by $\dim\mathscr{S}$ , is called the dimension of the CRN.

Definition 4 (Stoichiometric compatibility class).

For a CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ , let $x_{0}\in\mathbb{R}_{\geq 0}^{n}$ , the set $x_{0}+\mathscr{S}=\{x_{0}+x:x\in\mathscr{S}\}$ is called the stoichiometric compatibility class of $x_{0}$ . Further, $(x_{0}+\mathscr{S})\bigcap\mathbb{R}_{\geq 0}^{n}$ and $(x_{0}+\mathscr{S})\bigcap\mathbb{R}_{>0}^{n}$ are called the nonnegative stoichiometric compatibility class and the positive stoichiometric compatibility class of $x_{0}$ , respectively.

Based on these notions, we further introduce the concept of deficiency, which is often used to help characterize dynamical behaviors.

Definition 5 (Deficiency).

For a CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ , let $c$ denote the number of complexes and $\ell$ denote the number of linkage classes. The deficiency of the CRN is defined by $\delta=c-\ell-\dim\mathscr{S}$ .

The deficiency of a CRN is usually nonnegative because it can be seen as the dimension of a certain linear subspace [21]. To better understand the above abstract concepts, we illustrate them with a concrete example.

Example 6.

For the following CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ shown on the left-hand side

we get $\mathcal{S}=\{X_{1},X_{2},X_{3}\},~\mathcal{C}=\{(2,0,0)^{\top},(0,1,0)^{\top},(0,1,1)^{\top}\},~\mathcal{R}=\{(2,0,0)^{\top}\to(0,1,0)^{\top},(0,1,0)^{\top}\to(0,1,1)^{\top},(0,1,1)^{\top}\to(2,0,0)^{\top}\}$ , $\mathscr{S}=\mathrm{span}\{(-2,1,0)^{\top},(0,0,1)^{\top},\\ (2,-1,-1)^{\top}\}$ , $\dim\mathscr{S}=2$ , and $\delta=3-1-2=0$ . Clearly, it is weakly reversible, but not reversible.

2.2 MAS

When a CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ is equipped with mass-action kinetics, the rate of reaction $v_{\cdot j}\to v_{\cdot j}^{\prime}$ is measured by $\kappa_{j}x^{v_{\cdot j}}$ , where $\kappa_{j}>0$ represents the rate constant, and $x\in\mathbb{R}_{\geq 0}^{n}$ with each element $x_{i}~(i=1,...,n)$ to represent the concentration of the species $X_{i}$ .

Definition 7 (MAS).

A MAS is a CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ equipped with mass-action kinetics $\kappa=(\kappa_{1},...,\kappa_{r})^{\top}$ , often labeled by $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ or $(\mathcal{N},\kappa)$ .

The dynamics of $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ describes the change of concentrations of all species over time $t$ , and thus follows

\dot{x}=\sum_{j=1}^{r}\kappa_{j}x^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right),

(1)

which are essentially polynomial ODEs. By integrating (1) from $0$ to $t$ , we get

x(t)=x_{0}+\sum_{j=1}^{r}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)\int_{0}^{t}\kappa_{j}x^{v_{\cdot j}}(\tau)\mathrm{d}\tau,

(2)

where $x_{0}=x(0)$ is the initial state of $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ . This suggests that the state of $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ will evolve in the nonnegative stoichiometric compatibility class of $x_{0}$ , i.e., in $(x_{0}+\mathscr{S})\bigcap\mathbb{R}_{\geq 0}^{n}$ . We can also write the equivalent expression of (1) as a sum over complexes

\dot{x}=\sum_{z\in\mathcal{C}}x^{z}\sum_{\{j|v_{.j}=z\}}\kappa_{j}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right).

(3)

Revisiting Example 6 reveals that the dynamics of the MAS (on the right-hand side of the diagram in Example 6) is given by

\left(\begin{array}[]{c}\dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\end{array}\right)=\left(\begin{array}[]{c}-2\kappa_{1}x_{1}^{2}+2\kappa_{3}x_{2}x_{3}\\ \kappa_{1}x_{1}^{2}-\kappa_{3}x_{2}x_{3}\\ \kappa_{2}x_{2}-\kappa_{3}x_{2}x_{3}\end{array}\right).

Definition 8 (Equilibrium).

For a MAS $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ governed by (1), a constant vector $x^{*}\in\mathbb{R}^{n}_{>0}$ is called a positive equilibrium of the system if

\sum_{j=1}^{r}\kappa_{j}(x^{*})^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=0.

(4)

A MAS that admits a positive equilibrium is said to be balanced.

Definition 9 (Complex balanced equilibrium).

For a MAS $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ governed by (1), a positive equilibrium $x^{*}\in\mathbb{R}^{n}_{>0}$ is called a complex balanced equilibrium of the system if

\sum_{\{j|v_{\cdot j}=z\}}\kappa_{j}(x^{*})^{v_{\cdot j}}=\sum_{\{j|v_{\cdot j}^{\prime}=z\}}\kappa_{j}(x^{*})^{v_{\cdot j}},\quad\forall z\in\mathcal{C}.

(5)

A MAS that admits a complex balanced equilibrium is said to be complex balanced.

The complex balanced equilibrium is certainly an equilibrium, but not vice versa. There are some well-known results about complex balanced MAS [25, 20], including

(1)

a complex balanced MAS must be weakly reversible;
(2)

if a MAS has a complex balanced equilibrium, then all the other equilibria (if any) are complex balanced;
(3)

Deficiency Zero Theorem, which states that for any $\kappa\in\mathbb{R}^{r}_{>0}$ , the MAS $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ is complex balanced if it is weakly reversible and has zero deficiency;
(4)

within each positive stoichiometric compatibility class, there is precisely one complex balanced equilibrium, and moreover, each complex balanced equilibrium is locally asymptotically stable.

For the last result, the pseudo-Helmholtz free energy function

V(x)=\sum_{i=1}^{n}\left(x_{i}(\ln x_{i}-\ln x_{i}^{*}-1)+x^{*}_{i}\right),~x\in\mathbb{R}^{n}_{>0}

(6)

is suggested as a Lyapunov function to establish the local asymptotic stability of each complex balanced equilibrium $x^{*}$ .

2.3 Motivation

In the theoretical analysis of an MAS $(\mathcal{S},\mathcal{C},\mathcal{R},\kappa)$ , the reaction rate constants $\kappa$ are typically treated as fixed parameters [12, 20, 25, 30, 31]. However, under realistic biological conditions, these parameters often vary over time in response to external inputs, like temperature fluctuations, external electric-field effects [4], uncertainty arising from the dynamics of unmodeled species [36], or the encoding of species concentrations as the inputs for molecular computation [26, 27]. We thus consider a time-varying version of MAS governed by

\dot{x}=f(x;u)=\sum_{j=1}^{r}u_{j}(t)x^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right),

(7)

which is also called rate-controlled MAS in ISS capture [8]. Here, $x\in\mathbb{R}^{n}_{\geq 0},~u:\mathbb{R}_{\geq 0}\to\mathbb{U}$ with $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ is a piecewise locally Lipschitz function. Notably, if the system starts from an initial point $x_{0}\in\mathbb{R}^{n}_{\geq 0}$ , then the state of (7) will still evolve in $(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{\geq 0}$ .

For the purpose of capturing ISS and its application to molecular computation, there are two questions of concern about (7):

Question 1: For a bounded input $u(t)$ , can the solution $x(t)$ of (7) be bounded by an increasing function of $u$ ?

Question 2: If $\lim_{t\to\infty}u(t)=u^{*}$ , can we have $\lim_{t\to\infty}x(t)=x^{*}$ , where $x^{*}$ satisfies $f(x^{*};u^{*})=0$ ?

For Question 1, Chaves and Sontag [8, 7] have addressed it through formalizing the concept of semiglobal ISS (see definition 11 below). They also proved a class of special weakly reversible CRNs with zero deficiency and single linkage class to be semiglobal ISS. Building on this concept, we further addressed Question 2, and applied the result to enable serial computing by parallel chemical reactions [26, 27]. The strong application potential of the notion of semiglobal ISS in molecular computations motivates us to find more CRNs (not limited to single-linkage-class, deficiency-zero, weakly reversible CRNs) that can exhibit semiglobal ISS.

In what follows, we proceed by reviewing the existing result for the case of deficiency-zero, single-linkage-class, weakly reversible CRNs [8, 7], and then extending it to the cases of nonzero-deficiency weakly reversible CRNs and non-weakly reversible CRNs.

3 ISS for weakly reversible CRNs

In this section, we focus on exploring the ISS property for weakly reversible CRNs, including revisiting existing result for zero-deficiency weakly reversible CRNs, and then extending them to weakly reversible CRNs with nonzero deficiency and multiple linkage classes.

3.1 Existing result for zero-deficiency, weakly reversible CRNs

We begin by recalling several definitions [8, 7] related to ISS for MASs.

Definition 10.

For a time-varying MAS of (7) with input-value set $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ , if for each initial state $x(0)\in\mathbb{R}^{n}_{>0}$ and each input $u(\cdot)\in\mathbb{U}$ , the solution of (7) at any time $t\in J_{x(0),u}=[0,t_{max})$ with $t_{max}$ to represent the blowing-up time, is $x(t)\in\mathbb{R}^{n}_{>0}$ , then the system is $\mathbb{R}^{n}_{>0}$ -forward invariant. Further, the system is $\mathbb{R}^{n}_{>0}$ -forward complete if it is $\mathbb{R}^{n}_{>0}$ -forward invariant with $J_{x(0),u}=[0,+\infty)$ .

Definition 11.

The time-varying MAS of (7) with input-value set $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ is called semiglobally input-to-state stable with respect to $(x^{*},u^{*})$ with fixed points $x^{*}\in\mathbb{R}^{n}_{>0}$ and $u^{*}\in\mathbb{U}$ , if

(i)

the system is $\mathbb{R}^{n}_{>0}$ -forward complete;
(ii)

for every compact set $F\subseteq\mathbb{R}^{n}_{\geq 0}$ containing $x^{*}$ , and moreover, $\forall t\geq 0$ we have $x(s)\in F$ for any $s\in[0,t]$ , there exist a class $\mathcal{KL}$ function $\beta=\beta_{F}$ and a class $\mathcal{K}_{\infty}$ function $\gamma=\gamma_{F}$ such that

$|x(t)-x^{*}|\leq\beta(|x_{0}-x^{*}|,t)+\gamma(\|u-u^{*}\|)$ (8)

for each initial condition $x_{0}=x(0)\in F\cap\{x^{*}+\mathscr{S}\}\cap\mathbb{R}^{n}_{>0}$ and each input $u(\cdot)\in\mathbb{U}$ .

Definition 12.

Given a time-varying MAS of (7) with input-value set $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ , a continuous function $V:\mathbb{R}^{n}_{\geq 0}\to\mathbb{R}_{\geq 0}$ is a semiglobal ISS-Lyapunov function with respect to the fixed points $x^{*}\in\mathbb{R}^{n}_{>0}$ and $u^{*}\in\mathbb{U}$ for the system, if

(i)

the restriction of $V$ to $\mathbb{R}^{n}_{>0}$ is continuously differentiable;
(ii)

there exist class $\mathcal{K}_{\infty}$ functions $\underline{\alpha},\overline{\alpha}$ such that

$\underline{\alpha}(|x-x^{*}|)\leq V(x)\leq\overline{\alpha}(|x-x^{*}|)$

for each $x\in\mathbb{R}^{n}_{\geq 0}$ ;
(iii)

for every compact set $F\subseteq\mathbb{R}^{n}_{\geq 0}$ containing $x^{*}$ , there exist $\mathcal{K}_{\infty}$ functions $\alpha=\alpha_{F},\rho=\rho_{F}$ such that

$\nabla^{\top}V(x)f(x;u)\leq-\alpha(|x-x^{*}|)+\rho(|u-u^{*}|)$

for every $x\in F\cap\{x^{*}+\mathscr{S}\}\cap\mathbb{R}^{n}_{>0}$ and every $u\in\mathbb{U}$ .

Utilizing these concepts, Chaves and Sontag [7] proved that for a $\mathbb{R}^{n}_{>0}$ -forward complete time-varying MAS of (7) with input-value set $\mathbb{U}$ , if there is a semiglobal ISS-Lyapunov function $V$ with respect to the fixed point pair $(x^{*},u^{*})$ for the system, then it is semiglobal ISS with respect to $(x^{*},u^{*})$ . As an application, they also showed that a zero-deficiency, single-linkage-class, weakly reversible CRN is semiglobal ISS under the conditions of $f(x^{*};u^{*})=0$ and $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ being a specific set.

3.2 ISS for nonzero deficiency weakly reversible CRNs

Building on the above existing result, we try to capture semiglobal ISS for nonzero deficiency weakly reversible CRNs.

Lemma 13.

For a time-varying MAS of (7) with input-value set $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ to be a compact set, suppose that for any $u^{*}\in\mathbb{U}$ , the MAS with $u(t)\equiv u^{*}$ has a unique equilibrium $x^{*}=x^{*}(u^{*})$ in the positive stoichiometric compatibility class of $x^{*}$ . Let $V_{0}:\mathbb{R}^{n}_{\geq 0}\times\mathbb{R}^{r}_{>0}\to\mathbb{R}_{\geq 0}$ be a continuous function satisfying:

(1)

The restriction of $V_{0}$ to $\mathbb{R}^{n}_{>0}\times\mathbb{U}$ is continuously differentiable;

(2)

For any $u^{*}\in\mathbb{U}$ , there exist class $\mathcal{K}_{\infty}$ functions $\underline{\alpha}_{u^{*}},\overline{\alpha}_{u^{*}}$ such that

\underline{\alpha}_{u^{*}}(|x-x^{*}|)\leq V_{0}(x,u^{*})\leq\overline{\alpha}_{u^{*}}(|x-x^{*}|),~x\in\mathbb{R}^{n}_{\geq 0};

(3)

For any $u^{*}\in\mathbb{U}$ , it holds that

$\nabla_{x}^{\top}V_{0}(x,u^{*})f(x;u^{*})\leq 0,~x\in\mathbb{R}^{n}_{>0}$

with the equality holding if and only if $x=x^{*}$ .

Then $V(x)\triangleq V_{0}(x,u^{*})$ is a semiglobal ISS-Lyapunov function for the system (7) with respect to $(x^{*},u^{*})$ , where $x^{*}$ supports $f(x^{*};u^{*})=0$ .

Proof.

According to condition (1) and (2), we know that $V(x)$ satisfies condition (i) and (ii) in definition 12. Therefore, we just need to show that $V(x)$ satisfies condition (iii). For any compact set $F\subseteq\mathbb{R}^{n}_{\geq 0}$ containing $x^{*}$ , it holds

\begin{split}\nabla^{\top}V(x)f(x;u)&=\nabla_{x}^{\top}V_{0}(x,u^{*})f(x;u^{*})+\nabla_{x}^{\top}V_{0}(x,u^{*})f(x;u-u^{*})\\ &\triangleq I_{1}+I_{2}\end{split}

(9)

By condition (3) we can get $-I_{1}=-\nabla_{x}^{\top}V_{0}(x,u^{*})f(x;u^{*})\geq 0,~x\in F\cap(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ , with the equality holding if and only if $x=x^{*}$ . Denote $s_{\mathrm{max}}=\sup\{s:\exists~x\in F\cap(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ such that $|x-x^{*}|\geq s\}$ , and let

\tilde{\alpha}_{F}(s)=\inf_{\begin{subarray}{c}|x-x^{*}|\geq s\\ x\in F\cap(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}\end{subarray}}\left(-\nabla_{x}^{\top}V_{0}(x,u^{*})f(x;u^{*})\right),\quad 0\leq s\leq s_{\mathrm{max}}.

Since $\tilde{\alpha}_{F}(0)=0$ and $\tilde{\alpha}_{F}(s)>0$ for all $s\in(0,s_{\max}]$ , and $\tilde{\alpha}_{F}$ is continuous and nondecreasing on $[0,s_{\max}]$ , there exists a class $\mathcal{K}$ function $\alpha_{F}$ such that $\alpha_{F}(s)\leq\tilde{\alpha}_{F}(s)/2$ for all $s\in[0,s_{\max}]$ . Moreover, $\alpha_{F}$ can be extended to $\mathbb{R}_{\geq 0}$ as a class $\mathcal{K}_{\infty}$ function, for instance by defining $\alpha_{F}(s)=\alpha_{F}(s_{\max})+s-s_{\max}$ for $s>s_{\max}$ . Thus the following equation holds

-I_{1}\geq\tilde{\alpha}_{F}(|x-x^{*}|)\geq\alpha_{F}(|x-x^{*}|),~\forall~x\in F\cap(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}.

(10)

Since $F$ is a compact set and $\nabla_{x}^{\top}V_{0}(x,u^{*})$ is continuous, there is a constant $C_{F}>0$ such that

\begin{split}I_{2}&=\nabla_{x}^{\top}V_{0}(x,u^{*})\sum_{j=1}^{r}(u_{j}-u_{j}^{*})x^{v_{\cdot j}}(v_{\cdot j}^{\prime}-v_{\cdot j})\\ &\leq C_{F}\|u-u^{*}\|=\rho_{F}(\|u-u^{*}\|)\end{split}

(11)

where $x\in F\cap(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ , $\rho_{F}(s)=C_{F}s$ is a class $\mathcal{K}_{\infty}$ function. Substituting (10) and (11) into (9) we know that the function $V(x)$ satisfies condition (iii) in definition 12, and the proof is complete.

Proposition 14.

Consider a time-varying MAS of (7) with input-value set $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ to be a compact set, and suppose that it is $\mathbb{R}^{n}_{>0}$ -forward complete. Then for any $u^{*}\in\mathbb{U}$ , under the conditions given in lemma 13 the system (7) starting from $x_{0}$ is semiglobal ISS with respect to $(x^{*},u^{*})$ , where $x^{*}$ supports $f(x^{*};u^{*})=0$ and is contained in a compact set $F\subseteq\mathbb{R}^{n}_{\geq 0}$ , and moreover, $x_{0}\in F\cap\{x^{*}+\mathscr{S}\}\cap\mathbb{R}^{n}_{>0}$ .

Proof.

The proof follows directly by applying Lemma 3.7 in [7] and lemma 13.

proposition 14 provides some conditions to support the MAS of (7) to be semiglobal ISS. In particular, if the network is weakly reversible, we can get the main result in this section. For this purpose, we introduce an important concept that characterizes the set $\mathbb{U}$ .

Definition 15 (Toric Locus [12]).

For a CRN $\mathcal{N}=(\mathcal{S},\mathcal{C},\mathcal{R})$ , its toric locus $\mathcal{T}(\mathcal{N})\subseteq\mathbb{R}^{r}_{>0}$ denotes the set of parameters $\kappa\in\mathbb{R}^{r}_{>0}$ , for which the corresponding MAS $(\mathcal{N},\kappa)$ is complex balanced.

As discussed in [12], for a given CRN $\mathcal{N}$ , its toric locus $\mathcal{T}(\mathcal{N})\neq\varnothing$ if and only if $\mathcal{N}$ is weakly reversible.

Theorem 16.

Consider a time-varying MAS of (7) with weakly reversible network structure. Denote its toric locus by $\mathcal{T}$ , then for any compact input-value set $\mathbb{U}\subseteq\mathcal{T}$ and any $u^{*}\in\mathbb{U}$ , the system (7) starting from $x_{0}$ is semiglobal ISS with respect to $(x^{*},u^{*})$ , where $x^{*}$ supports $f(x^{*};u^{*})=0$ and is contained in a compact set $F\subseteq\mathbb{R}^{n}_{\geq 0}$ , and moreover, $x_{0}\in F\cap\{x^{*}+\mathscr{S}\}\cap\mathbb{R}^{n}_{>0}$ .

Proof.

We complete the proof in three steps.

Step 1. We prove that system (7), with input-value set $\mathbb{U}$ , is $\mathbb{R}^{n}_{>0}$ -forward invariant, i.e., for each initial state $x(0)=x_{0}\in\mathbb{R}^{n}_{>0}$ and each $\mathbb{U}$ -valued input $u(\cdot)$ , the corresponding maximal solution of (7), which is defined on $[0,t_{max})$ , satisfies $x(t)\in(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{>0},~0\leq t<t_{max}$ . This proof is very similar to the proof of Proposition 5.1 in [8], so we will not reproduce it here.

Step 2. We prove that system (7), with input-value set $\mathbb{U}$ , is $\mathbb{R}^{n}_{>0}$ -forward complete, i.e., $t_{max}=+\infty$ . To prove this claim, we argue by contradiction. Suppose $t_{max}<+\infty$ , then it holds $\lim_{t\to t_{max}}|x(t)|=+\infty$ . Since $\mathbb{U}\subseteq\mathcal{T}$ , for any $u\in\mathbb{U}$ , the MAS is complex balanced, and thus there is precisely one equilibrium in each positive stoichiometric compatibility class. Now consider the function

W(t)=V(x(t),x^{*}(u(t)))=\sum_{i=1}^{n}\left(x_{i}(t)(\ln x_{i}(t)-\ln x_{i}^{*}(u(t))-1)+x^{*}_{i}(u(t))\right),

where $x^{*}(u(t))$ is the unique equilibrium of $\dot{x}=f(x,u(t))$ in $(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ . Clearly, it holds $\lim_{t\to t_{max}}W(t)=+\infty$ , and for $t\in[0,t_{\max})$ we have

\begin{split}\frac{dW}{dt}&=\nabla_{x}^{\top}V(x(t),x^{*}(u(t)))\dot{x}+\nabla^{\top}_{x^{*}}V(x(t),x^{*}(u(t)))\dot{x}^{*}\\ &=\left(\mathrm{Ln}(x)-\mathrm{Ln}(x^{*})\right)^{\top}f(x,u)+\sum_{i=1}^{n}\frac{1}{x_{i}^{*}}(x_{i}^{*}-x_{i})\dot{x}_{i}^{*}\end{split}

(12)

where $\mathrm{Ln}(x)=(\ln x_{1},...,\ln x_{n})^{\top},~\dot{x}^{*}=\frac{d}{dt}(x^{*}(u(t)))$ . Since (6) is the Lyapunov function of complex balanced MAS, we have

\left(\mathrm{Ln}(x)-\mathrm{Ln}(x^{*})\right)^{\top}f(x,u)\leq 0.

(13)

In addition, according to the Corollary 3.13 in [12], the function $x^{*}=x^{*}(u)$ is continuously differentiable. As $u(t)\in\mathbb{U}$ , where $u(t)$ is piecewise locally Lipschitz function, and $\mathbb{U}$ is a compact set, we have

\dot{x}_{i}^{*}=\sum_{j=1}^{r}\frac{dx_{i}^{*}}{du_{j}}\frac{du_{j}}{dt}\leq\sum_{j=1}^{r}\left|\frac{dx_{i}^{*}}{du_{j}}\right|\cdot\left|\frac{du_{j}}{dt}\right|\leq\sum_{j=1}^{r}C_{1}\cdot L=C_{1}Lr,

(14)

and

1/x_{i}^{*}\leq C_{2},\quad x_{i}^{*}\leq C_{2}

(15)

where $C_{1},C_{2},L>0$ are constants. Finally, to estimate $x_{i}^{*}-x_{i}$ , we define the function

\varphi(r)=r(\ln r-\ln x_{i}^{*}-1)+2x_{i}^{*}-|x_{i}^{*}-r|.

For $r\geq x_{i}^{*}$ , we have $\varphi^{\prime}(r)=\ln r-\ln x_{i}^{*}-1$ , and $\varphi^{\prime}(r)>0\Longleftrightarrow r>ex_{i}^{*}$ , therefore it holds $\varphi(r)\geq\varphi(ex_{i}^{*})=(3-e)x_{i}^{*}>0$ ; For $0<r<x_{i}^{*}$ , we have $\varphi^{\prime}(r)=\ln r-\ln x_{i}^{*}+1$ , and $\varphi^{\prime}(r)>0\Longleftrightarrow x_{i}^{*}/e<r<x_{i}^{*}$ , therefore it holds $\varphi(r)\geq\varphi(x_{i}^{*}/e)=(1-1/e)x_{i}^{*}>0$ . In a word, we have $\varphi(r)>0,~\forall r>0$ . Specially, $\varphi(x_{i})>0$ , which means

|x_{i}^{*}-x_{i}|<x_{i}(\ln x_{i}-\ln x_{i}^{*}-1)+2x_{i}^{*}

(16)

Substituting eqs. 13, 14, 15, and 16 into (12) we obtain

\begin{split}\frac{dW}{dt}&\leq C_{1}C_{2}Lr\sum_{j=1}^{n}\left(x_{i}(\ln x_{i}-\ln x_{i}^{*}-1)+2x_{i}^{*}\right)\\ &\leq C_{1}C_{2}LrW+C_{1}C_{2}^{2}Lrn=C_{3}W+C_{4},\end{split}

(17)

where $C_{3}=C_{1}C_{2}Lr,~C_{4}=2C_{1}C_{2}^{2}Lrn>0$ are constants. Applying Gronwall’s inequality‌ to (17) it yields

W(t)\leq\left(W_{0}+\frac{C_{4}}{C_{3}}\right)e^{C_{3}t}-\frac{C_{4}}{C_{3}},\quad t\in[0,t_{max}),

(18)

which is a contradiction to $\lim_{t\to t_{max}}W(t)=+\infty$ . Therefore, we have $t_{max}=+\infty$ .

Step 3. We prove the result of this theorem. Since for any $u^{*}\in\mathbb{U}$ , the MAS is complex balanced, and thus the function (6) satisfies all the conditions in lemma 13. By proposition 14, the conclusion is true.

Remark 17.

theorem 16 imposes no restriction on the deficiency or on the number of linkage classes. It can therefore be viewed as a generalization of the result in [8].

Example 18.

Consider the p53 signaling and dimerization network [32, 33]

where $X_{1},X_{2},X_{3}$ represent p53 tumor suppressor, its negative regulator Mdm2, and its dimerization product, respectively. The first four reactions (left hand side) represent the negative feedback regulation process of the p53–Mdm2 system, while the last two reactions (right hand side) represent the dimerization process of p53. The dynamics of this time-varying MAS is

\left(\begin{array}[]{c}\dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\end{array}\right)=\left(\begin{array}[]{c}u_{1}-u_{3}x_{1}x_{2}-2u_{5}x_{1}^{2}+2u_{6}x_{3}\\ u_{2}x_{1}-u_{4}x_{2}\\ u_{5}x_{1}^{2}-u_{6}x_{3}\end{array}\right).

(19)

Clearly, the CRN is weakly reversible, and has 6 complexes, 2 linkage classes and $\dim\mathscr{S}=3$ , which means its deficiency $\delta=6-2-3=1$ . Moreover, the MAS is complex balanced if and only if there exists a $x^{*}\in\mathbb{R}^{3}_{>0}$ such that

\begin{cases}u_{2}x_{1}^{*}=u_{1}\\ u_{3}x_{1}^{*}x_{2}^{*}=u_{2}x_{1}^{*}\\ u_{4}x_{2}^{*}=u_{3}x_{1}^{*}x_{2}^{*}\\ u_{5}(x_{1}^{*})^{2}=u_{6}x_{3}^{*}\end{cases}\Longleftrightarrow~u_{1}u_{3}=u_{2}u_{4},~u_{5},u_{6}>0,

i.e., the toric locus $\mathcal{T}=\{u\in\mathbb{R}^{6}_{>0}:u_{1}u_{3}=u_{2}u_{4},~u_{5},u_{6}>0\}$ . By setting $u^{*}=(2,1,2,4,1,4)^{\top}\in\mathcal{T}$ , we get $x^{*}=(2,0.5,1)^{\top}$ . Then according to theorem 16, for any compact set $\mathbb{U}\subseteq\mathcal{T}$ containing $u^{*}$ , this MAS, with input-value set $\mathbb{U}$ , is semiglobal ISS with respect to $(x^{*},u^{*})$ . To observe the ISS behavior more intuitively, we further make numerical simulations for this system with two kinds of different inputs starting from the same initial point $x_{0}=(3.7,1.4,0.2)^{\top}$ , shown in fig. 1. As can be seen, for the first kind of input (converging to $u^{*}$ ), the state will converge to $x^{*}$ , while for the second kind of input (oscillating along $u^{*}$ ), the long-term behavior of the state is controlled by a function of $|u-u^{*}|$ .

Refer to caption — Figure 1: ISS exhibition of (19) with different inputs: (I) $u_{1}(t)=2-e^{-t},u_{2}(t)=1,u_{3}(t)=2,u_{4}(t)=4-2e^{-t},u_{5}(t)=1+\frac{\sin(2t)}{2+7t},u_{6}(t)=4-\frac{1}{1+3t}$ , which converge to $u^{*}$ ; (II) $u_{1}(t)=2-0.5\cos t,u_{2}(t)=1,u_{3}(t)=2,u_{4}(t)=4-\cos t,u_{5}(t)=1+0.3\sin(2t),u_{6}(t)=4-\frac{1}{1+3t}$ , which oscillate along $u^{*}$ .

4 ISS for non-weakly reversible CRNs

In this section, we further extend the analysis of ISS to non-weakly reversible CRNs.

In general, it is difficult to capture the ISS property of a MAS only based on the dynamics. The above analysis indicates that the structure information of network is quite crucial in catching ISS. To fully utilize this point, we try to make a bridge between the dynamics of the concerned MAS and that of weakly reversible MAS by leveraging the notions of linear conjugacy [30] and reconstruction [31].

4.1 Dynamic bridge through linear conjugacy

We first give the definition of linear conjugacy.

Definition 19 (Linear conjugacy [30]).

Two MASs $(\mathcal{N},u)$ and $(\tilde{\mathcal{N}},\tilde{u})$ are linearly conjugate if there exists a linear, bijective mapping $h:\mathbb{R}^{n}_{>0}\to\mathbb{R}^{n}_{>0}$ such that $h\left(\Phi(x_{0},t)\right)=\tilde{\Phi}\left(\tilde{x}_{0},t\right)$ for all $x_{0}\in\mathbb{R}^{n}_{>0}$ and $\tilde{x}_{0}=h(x_{0})$ , where $\Phi(x_{0},t)$ and $\tilde{\Phi}(\tilde{x}_{0},t)$ denote the system orbits of the two MASs with initial values $x_{0}$ and $\tilde{x}_{0}$ , respectively.

As discussed in [30], linear conjugacy also means that there exists a positive, diagonal matrix $D=\mathrm{diag}(d_{1},...,d_{n})$ to link the dynamics of two MASs through $D\dot{x}=\dot{\tilde{x}}$ , i.e.,

D\sum_{j=1}^{r}u_{j}x^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=\sum_{j=1}^{\tilde{r}}\left(\tilde{u}_{j}\prod_{i=1}^{n}d_{i}^{\tilde{v}_{ij}}\right)x^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right),\quad\forall~x\in\mathbb{R}^{n}_{\geq 0}.

(20)

Specially, when $D$ is the identity matrix, $(\mathcal{N},u)$ and $(\tilde{\mathcal{N}},\tilde{u})$ are said to be dynamically equivalent [13, 17].

Utilizing the equivalent expression of dynamics (see (3)) and denoting the set of reactant complexes by $\mathcal{C}_{react}=\{v_{\cdot 1},...,v_{\cdot r}\}\cup\{\tilde{v}_{\cdot 1},...,\tilde{v}_{\cdot\tilde{r}}\}$ , we can rewrite (20) as

D\sum_{z\in\mathcal{C}_{react}}x^{z}\sum_{\{j|v_{\cdot j}=z\}}u_{j}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=\sum_{z\in\mathcal{C}_{react}}x^{z}\sum_{\{j|\tilde{v}_{\cdot j}=z\}}\left(\tilde{u}_{j}\prod_{i=1}^{n}d_{i}^{z_{i}}\right)\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right).

Note that $z$ is any complex in $\mathcal{C}_{react}$ , so the above equation is equivalent to

D\sum_{\{j|v_{\cdot j}=z\}}u_{j}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=\sum_{\{j|\tilde{v}_{\cdot j}=z\}}\left(\tilde{u}_{j}\prod_{i=1}^{n}d_{i}^{z_{i}}\right)\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right),\quad\forall~z\in\mathcal{C}_{react}.

(21)

Then we present the main theorem of this subsection, which characterizes the input-value set $\mathbb{U}$ using the concept of linear conjugacy.

Theorem 20.

For a time-varying MAS $(\mathcal{S},\mathcal{C},\mathcal{R},u)$ of (7), suppose there is another time-varying MAS $(\mathcal{S},\tilde{\mathcal{C}},\tilde{\mathcal{R}},\tilde{u})$ , which is weakly reversible and has toric locus $\mathcal{T}(\tilde{\mathcal{N}})$ with $\tilde{\mathcal{N}}=(\mathcal{S},\tilde{\mathcal{C}},\tilde{\mathcal{R}})$ , to be linear conjugate to it, i.e., there exists a positive, diagonal matrix $D=\mathrm{diag}(d_{1},...,d_{n})$ $(d_{i}>0)$ to support (21). Further, suppose that $\tilde{u}$ depends smoothly on $u\in\mathcal{U}$ with

\mathcal{U}=\{u\in\mathbb{R}^{r}_{>0}:~\exists~\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N}})~such~that~\forall~z\in\mathcal{C}_{react}~it~holds~(\ref{eq_def_LC})\}.

(22)

Then for any compact set $\mathbb{U}\subseteq\mathcal{U}$ and any $u^{*}\in\mathbb{U}$ , the system (7), with input-value set $\mathbb{U}$ , is semiglobal ISS with respect to $(x^{*},u^{*})$ , where $x^{*}$ satisfies $f(x^{*};u^{*})=0$ .

To prove theorem 20, we introduce two lemmas with the proofs given in the Appendix.

Lemma 21.

Suppose all the conditions in theorem 20 are satisfied. Then for any $u\in\mathcal{U}$ given in (22), the MAS $(\mathcal{S},\mathcal{C},\mathcal{R},u)$ has a unique equilibrium $x^{*}=x^{*}(u)$ in each positive stoichiometric compatibility class. Moreover, $x^{*}=x^{*}(u)$ depends smoothly on the parameter values $u\in\mathcal{U}$ .

Lemma 22.

Suppose all the conditions in theorem 20 are satisfied. Then for any compact set $\mathbb{U}\subseteq\mathcal{U}$ given in (22), there exists a function $V_{0}(x,u)$ satisfying condition (1)-(3) in lemma 13.

Proof of theorem 20.

According to lemma 21 and lemma 22, we know that all the conditions in lemma 13 are satisfied. Then by proposition 14, the conclusion is true.

Example 23.

We use the MAS $(\mathcal{N},u)$ discussed in [30] to exhibit ISS, which takes

X_{1}+2X_{2}\overset{u_{1}}{\longrightarrow}X_{1}+3X_{2}\overset{u_{2}}{\longrightarrow}X_{1}+X_{2}\overset{u_{3}}{\longrightarrow}3X_{1},~2X_{1}\overset{u_{4}}{\longrightarrow}X_{2}

with dynamics

\left(\begin{array}[]{c}\dot{x}_{1}\\ \dot{x}_{2}\end{array}\right)=\left(\begin{array}[]{c}2u_{3}x_{1}x_{2}-2u_{4}x_{1}^{2}\\ u_{1}x_{1}x_{2}^{2}-2u_{2}x_{1}x_{2}^{3}-u_{3}x_{1}x_{2}+u_{4}x_{1}^{2}\end{array}\right).

(23)

This MAS was said to be linear conjugate to the following MAS $(\tilde{\mathcal{N}},\tilde{u})$

X_{1}+2X_{2}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{1}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{2}}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{1}}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{2}}}\mkern 9.0mu$}}}}$}}{}X_{1}+3X_{2},~X_{1}+X_{2}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{3}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{4}}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{3}}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{4}}}\mkern 9.0mu$}}}}$}}{}2X_{1}

with dynamics

\left(\begin{array}[]{c}\dot{x}_{1}\\ \dot{x}_{2}\end{array}\right)=\left(\begin{array}[]{c}\tilde{u}_{3}x_{1}x_{2}-\tilde{u}_{4}x_{1}^{2}\\ \tilde{u}_{1}x_{1}x_{2}^{2}-\tilde{u}_{2}x_{1}x_{2}^{3}-\tilde{u}_{3}x_{1}x_{2}+\tilde{u}_{4}x_{1}^{2}\end{array}\right),

(24)

linked by $D=\mathrm{diag}(1,2)$ . Note that the second MAS is weakly reversible, and has deficiency $\delta=4-2-2=0$ , so we have $\mathcal{T}(\tilde{\mathcal{N}})=\mathbb{R}^{4}_{>0}$ . Then (21) follows

\begin{cases}\left(\begin{array}[]{cc}1&\\ &2\\ \end{array}\right)u_{1}\left(\begin{array}[]{c}0\\ 1\\ \end{array}\right)=\tilde{u}_{1}\cdot 4\left(\begin{array}[]{c}0\\ 1\\ \end{array}\right),\\ \left(\begin{array}[]{cc}1&\\ &2\\ \end{array}\right)u_{2}\left(\begin{array}[]{c}0\\ -2\\ \end{array}\right)=\tilde{u}_{2}\cdot 8\left(\begin{array}[]{c}0\\ -1\\ \end{array}\right),\\ \left(\begin{array}[]{cc}1&\\ &2\\ \end{array}\right)u_{3}\left(\begin{array}[]{c}2\\ -1\\ \end{array}\right)=\tilde{u}_{3}\cdot 2\left(\begin{array}[]{c}1\\ -1\\ \end{array}\right),\\ \left(\begin{array}[]{cc}1&\\ &2\\ \end{array}\right)u_{4}\left(\begin{array}[]{c}-2\\ 1\\ \end{array}\right)=\tilde{u}_{4}\left(\begin{array}[]{c}-1\\ 1\\ \end{array}\right).\end{cases}\Longleftrightarrow~\begin{cases}u_{1}=2\tilde{u}_{1},\\ u_{2}=2\tilde{u}_{2},\\ u_{3}=\tilde{u}_{3},\\ u_{4}=\frac{1}{2}\tilde{u}_{4}.\end{cases}

(25)

Therefore, according to (22) we have

\mathcal{U}=\left\{u\in\mathbb{R}^{4}_{>0}:~(u_{1},u_{2},u_{3},u_{4})=\left(2\tilde{u}_{1},2\tilde{u}_{2},\tilde{u}_{3},\frac{1}{2}\tilde{u}_{4}\right),\tilde{u}\in\mathbb{R}^{4}_{>0}\right\}=\mathbb{R}^{4}_{>0}.

By setting $u^{*}=(2,1,3,1)^{\top}\in\mathcal{U}$ , we get $x^{*}=(1,3)^{\top}$ . Then according to theorem 20, for any compact set $\mathbb{U}\subseteq\mathbb{R}^{4}_{>0}$ containing $u^{*}$ , the MAS (23), with input-value set $\mathbb{U}$ , is semiglobal ISS with respect to $(x^{*},u^{*})$ .

The network in Example 23 is an artificially constructed abstract network and does not necessarily have a direct biological interpretation. In Appendix B, we establish the ISS of a biological network with practical significance.

As shown in Example 23, once the complex balanced MAS $(\tilde{\mathcal{N}},\tilde{u})$ and the linking matrix $D$ are specified, determining the set $\mathcal{U}$ reduces to solving a system of linear equations (like (25)). However, this system of linear equations may have no solution. That is, the set $\mathcal{U}$ in theorem 20 may be empty. To avoid this possibility, it is necessary to investigate the conditions under which $\mathcal{U}$ is nonempty. For this purpose, similar to the methods in [16] and [30], we define two sets

C_{\mathcal{R},\mathbb{K}}(z)=\left\{\sum_{\{j|v_{\cdot j}=z\}}u_{j}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right):u\in\mathbb{K}\right\}

(26)

and

C_{\tilde{\mathcal{R}},\tilde{\mathbb{K}}}(z)=\left\{\sum_{\{j|\tilde{v}_{\cdot j}=z\}}\left(\tilde{u}_{j}\prod_{i=1}^{n}d_{i}^{z_{i}}\right)\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right):\tilde{u}\in\tilde{\mathbb{K}}\right\},

(27)

where $z\in\mathcal{C}_{react},~\mathbb{K}\subseteq\mathbb{R}^{r}_{>0}$ and $\tilde{\mathbb{K}}\subseteq\mathbb{R}^{\tilde{r}}_{>0}$ .

Proposition 24.

The set $\mathcal{U}$ defined by (22) is not empty if and only if for every $z\in\mathcal{C}_{react}$ we have $D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)\cap C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)\neq\varnothing$ , where $D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)=\{D\xi:\xi\in C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)\}$ .

Proof.

The proof is trivial, since $D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)\cap C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)\neq\varnothing,~\forall~z\in\mathcal{C}_{react}$ if and only if $\exists~u\in\mathbb{R}^{r}_{>0},~\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N}})$ such that (21) holds.

Another point worth noting is the condition under which $\mathcal{U}=\mathbb{R}^{r}_{>0}$ , which implies that $\mathbb{U}$ can be any compact subset of positive orthant $\mathbb{R}^{r}_{>0}$ .

Proposition 25.

The set defined by (22) satisfies $\mathcal{U}=\mathbb{R}^{r}_{>0}$ if and only if for every $z\in\mathcal{C}_{react}$ we have $D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)\subseteq C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)$ .

Proof.

“if”: $\forall~u\in\mathbb{R}^{r}_{>0}$ , by the given condition we know that $\exists~\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N}})$ such that for every $z\in\mathcal{C}_{react}$ , it holds (21), which means $u\in\mathcal{U}$ . Since $u$ is arbitrary, it follows that $\mathbb{R}^{r}_{>0}\subseteq\mathcal{U}$ . It is obvious that $\mathcal{U}\subseteq\mathbb{R}^{r}_{>0}$ , we thus have $\mathcal{U}=\mathbb{R}^{r}_{>0}$ .

“only if”: For every $z\in\mathcal{C}_{react}$ , $\forall~\xi\in D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)$ , $\exists~u\in\mathbb{R}^{r}_{>0}$ such that $\xi=D\sum_{\{j|v_{\cdot j}=z\}}u_{j}(v_{\cdot j}^{\prime}-v_{\cdot j})$ . Since $\mathcal{U}=\mathbb{R}^{r}_{>0}$ , $\exists~\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N}})$ such that

\xi=D\sum_{\{j|v_{\cdot j}=z\}}u_{j}(v_{\cdot j}^{\prime}-v_{\cdot j})=\sum_{\{j|\tilde{v}_{\cdot j}=z\}}\left(\tilde{u}_{j}\prod_{i=1}^{n}d_{i}^{z_{i}}\right)\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right),

which means $\xi\in C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)$ . Since $\xi$ is arbitrary, it follows that $D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)\subseteq C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)$ .

Revisit Example 23, for $z=(1,2)^{\top}\in\mathcal{C}_{react}$ , we can get

\begin{split}D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)&=\left\{u_{1}\left(\begin{array}[]{c}0\\ 2\\ \end{array}\right):u\in\mathbb{R}^{4}_{>0}\right\}\\ C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)&=\left\{\tilde{u}_{1}\left(\begin{array}[]{c}0\\ 4\\ \end{array}\right):\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N})}=\mathbb{R}^{4}_{>0}\right\}\end{split}

Proposition 26.

The set defined by (22) satisfies $\mathcal{U}=\mathbb{R}^{r}_{>0}$ if the following conditions hold:

(1)

$\tilde{\mathcal{N}}$ is weakly reversible and has zero deficiency;
(2)

$r=\tilde{r}$ , and $v_{\cdot i}\neq v_{\cdot j},~\forall~1\leq i<j\leq r$ ;
(3)

$\exists~u,\tilde{u}\in\mathbb{R}^{r}_{>0}$ such that (21) holds.

Proof.

We complete the proof in three steps.

Step 1. Denote $\mathcal{C}_{1}=\{v_{\cdot 1},...,v_{\cdot r}\},~\tilde{\mathcal{C}}_{1}=\{\tilde{v}_{\cdot 1},...,\tilde{v}_{\cdot j}\}$ , then we prove that $\mathcal{C}_{1}\subseteq\tilde{\mathcal{C}}_{1}$ . This proof is proceeded by contradiction. Suppose $\exists~1\leq j\leq r$ such that $v_{\cdot j}\notin\tilde{\mathcal{C}}_{1}$ . According to condition (3) and $v_{\cdot i}\neq v_{\cdot j},~\forall~i\neq j$ , we have

Du_{j}(v_{\cdot j}^{\prime}-v_{\cdot j})=\sum_{\{l|\tilde{v}_{\cdot l}=v_{\cdot j}\}}\tilde{u}_{l}\left(\prod_{i=1}^{n}d_{i}^{v_{ij}}\right)(\tilde{v}_{\cdot l}^{\prime}-\tilde{v}_{\cdot l})=0.

Since $D$ is invertible and $v_{\cdot j}\neq v_{\cdot j}^{\prime}$ , we can get $u_{j}=0$ . This is in contradiction to $u\in\mathbb{R}^{r}_{>0}$ . Therefore, we have $v_{\cdot j}\in\tilde{\mathcal{C}}_{1},~\forall~1\leq j\leq r$ , i.e., $\mathcal{C}_{1}\subseteq\tilde{\mathcal{C}}_{1}$ .

Step 2. We prove that $\mathcal{C}_{1}=\tilde{\mathcal{C}}_{1}$ , and $\tilde{v}_{\cdot i}\neq\tilde{v}_{\cdot j},~\forall~1\leq i<j\leq\tilde{r}$ . The proof is trivial, since $\mathcal{C}_{1}\subseteq\tilde{\mathcal{C}}_{1}$ and $r=\tilde{r}$ .

Step 3. We prove that $\mathcal{U}=\mathbb{R}^{r}_{>0}$ . According to Step 2 we have $\mathcal{C}_{react}=\mathcal{C}_{1}=\tilde{\mathcal{C}}_{1}$ . For simplicity and without loss of generality, we set $v_{\cdot j}=\tilde{v}_{\cdot j},~1\leq j\leq r$ . For any $\mu\in\mathbb{R}^{r}_{>0}$ , denote $c=(\frac{\mu_{1}}{u_{1}},...,\frac{\mu_{r}}{u_{r}})^{\top}\in\mathbb{R}^{r}_{>0}$ , then $\forall~1\leq j\leq r$ we have

\begin{split}D\mu_{j}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)&=c_{j}Du_{j}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)\\ &=c_{j}\tilde{u}_{j}\left(\prod_{i=1}^{n}d_{i}^{v_{ij}}\right)(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j})\\ &=\tilde{\mu}_{j}\left(\prod_{i=1}^{n}d_{i}^{v_{ij}}\right)(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j})\end{split}

where $\tilde{\mu}=(c_{1}\tilde{u}_{1},...,c_{r}\tilde{u}_{r})^{\top}\in\mathbb{R}^{r}_{>0}$ . By condition (1) it follows $\mathbb{R}^{r}_{>0}=\mathcal{T}(\tilde{\mathcal{N}})$ . Therefore, we have $\mu\in\mathcal{U}$ . Since $\mu$ is arbitrary, it follows that $\mathbb{R}^{r}_{>0}\subseteq\mathcal{U}$ . It is obvious that $\mathcal{U}\subseteq\mathbb{R}^{r}_{>0}$ , we thus have $\mathcal{U}=\mathbb{R}^{r}_{>0}$ .

Recall that $\tilde{\mathcal{N}}$ in Example 23 is weakly reversible and has zero deficiency, and $r=\tilde{r}=4,~v_{\cdot i}\neq v_{\cdot j},~1\leq i<j\leq 4$ . In addition, let $u=(2,2,1,1)^{\top},~\tilde{u}=(1,1,1,2)^{\top}$ , and then (21) holds. Then by proposition 26 we have $\mathcal{U}=\mathbb{R}^{4}_{>0}$ .

It should be noted that all the above propositions are derived under the condition that $D$ and $\tilde{\mathcal{N}}$ are given. When only $\mathcal{N}$ is provided, finding a weakly reversible network $\tilde{\mathcal{N}}$ that is linearly conjugate to it under $D$ is not an easy task. This task is called “linearly conjugate weakly reversible realizations”, and in the case $D=I$ , it is also known as “weakly reversible realizations”. The main approach to solve these problems is to reformulate it as a mixed-integer linear programming problem [35, 28, 29, 14, 11].

4.2 Dynamic bridge through reconstruction

In this subsection, we utilize another notion of reconstruction [31] to capture ISS for a broader class of MASs but assisted with conservation laws. We first present the related concepts.

Definition 27 (Conservative CRN).

A CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ is conservative if there exists a vector $\rho\in\mathscr{S}^{\bot}\cap\mathbb{R}^{n}_{>0}$ , where $\mathscr{S}^{\bot}\triangleq\{\xi\in\mathbb{R}^{n}:\xi^{\top}s=0,~\forall s\in\mathscr{S}\}$ , and $\rho$ is called a conserved vector.

Definition 28 (Conserved matrix).

For a conservative CRN $(\mathcal{S},\mathcal{C},\mathcal{R})$ , assume that there exist $q$ linearly independent conserved vectors $\xi_{j}\in\mathscr{S}^{\bot},~j=1,...,q$ , where $1\leq q\leq\dim\left(\mathscr{S}^{\bot}\right)$ , then the matrix $C\in\mathbb{R}^{n\times q}_{>0}$ equipped with $\xi_{j}$ as column vectors is called a conserved matrix for the network.

Note that the conserved matrix $C$ for a CRN is not unique, but all of them are full column rank. For a CRN that is not conservative, we have $q=0$ , i.e., the conserved matrix does not exist. The following lemma characterizes the number of nonfree variables in a conservative CRN.

Lemma 29 ([31], Theorem 2).

For a conservative MAS $(\mathcal{S},\mathcal{C},\mathcal{R},u)$ , suppose $C\in\mathbb{R}^{n\times q}_{>0}$ is a conserved matrix and $x^{*}$ is an equilibrium. Then in $(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{\geq 0}$ there are at least $q$ nonfree state variables, and moreover, these variables are uniquely determined by other $n-q$ state variables and $x^{*}$ , but independent of $C$ .

For simplicity and without loss of generality, we assume throughout the remainder of the paper that for a conservative MAS $(\mathcal{S},\mathcal{C},\mathcal{R},u)$ , the last $q$ state variables, from $x_{n-q+1}$ to $x_{n}$ , are nonfree while others (from $x_{1}$ to $x_{n-q}$ ) are free. We then give the definition of reconstruction, which will play a crucial role in the subsequent analysis.

Definition 30 (Reconstruction [31]).

For a conservative MAS $(\mathcal{N},u)$ of (7), let $x^{*}\in\mathbb{R}^{n}_{>0}$ be an equilibrium and $C\in\mathbb{R}^{n\times q}_{>0}$ be a conserved matrix. If there exists a positive matrix $D_{1}=\mathrm{diag}(d_{1},...,d_{n-q})$ and a MAS $(\tilde{\mathcal{N}},\tilde{u})$ given by $\dot{\tilde{x}}=\sum_{j=1}^{\tilde{r}}\tilde{u}_{j}\tilde{x}^{\tilde{v}_{\cdot j}}(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j})$ such that $\forall~x\in(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{\geq 0}$ it holds

D\sum_{j=1}^{r}u_{j}x^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=\left(\begin{array}[]{c}\sum_{j=1}^{\tilde{r}}\tilde{u}_{j}\tilde{x}^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ 0_{q}\end{array}\right),

(28)

where

D=\left(\begin{array}[]{cc}D_{1}&0_{(n-q)\times q}\\ C_{l}^{\top}&C_{r}^{\top}\end{array}\right),~~\left(\begin{array}[]{cc}C_{l}^{\top}&C_{r}^{\top}\end{array}\right)=C^{\top},~~\tilde{x}=(x_{1},...,x_{n-q})^{\top},

(29)

then $(\tilde{\mathcal{N}},\tilde{u})$ is called a reconstruction of $(\mathcal{N},u)$ , and $D$ is termed as the reconstructing matrix from $(\mathcal{N},u)$ to $(\tilde{\mathcal{N}},\tilde{u})$ .

Remark 31.

If the MAS is not conservative, i.e., $q=0$ , then $D=D_{1}=\mathrm{diag}(d_{1},...,d_{n})$ , and (28) and (20) are equivalent. This means that reconstruction relation between two MASs will degenerate to the linear conjugacy relation in this case.

We then present the main theorem of this subsection, which characterizes the input-value set $\mathbb{U}$ using the concept of reconstruction.

Theorem 32.

For a time-varying MAS $(\mathcal{S},\mathcal{C},\mathcal{R},u)$ of (7), suppose there is another time-varying MAS $(\tilde{\mathcal{S}},\tilde{\mathcal{C}},\tilde{\mathcal{R}},\tilde{u})$ , which is weakly reversible and has toric locus $\mathcal{T}(\tilde{\mathcal{N}})$ with $\tilde{\mathcal{N}}=(\tilde{\mathcal{S}},\tilde{\mathcal{C}},\tilde{\mathcal{R}})$ , to be its reconstruction, i.e., there exists a positive, diagonal matrix $D_{1}=\mathrm{diag}(d_{1},...,d_{n-q})$ $(d_{i}>0)$ to support (28). Further, suppose that $\tilde{u}$ depends smoothly on $u\in\mathcal{U}$ with

\mathcal{U}=\{u\in\mathbb{R}^{r}_{>0}:~\exists~\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N}})~such~that~(\ref{eq_def_reconstruction})~holds\}.

(30)

To prove theorem 32, we need the concept of reverse reconstruction and some lemmas.

Definition 33 (Reverse reconstruction [31]).

For a conservative MAS $(\mathcal{N},\kappa)$ , let $(\tilde{\mathcal{N}},\tilde{\kappa})$ denote its reconstruction under the reconstructing matrix $D$ . A MAS $(\hat{\mathcal{N}},\hat{\kappa})$ is the reverse reconstruction of $(\mathcal{N},\kappa)$ with respect to $(\tilde{\mathcal{N}},\tilde{\kappa})$ if its complexes set $\hat{\mathcal{C}}=\cup_{j=1}^{\hat{r}}\{\hat{Z}_{\cdot j},\hat{Z}_{\cdot j}^{\prime}\}$ and reaction set $\hat{\mathcal{R}}=\cup_{j=1}^{\hat{r}}\{\hat{Z}_{\cdot j}\overset{\hat{\kappa}_{j}}{\longrightarrow}\hat{Z}_{\cdot j}^{\prime}\}$ satisfy

(i)

$\hat{r}=\tilde{r}$ ;
(ii)

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

As discussed in [31], although the underlying CRN in the reverse reconstruction does not fit into the CRN class defined in definition 1 (the elements in $\hat{Z}_{\cdot j}^{\prime}$ are not necessarily nonnegative integers), it does not invalidate subsequent results. In addition, reverse reconstruction has an important property: its dynamics are equivalent to those of the first $n-q$ species in the original network. This is stated precisely in the following lemma.

Lemma 34 ([31], Proposition 2).

For a conservative MAS $(\mathcal{N},\kappa)$ with a conserved matrix $C\in\mathbb{R}^{n\times q}_{>0}$ and an equilibrium $x^{*}\in\mathbb{R}^{n}_{>0}$ , let $(\tilde{\mathcal{N}},\tilde{\kappa})$ be its reconstruction bridged by the reconstructing matrix $D$ defined in (29), and let $(\hat{\mathcal{N}},\hat{\kappa})$ be its reverse reconstruction with respect to $(\tilde{\mathcal{N}},\tilde{\kappa})$ . Then $\hat{x}^{*}=(x^{*}_{1},...,x^{*}_{n-q})\in\mathbb{R}^{n-q}_{>0}$ is an equilibrium in $(\hat{\mathcal{N}},\hat{\kappa})$ , and the dynamics of $(\hat{\mathcal{N}},\hat{\kappa})$ is equivalent to that of $(\mathcal{N},\mathcal{\kappa})$ at the first $n-q$ species.

The following two lemmas are essential for the proof of theorem 32. We also put the proofs for these two lemmas in the Appendix to improve readability.

Lemma 35.

Suppose all the conditions in theorem 32 are satisfied. Denote $(\hat{\mathcal{N}},\hat{u})$ as the reverse reconstruction of $(\mathcal{N},u)$ with respect to $(\tilde{\mathcal{N}},\tilde{u})$ , and $\hat{\mathcal{U}}=\{\hat{u}=(u_{1},...,u_{\hat{r}}):u\in\mathcal{U}\}$ . Then for any $\hat{u}\in\hat{\mathcal{U}}$ , the MAS $(\hat{\mathcal{N}},\hat{u})$ has a unique equilibrium $\hat{x}^{*}=\hat{x}^{*}(\hat{u})$ in $(\hat{x}^{*}+\hat{\mathscr{S}})\cap\mathbb{R}^{n-q}_{>0}$ . Moreover, $\hat{x}^{*}=\hat{x}^{*}(\hat{u})$ depends smoothly on the parameter values $\hat{u}\in\hat{\mathcal{U}}$ .

Lemma 36.

Suppose all the conditions in lemma 35 are satisfied. Then for any compact set $\hat{\mathbb{U}}\subseteq\hat{\mathcal{U}}$ , there exists a function $V_{0}(\hat{x},\hat{u})$ satisfying the condition (1)-(3) in lemma 13.

Proof of theorem 32.

For any compact set $\mathbb{U}\subseteq\mathcal{U}$ and any $u^{*}\in\mathbb{U}$ , suppose that $(\tilde{\mathcal{N}},\tilde{u}^{*})$ is the reconstruction of $(\mathcal{N},u^{*})$ bridged by the reconstructing matrix $D$ , and $(\hat{\mathcal{N}},\hat{u}^{*})$ is the reverse reconstruction of $(\mathcal{N},u^{*})$ with respect to $(\tilde{\mathcal{N}},\tilde{u}^{*})$ . By lemma 35, lemma 36 and proposition 14, we know that for compact set $\hat{\mathbb{U}}=\{(u_{1},...,u_{\hat{r}}):u\in\mathbb{U}\}$ and any $\hat{u}^{*}\in\hat{\mathbb{U}}$ , the time-varying system

\dot{\hat{x}}=\hat{f}(\hat{x};\hat{u})=\sum_{j=1}^{\hat{r}}\hat{u}_{j}(t)\hat{x}^{\hat{v}_{\cdot j}}\left(\hat{v}_{\cdot j}^{\prime}-\hat{v}_{\cdot j}\right),

(31)

with input-value set $\hat{\mathbb{U}}$ , is semiglobal ISS with respect to $(\hat{x}^{*},\hat{u}^{*})$ , where $\hat{x}^{*}$ satisfies $\hat{f}(\hat{x}^{*};\hat{u}^{*})=0$ . For system (7) corresponding to $\mathcal{N}$ , denote its species concentrations $x=(x_{\bot}^{\top},x_{\top}^{\top})$ , where $x_{\bot}\in\mathbb{R}^{n-q}_{\geq 0},~x_{\top}\in\mathbb{R}^{q}_{\geq 0}$ . According to lemma 34 we know that the dynamics of $\hat{x}$ is equivalent to that of $x_{\bot}$ . Therefore, $x_{\bot}$ is $\mathbb{R}^{n-q}_{>0}$ -forward complete, and for any compact set $F\subseteq\mathbb{R}^{n}_{\geq 0}$ containing $x^{*}=\left((x_{\bot}^{*})^{\top},(x_{\top}^{*})^{\top}\right)^{\top}$ (it satisfies $f(x^{*};u^{*})=0$ and $x_{\bot}^{*}=\hat{x}^{*}$ ), there exist class $\mathcal{KL}$ function $\hat{\beta}$ and class $\mathcal{K}_{\infty}$ function $\hat{\gamma}$ such that

\begin{split}|x_{\bot}-x_{\bot}^{*}|&\leq\hat{\beta}\left(|x_{\bot}(0)-x_{\bot}^{*}|,t\right)+\hat{\gamma}\left(\|\hat{u}-\hat{u}^{*}\|\right)\\ &\leq\hat{\beta}\left(|x_{0}-x^{*}|,t\right)+\hat{\gamma}\left(\|u-u^{*}\|\right)\end{split}

(32)

for each $\mathbb{U}$ -value input $u(\cdot)$ and every initial condition $x_{0}\in F\cap(x^{*}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ and $\forall~t\geq 0$ such that $x(s)\in F~\forall s\in[0,t]$ . According to lemma 29, $x_{\top}$ is uniquely determined by $x_{\bot}$ and $x^{*}$ . Specifically, it satisfies $C_{l}^{\top}x_{\bot}+C_{r}^{\top}x_{\top}=C^{\top}x=C^{\top}x^{*}=C_{l}^{\top}x_{\bot}^{*}+C_{r}^{\top}x_{\top}^{*}$ , i.e.,

x_{\top}=C_{r}^{-\top}C_{l}^{\top}(x_{\bot}^{*}-x_{\bot})+x_{\top}^{*}.

(33)

Since $x_{\bot}$ is $\mathbb{R}^{n-q}_{>0}$ -forward complete, it follows that $x_{\top}$ is $\mathbb{R}^{q}_{>0}$ -forward complete, and thus $x$ is $\mathbb{R}^{n}_{>0}$ -forward complete. Moreover, by (33) we have $|x_{\top}-x_{\top}^{*}|\leq\left|C_{r}^{-\top}C_{l}^{\top}\right||x_{\bot}-x_{\bot}^{*}|$ . Combining this with (32), we obtain

\begin{split}|x-x^{*}|&\leq\sqrt{1+\left|C_{r}^{-\top}C_{l}^{\top}\right|^{2}}~|x_{\bot}-x_{\bot}^{*}|\\ &\leq\sqrt{1+\left|C_{r}^{-\top}C_{l}^{\top}\right|^{2}}\left(\hat{\beta}\left(|x_{0}-x^{*}|,t\right)+\hat{\gamma}\left(\|u-u^{*}\|\right)\right)\end{split}

(34)

By setting $\beta(s,t)=\sqrt{1+\left|C_{r}^{-\top}C_{l}^{\top}\right|^{2}}~\hat{\beta}\left(s,t\right)$ and $\gamma(s)=\sqrt{1+\left|C_{r}^{-\top}C_{l}^{\top}\right|^{2}}~\hat{\gamma}\left(s\right)$ , the proof is complete.

Example 37.

We use the MAS $(\mathcal{N},u)$ discussed in [31] to exhibit ISS, which takes

2X_{1}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\mathrm{u}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{1}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\mathrm{u}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{2}}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\mathrm{u}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{1}}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\mathrm{u}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{2}}}\mkern 9.0mu$}}}}$}}{}X_{1}+X_{2}\overset{u_{3}}{\longrightarrow}2X_{2}

with dynamics

\left(\begin{array}[]{c}\dot{x}_{1}\\ \dot{x}_{2}\end{array}\right)=\left(\begin{array}[]{ccc}-1&1&-1\\ 1&-1&1\end{array}\right)\left(\begin{array}[]{c}u_{1}x_{1}^{2}\\ u_{2}x_{1}x_{2}\\ u_{3}x_{1}x_{2}\end{array}\right)=\left(\begin{array}[]{c}-u_{1}x_{1}^{2}+(u_{2}-u_{3})x_{1}x_{2}\\ u_{1}x_{1}^{2}-(u_{2}-u_{3})x_{1}x_{2}\end{array}\right)

(35)

and a conserved matrix $C=(1,1)^{\top}$ . It is proved that the MAS

2X_{1}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{1}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{2}}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{1}}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\tilde{u}\mathrm{}{\vphantom{\mathrm{X}}}_{\smash[t]{\mathrm{2}}}\mkern 9.0mu$}}}}$}}{}X_{1}

is the reconstruction of $(\mathcal{N},u)$ under the reconstructing matrix

D=\left(\begin{array}[]{cc}2&0\\ 1&1\end{array}\right)=\left(\begin{array}[]{cc}D_{1}&0\\ C_{l}^{\top}&C_{r}^{\top}\end{array}\right).

Note that the second MAS is weakly reversible, and has deficiency $\delta=2-1-1=0$ , so it follows that $\mathcal{T}(\tilde{\mathcal{N}})=\mathbb{R}^{2}_{>0}$ . Let $x^{*}\in\mathbb{R}^{2}_{>0}$ be an equilibrium of $(\mathcal{N},u)$ , then (28) is given by

\begin{split}&2\left(-u_{1}x_{1}^{2}+u_{2}x_{1}(x_{1}^{*}+x_{2}^{*}-x_{1})-u_{3}x_{1}(x_{1}^{*}+x_{2}^{*}-x_{1})\right)=-\tilde{u}_{1}x_{1}^{2}+\tilde{u}_{2}x_{1}.\\ \Longleftrightarrow~&2\left(-u_{1}-u_{2}+u_{3}\right)x_{1}^{2}+C^{*}\left(u_{2}-u_{3}\right)x_{1}=-\tilde{u}_{1}x_{1}^{2}+\tilde{u}_{2}x_{1},~\forall x\in(x^{*}+\mathscr{S})\cap\mathbb{R}^{2}_{\geq 0}.\\ \Longleftrightarrow~&\begin{cases}2\left(-u_{1}-u_{2}+u_{3}\right)=-\tilde{u}_{1},\\ C^{*}\left(u_{2}-u_{3}\right)=\tilde{u}_{2}.\end{cases}\Longleftrightarrow~\begin{cases}\tilde{u}_{1}=2(u_{1}+u_{2}-u_{3}),\\ \tilde{u}_{2}=C^{*}(u_{2}-u_{3}).\end{cases}\end{split}

where $C^{*}=2(x_{1}^{*}+x_{2}^{*})$ is a constant. Therefore, whenever $u\in\mathbb{R}^{3}_{>0}$ and $u_{2}>u_{3}$ , there exists a $\tilde{u}\in\mathbb{R}^{2}_{>0}=\mathcal{T}(\tilde{\mathcal{N}})$ such that (28) holds, which implies that $\mathcal{U}=\{u\in\mathbb{R}^{3}_{>0}:u_{2}>u_{3}\}$ . Set $u^{*}=(1,4,3)^{\top}\in\mathcal{U}$ and $x^{*}=(1,1)^{\top}$ , then by theorem 32 we know that for any compact set $\mathbb{U}\subseteq\mathcal{U}$ containing $u^{*}$ , the system (35), with input-value set $\mathbb{U}$ , is semiglobal ISS with respect to $(x^{*},u^{*})$ .

Unlike linear conjugacy, reconstruction generally involves MASs with different dimensions before and after the transformation. Consequently, analyzing (30) becomes highly challenging, particularly in characterizing the conditions under which $\mathcal{U}$ is not empty or is $\mathbb{R}^{r}_{>0}$ . This remains an important direction for future research.

As in the previous subsection, the results in this subsection are established under the assumption that $D$ and $\tilde{\mathcal{N}}$ are given. If only $\mathcal{N}$ is given, finding $\tilde{\mathcal{N}}$ typically requires solving a mixed-integer programming problem. A feasible algorithm can be found in [31].

5 Application

Recently, CRNs have been widely used as a mathematical framework for studying molecular computations [2, 18]. In this framework, the initial concentrations of certain species (i.e., $x_{i}(0)$ ) are treated as inputs, while the limiting steady state concentrations of other species (i.e., $\lim_{t\to\infty}x_{j}(t)$ ), are interpreted as outputs. By designing the resulting input-output relation to match a prescribed target function, one can theoretically implement the molecular computation of that function. Among these studies, an important question is: if two MASs $(\mathcal{N}_{1},\kappa_{1})$ and $(\mathcal{N}_{2},\kappa_{2})$ serve to compute functions $u=\sigma_{1}(w)$ and $x=\sigma_{2}(u)$ , respectively, can their coupling compute the corresponding composite function $x=\sigma_{2}\circ\sigma_{1}(w)$ ? If the answer is affirmative, then $(\mathcal{N}_{1},\kappa_{1})$ and $(\mathcal{N}_{2},\kappa_{2})$ are said to be “dynamically composable” [27].

Using the assumptions and techniques developed in [27], we can view the coupling dynamics as a time-varying system (7), where the input $u(t)$ is generated by the dynamics of $(\mathcal{N}_{1},\kappa_{1})$ , and the key issue in establishing dynamical composability reduces to Question 2 in Section 2. The following theorem provides a partial answer to this question.

Theorem 38.

Consider a MAS $(\mathcal{S},\mathcal{C},\mathcal{R},u)$ of (7) starting from $x(0)=x_{0}$ , and suppose $x^{*}$ is its unique equilibrium in $(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ under $u(t)\equiv u^{*}\in\mathbb{R}^{r}_{>0}$ . Then the limiting state of (7) satisfies $\lim_{t\to\infty}x(t)=x^{*}$ if

(1)

$u(t)\in\mathbb{U},~\forall~t\geq 0$ , where $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ is a compact set, and $\lim_{t\to\infty}u(t)=u^{*}$ ;
(2)

the MAS of (7), with input-value set $\mathbb{U}$ , is semiglobal ISS with respect to $(x^{*},u^{*})$ ;
(3)

the MAS of (7) is permanent, that is, there exists $\varepsilon>0$ such that for any initial value $x_{0}\in\mathbb{R}^{n}_{>0}$ , the solution of (7) satisfies $\liminf_{t\to\infty}x_{i}(t)>\varepsilon$ and $\limsup_{t\to\infty}x_{i}(t)<1/\varepsilon$ for all $i=1,2,...,n$ .

Proof.

By condition (3), there exists a compact set $F\subseteq\mathbb{R}^{n}_{>0}$ such that $x(t)\in F,~\forall~t\geq 0$ . Using condition (2) and treating $s\geq 0$ as the initial time, there exist a class $\mathcal{KL}$ function $\beta$ and a class $\mathcal{K}_{\infty}$ function $\gamma$ such that

|x(t)-x^{*}|\leq\beta(|x(s)-x^{*}|,t-s)+\gamma\left(\mathrm{ess.sup}\{|u(\tau)-u^{*}|:\tau\geq s\}\right),

(36)

where $0\leq s\leq t$ . Substituting $s=\frac{t}{2}$ into (36) yields

|x(t)-x^{*}|\leq\beta\left(|x(\frac{t}{2})-x^{*}|,\frac{t}{2}\right)+\gamma\left(\mathrm{ess.sup}\{|u(\tau)-u^{*}|:\tau\geq\frac{t}{2}\}\right).

(37)

By $u(t)\to u^{*}$ we obtain $\gamma\left(\mathrm{ess.sup}\{|u(t)-u^{*}|:\tau\geq\frac{t}{2}\}\right)\to 0$ as $t\to\infty$ . To estimate $|x(\frac{t}{2})-x^{*}|$ , apply (36) with $s=0$ and $t$ replaced by $\frac{t}{2}$ , then we have

|x(\frac{t}{2})-x^{*}|\leq\beta\left(|x_{0}-x^{*}|,\frac{t}{2}\right)+\gamma\left(\mathrm{ess.sup}\{|u(\tau)-u^{*}|:\tau\geq 0\}\right).

(38)

Since $u(t)\in\mathbb{U}$ is bounded, there exists a constant $M>0$ such that

\gamma\left(\mathrm{ess.sup}\{|u(\tau)-u^{*}|:\tau\geq 0\}\right)\leq M.

(39)

In addition, since $\beta$ is a class $\mathcal{KL}$ function, for sufficiently large $t$ , it holds

\beta\left(|x_{0}-x^{*}|,\frac{t}{2}\right)\leq M.

(40)

Substituting eqs. 38, 39, and 40 into (37), we obtain

\begin{split}|x(t)-x^{*}|&\leq\beta\left(2M,\frac{t}{2}\right)+\gamma\left(\mathrm{ess.sup}\{|u(\tau)-u^{*}|:\tau\geq\frac{t}{2}\}\right)\\ &\to 0\quad(t\to\infty)\end{split}

(41)

that is, $\lim_{t\to\infty}x(t)=x^{*}$ .

theorem 38 provides a sufficient condition for the global asymptotic stability of the time-varying system. The permanence condition involved in the assumptions has been widely used in the analysis of global asymptotic stability for MASs [15, 22].

Example 39.

Revisit the CRNs in Example 23. Consider an input $u:\mathbb{R}_{\geq 0}\to\mathbb{R}^{r}_{>0}$ that satisfies $\lim_{t\to\infty}u(t)=u^{*}$ . Then $u(t)$ is bounded, and there exists a compact set $\mathbb{U}\subseteq\mathbb{R}^{r}_{>0}$ such that $u(t)\in\mathbb{U},~\forall~t\geq 0$ . That is, the condition (1) in theorem 38 holds. According to the discussion in Example 23, the condition (2) in theorem 38 holds.

To obtain the permanence of time-varying system (23) with initial value $x(0)=x_{0}$ and $u(t)\in\mathbb{U}$ , consider the CRN $\tilde{\mathcal{N}}$ corresponding to time-varying system (24) with initial value $\tilde{x}(0)=\mathrm{diag}(1,2)\cdot x_{0}$ , where $(\tilde{u}_{1},\tilde{u}_{2},\tilde{u}_{3},\tilde{u}_{4})=\left(\frac{1}{2}u_{1},\frac{1}{2}u_{2},u_{3},2u_{4}\right)$ . Then it can be directly verified that $x(t)=\mathrm{diag}(1,1/2)\cdot\tilde{x}(t),~\forall~t\geq 0$ , and $\tilde{u}(t)\in\tilde{\mathbb{U}}$ , where $\tilde{\mathbb{U}}$ is a compact subset of $\mathbb{R}^{4}_{>0}$ . In addition, since $\tilde{\mathcal{N}}$ is weakly reversible, by Theorem 6.4 in [15] we know that $\tilde{x}(t)$ is permanent. Therefore, $x(t)$ is also permanent. That is, the condition (3) in theorem 38 holds.

Since all the conditions in theorem 38 are satisfied, we have $\lim_{t\to\infty}x(t)=x^{*}$ .

Example 39 shows that the time-varying system (23) satisfies $x(t)\to x^{*}$ if $u(t)\to u^{*}$ . Seeing $u$ and $x$ as the upstream and downstream computational modules in molecular computations, respectively, this property indicates that the coupled system can implement a layer-by-layer computation through parallel chemical reactions. Specifically, consider an upstream computational module

U_{2}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\mathrm{\text{$\text{$1$}$}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\mathrm{\text{$\text{$1$}$}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\mathrm{\text{$\text{$1$}$}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\mathrm{\text{$\text{$1$}$}}\mkern 9.0mu$}}}}$}}{}U_{3},~U_{2}\overset{1}{\longrightarrow}U_{2}+U_{4},~2U_{4}\overset{1}{\longrightarrow}\varnothing

(42)

with dynamics to be $(\dot{u}_{1},\dot{u}_{2},\dot{u}_{3},\dot{u}_{4})^{\top}=(0,u_{3}-u_{2},u_{2}-u_{3},u_{2}-2u_{4}^{2})^{\top}$ . This module implements the molecular computation of the function $v=\sigma_{1}(u)=\left(u_{1},\frac{u_{2}+u_{3}}{2},\frac{u_{2}+u_{3}}{2},\frac{\sqrt{u_{2}+u_{3}}}{2}\right)^{\top}$ . The downstream computational module

\begin{split}&X_{1}+2X_{2}+U_{1}\overset{1}{\longrightarrow}X_{1}+3X_{2}+U_{1},~2X_{1}+U_{4}\overset{1}{\longrightarrow}X_{2}+U_{4},\\ &X_{1}+3X_{2}+U_{2}\overset{1}{\longrightarrow}X_{1}+X_{2}+U_{2},~X_{1}+X_{2}+U_{3}\overset{1}{\longrightarrow}3X_{1}+U_{3},\end{split}

(43)

which has dynamics (23), implements the molecular computation of $x=\sigma_{2}(v)=\left(\frac{v_{1}v_{3}}{2v_{2}v_{4}},\frac{v_{1}}{2v_{2}}\right)^{\top}$ . By Example 39 it is concluded that the coupled system, which is given by

\begin{split}&U_{2}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\mathrm{\text{$\text{$1$}$}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\mathrm{\text{$\text{$1$}$}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\mathrm{\text{$\text{$1$}$}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\mathrm{\text{$\text{$1$}$}}\mkern 9.0mu$}}}}$}}{}U_{3},~U_{2}\overset{1}{\longrightarrow}U_{2}+U_{4},~2U_{4}\overset{1}{\longrightarrow}\varnothing,\\ &X_{1}+2X_{2}+U_{1}\overset{1}{\longrightarrow}X_{1}+3X_{2}+U_{1},~2X_{1}+U_{4}\overset{1}{\longrightarrow}X_{2}+U_{4},\\ &X_{1}+3X_{2}+U_{2}\overset{1}{\longrightarrow}X_{1}+X_{2}+U_{2},~X_{1}+X_{2}+U_{3}\overset{1}{\longrightarrow}3X_{1}+U_{3},\end{split}

(44)

with dynamics to be

\left(\begin{array}[]{c}\dot{u}_{1}\\ \dot{u}_{2}\\ \dot{u}_{3}\\ \dot{u}_{4}\\ \dot{x}_{1}\\ \dot{x}_{2}\end{array}\right)=\left(\begin{array}[]{c}0\\ u_{3}-u_{2}\\ u_{2}-u_{3}\\ u_{2}-2u_{4}^{2}\\ 2u_{3}x_{1}x_{2}-2u_{4}x_{1}^{2}\\ u_{1}x_{1}x_{2}^{2}-2u_{2}x_{1}x_{2}^{3}-u_{3}x_{1}x_{2}+u_{4}x_{1}^{2}\end{array}\right),

(45)

can implement the molecular computation of the composite function $\sigma_{2}\circ\sigma_{1}$ , i.e., $(x_{1},x_{2})^{\top}=\sigma_{2}(\sigma_{1}(u))=\left(\frac{u_{1}}{\sqrt{u_{2}+u_{3}}},\frac{u_{1}}{u_{2}+u_{3}}\right)^{\top}$ . This suggests the computation by chemical reactions in (44) follows the same result as given by cascading the computational result of reactions in (42) and that of reactions in (43). To observe the parallel computation result more visually, fig. 2 presents the numerical simulation results for this molecular computation. It can be observed that both $x_{1}$ and $x_{2}$ converge to their corresponding concentrations, verifying the effectiveness of the parallel molecular computation.

6 Conclusions

In this paper, we investigate the ISS of CRNs with time-varying reaction-rate inputs, building on the work of Chaves and Sontag. For mass-action kinetics, we first extend the established ISS results for weakly reversible networks with zero deficiency and single linkage class to a broader class of weakly reversible networks. In particular, by restricting the input-value set to a compact subset of the toric locus, we establish semiglobal ISS for weakly reversible networks that have arbitrary deficiency and linkage classes. Then, by exploiting the concepts of linear conjugacy and reconstruction, we further extend these ISS results to a class of networks that are not weakly reversible, even though their toric locus are empty. Finally, we apply the above results to molecular computations. We show that, when the reaction-rate input converges to a constant value and the corresponding trajectory is permanent, the semiglobal ISS property guarantees convergence of the state to the equilibrium associated with the limiting input, thereby providing a theoretical foundation for the implementation of parallel molecular computations. A concrete example illustrates the application of the proposed theory to molecular computations.

Several directions deserve further investigation. First, according to the results established in Section 4, the set $\mathcal{U}$ determines the admissible range of the input $u(\cdot)$ . Therefore, it is necessary to develop efficient computational methods for determining an admissible, or even the maximal admissible, input set $\mathcal{U}$ for a given CRN. Second, the notion of ISS can be extended to CRNs with time delays, and one may further investigate structural conditions under which delayed networks possess the ISS property. In addition, extending the analysis to other classes of kinetics, such as Michaelis–Menten kinetics, would also be of both theoretical and practical significance.

Appendix A Proofs

Proof of lemma 21.

We complete the proof in three steps.

Step 1. We prove the existence of the equilibrium in $(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ . For any $u\in\mathcal{U}$ , we know that there exists a $\tilde{u}$ such that $\forall~v\in\mathcal{C}_{react}$ it holds (21) and $(\tilde{\mathcal{N}},\tilde{u})$ is complex balanced, which means there exists an equilibrium in each positive stoichiometric compatibility class. Specifically, suppose $\tilde{x}^{*}\in(Dx_{0}+\tilde{\mathscr{S}})\cap\mathbb{R}^{n}_{>0}$ is the equilibrium of $(\tilde{\mathcal{N}},\tilde{u})$ and thus satisfies $\sum_{j=1}^{\tilde{r}}\tilde{u}_{j}(\tilde{x}^{*})^{\tilde{v}_{\cdot j}}(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j})=0$ . Then we have

\begin{split}&D\sum_{j=1}^{r}u_{j}\left(D^{-1}\tilde{x}^{*}\right)^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)\\ =&D\sum_{j=1}^{r}u_{j}\left(\tilde{x}^{*}\right)^{v_{\cdot j}}\left(\prod_{i=1}^{n}d_{i}^{-v_{ij}}\right)\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)\\ =&\sum_{v\in\mathcal{C}_{react}}(\tilde{x}^{*})^{v}\left(\prod_{i=1}^{n}d_{i}^{-v_{i}}\right)D\sum_{\{j~:~v_{\cdot j}=v\}}u_{j}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)\\ =&\sum_{v\in\mathcal{C}_{react}}(\tilde{x}^{*})^{v}\left(\prod_{i=1}^{n}d_{i}^{-v_{i}}\right)\sum_{\{j~:~\tilde{v}_{\cdot j}=v\}}\left(\tilde{u}_{j}\prod_{i=1}^{n}d_{i}^{v_{i}}\right)\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ =&\sum_{v\in\mathcal{C}_{react}}(\tilde{x}^{*})^{v}\sum_{\{j~:~\tilde{v}_{\cdot j}=v\}}\tilde{u}_{j}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ =&\sum_{j=1}^{\tilde{r}}\tilde{u}_{j}(\tilde{x}^{*})^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)=0\end{split}

(46)

with the third equality holding due to (21). Since $D$ is invertible, we can get $\sum_{j=1}^{r}u_{j}\left(x^{*}\right)^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=0$ , where $x^{*}=D^{-1}\tilde{x}^{*}$ . This means that $x^{*}$ is an equilibrium of $(\mathcal{N},u)$ . In addition, by the definition of linear conjugacy, we have $\Phi(x_{0},t)=D^{-1}\tilde{\Phi}(Dx_{0},t)$ . As $\tilde{x}^{*}\in(Dx_{0}+\tilde{\mathscr{S}})\cap\mathbb{R}^{n}_{>0}$ , we have $x^{*}\in(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ .

Step 2. We prove the uniqueness of the equilibrium in $(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ . Suppose that $x^{**}\in(x_{0}+\mathscr{S})\cap\mathbb{R}^{n}_{>0}$ satisfies $\sum_{j=1}^{r}u_{j}\left(x^{**}\right)^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=0$ . By an argument similar to that in Step 1, we obtain $\sum_{j=1}^{\tilde{r}}\tilde{u}_{j}\left(Dx^{**}\right)^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)=0$ , and $Dx^{**}\in(Dx_{0}+\tilde{\mathscr{S}})\cap\mathbb{R}^{n}_{>0}$ . Since $(\tilde{\mathcal{N}},\tilde{u})$ is complex balanced, there is a unique equilibrium in each positive stoichiometric compatibility class. Thus we have $Dx^{**}=\tilde{x}^{*}=Dx^{*}$ , i.e., $x^{**}=x^{*}$ .

Step 3. We establish the smooth dependence property. Since $x^{*}=D^{-1}\tilde{x}^{*}$ , we can obtain that $x^{*}$ depends smoothly on $\tilde{x}^{*}$ . According to the Corollary 3.13 in [12], we know that $\tilde{x}^{*}=\tilde{x}^{*}(\tilde{u})$ depends smoothly on $\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N}})$ . From the assumption in theorem 20 it follows that $\tilde{u}$ depends smoothly on $u\in\mathcal{U}$ . Therefore, we know that $x^{*}=x^{*}(u)$ depends smoothly on the parameter values $u\in\mathcal{U}$ .

Proof of lemma 22.

For any compact set $\mathbb{U}\subseteq\mathcal{U}$ , denote

V_{0}(x,u)=\sum_{i=1}^{n}d_{i}(x_{i}(\ln x_{i}-\ln x_{i}^{*}-1)+x_{i}^{*})

, where $x^{*}=x^{*}(u),~u\in\mathbb{U}$ is defined as in lemma 21.

According to lemma 21, we know that $x^{*}=x^{*}(u)$ depends smoothly on the parameter values $u\in\mathbb{U}$ . Therefore, the restriction of $V_{0}$ to $\mathbb{R}^{n}_{>0}\times\mathbb{U}$ is continuously differentiable, and the condition (1) in lemma 13 holds.

For any $u^{*}\in\mathbb{U}$ , it is straightforward to verify that $V_{0}(x,u^{*})\geq 0$ , with the equality holding if and only if $x=x^{*}(u^{*})$ . Denote

\underline{\alpha}_{u^{*}}(s)=\inf_{|x-x^{*}|\geq s}V_{0}(x,u^{*}),\quad\overline{\alpha}_{u^{*}}(s)=\sup_{|x-x^{*}|\leq s}V_{0}(x,u^{*}),\quad s\geq 0,

which satisfy $\underline{\alpha}_{u^{*}}(0)=\overline{\alpha}_{u^{*}}(0)=0$ , and are continuous, positive definite, and strictly increasing (and thus class $\mathcal{K}$ functions). Then it holds

\underline{\alpha}_{u^{*}}(|x-x^{*}|)\leq V_{0}(x,u^{*})\leq\overline{\alpha}_{u^{*}}(|x-x^{*}|),~x\in\mathbb{R}^{n}_{\geq 0}.

Moreover, since $V_{0}(x,u^{*})\to+\infty$ as $|x-x^{*}|\to+\infty$ , we have $\underline{\alpha}_{u^{*}}(s)\to+\infty,~\overline{\alpha}_{u^{*}}(s)\to+\infty$ . Therefore, $\underline{\alpha}_{u^{*}}$ and $\overline{\alpha}_{u^{*}}$ are class $\mathcal{K}_{\infty}$ functions, and the condition (2) in lemma 13 holds.

For any $u^{*}\in\mathbb{U}$ , by the proof of lemma 21, there exists a unique equilibrium $x^{*}=x^{*}(u^{*})$ for the MAS $\dot{x}=f(x;u^{*})$ in each positive stoichiometric compatibility class. Moreover, there exists a $\tilde{u}^{*}\in\mathcal{T}(\tilde{\mathcal{N}})$ such that $x^{*}=D^{-1}\tilde{x}^{*}$ , where $\tilde{x}^{*}=\tilde{x}^{*}(\tilde{u}^{*})$ is the complex balanced equilibrium of system $\dot{\tilde{x}}=\tilde{f}(\tilde{x};\tilde{u}^{*})=\sum_{j=1}^{\tilde{r}}\tilde{u}_{j}\tilde{x}^{\tilde{v}_{\cdot j}}(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j})$ . As $\tilde{x}^{*}$ is the complex balanced equilibrium, let $\tilde{V}_{0}(\tilde{x},\tilde{u}^{*})=\sum_{i=1}^{n}(\tilde{x}_{i}(\ln\tilde{x}_{i}-\ln\tilde{x}_{i}^{*}-1)+\tilde{x}_{i}^{*})$ , then it follows $\nabla^{\top}_{\tilde{x}}\tilde{V}_{0}(\tilde{x},\tilde{u}^{*})\tilde{f}(\tilde{x};\tilde{u}^{*})\leq 0$ , with the equality holding if and only if $\tilde{x}=\tilde{x}^{*}$ . Since

\begin{split}\nabla_{x}^{\top}V_{0}(x,u^{*})&=\left(d_{1}(\ln x_{1}-\ln x_{1}^{*}),...,d_{n}(\ln x_{n}-\ln x_{n}^{*})\right)\\ &=(\ln(d_{1}x_{1})-\ln\tilde{x}_{1}^{*},...,\ln(d_{n}x_{n})-\ln\tilde{x}_{n}^{*})~D\\ &=\nabla_{\tilde{x}}^{\top}\tilde{V}_{0}(Dx,\tilde{u}^{*})~D\end{split}

and

\begin{split}f(x;u^{*})&=\sum_{j=1}^{r}u_{j}^{*}x^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)\\ &=\sum_{v\in\mathcal{C}_{react}}x^{v}\sum_{\{j|v_{\cdot j}=v\}}u_{j}^{*}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)\\ &=\sum_{v\in\mathcal{C}_{react}}x^{v}D^{-1}\sum_{\{j|\tilde{v}_{\cdot j}=v\}}\left(\tilde{u}_{j}^{*}\prod_{i=1}^{n}d_{i}^{v_{i}}\right)\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ &=D^{-1}\sum_{v\in\mathcal{C}_{react}}(Dx)^{v}\sum_{\{j|\tilde{v}_{\cdot j}=v\}}\tilde{u}_{j}^{*}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ &=D^{-1}\sum_{j=1}^{\tilde{r}}\tilde{u}_{j}^{*}(Dx)^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ &=D^{-1}\tilde{f}(Dx;\tilde{u}^{*})\end{split}

we thus have $\nabla_{x}^{\top}V_{0}(x,u^{*})f(x;u^{*})=\nabla_{\tilde{x}}^{\top}\tilde{V}_{0}(Dx,\tilde{u}^{*})\tilde{f}(Dx;\tilde{u}^{*})\leq 0$ , with the equality holding if and only if $Dx=\tilde{x}^{*}$ , i.e., $x=D^{-1}\tilde{x}^{*}=x^{*}$ , and the condition (3) in lemma 13 holds.

Proof of lemma 35.

For any $\hat{u}\in\hat{\mathcal{U}}$ , since $(\hat{\mathcal{N}},\hat{u})$ is the reverse reconstruction of $(\mathcal{N},u)$ with respect to $(\tilde{\mathcal{N}},\tilde{u})$ , by (30) it follows that $(\tilde{\mathcal{N}},\tilde{u})$ is complex balanced. Then according to Theorem 5 in [31] we know that there exists a unique equilibrium $\hat{x}^{*}=\hat{x}^{*}(\hat{u})$ in $(\hat{x}_{0}+\hat{\mathscr{S}})\cap\mathbb{R}^{n-q}_{>0}$ , and $\hat{x}^{*}$ depends smoothly on $\tilde{x}^{*}$ , where $\tilde{x}^{*}$ is the complex balanced equilibrium of $(\tilde{\mathcal{N}},\tilde{u})$ . Further, by an argument similar to that used in the proof of lemma 21, we obtain that $\hat{x}^{*}=\hat{x}^{*}(\hat{u})$ depends smoothly on the parameter values $\hat{u}\in\hat{\mathcal{U}}$ .

Proof of lemma 36.

For any compact set $\hat{\mathbb{U}}\subseteq\hat{\mathcal{U}}$ , denote

V_{0}(\hat{x},\hat{u})=\sum_{i=1}^{n-q}d_{i}(\hat{x}_{i}(\ln\hat{x}_{i}-\ln\hat{x}_{i}^{*}-1)+\hat{x}_{i}^{*})

, where $\hat{x}^{*}=\hat{x}^{*}(\hat{u}),~\hat{u}\in\hat{\mathbb{U}}$ is defined as in lemma 35.

According to lemma 35, we know that $\hat{x}^{*}=\hat{x}^{*}(\hat{u})$ depends smoothly on the parameter values $\hat{u}\in\hat{\mathbb{U}}$ . Therefore, the restriction of $V_{0}$ to $\mathbb{R}^{n-q}_{>0}\times\hat{\mathbb{U}}$ is continuously differentiable, and the condition (1) in lemma 13 holds.

For any $\hat{u}^{*}\in\hat{\mathbb{U}}$ , it is straightforward to verify that $V_{0}(\hat{x},\hat{u}^{*})\geq 0$ , with the equality holding if and only if $\hat{x}=\hat{x}^{*}(\hat{u}^{*})$ . Similarly to the proof of lemma 22, one can show that the condition (2) in lemma 13 holds.

For any $\hat{u}^{*}\in\hat{\mathbb{U}}$ , there exist $u^{*}\in\mathcal{U}$ and $\tilde{u}^{*}\in\mathcal{T}(\tilde{\mathcal{N}})$ such that $\hat{u}^{*}=(u^{*}_{1},...,u^{*}_{\hat{r}})$ , and

D\sum_{j=1}^{r}u^{*}_{j}x^{v_{\cdot j}}\left(v_{\cdot j}^{\prime}-v_{\cdot j}\right)=\left(\begin{array}[]{c}\sum_{j=1}^{\tilde{r}}\tilde{u}^{*}_{j}\tilde{x}^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ 0_{q}\end{array}\right).

By lemma 34, the above expression can be rewritten as

D_{1}\sum_{j=1}^{\hat{r}}\hat{u}^{*}_{j}\hat{x}^{\hat{v}_{\cdot j}}\left(\hat{v}_{\cdot j}^{\prime}-\hat{v}_{\cdot j}\right)=\sum_{j=1}^{\tilde{r}}\tilde{u}^{*}_{j}\tilde{x}^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right),

(47)

and $\hat{x}^{*}=(x_{1}^{*},...,x_{n-q}^{*})\in\mathbb{R}^{n-q}_{>0}$ is an equilibrium of $(\hat{\mathcal{N}},\hat{u}^{*})$ , where $x^{*}\in\mathbb{R}^{n}_{>0}$ is an equilibrium of $(\mathcal{N},u^{*})$ . According to lemma 35, $\hat{x}^{*}$ is the unique equilibrium in $(\hat{x}^{*}+\hat{\mathscr{S}})\cap\mathbb{R}^{n-q}_{>0}$ . In addition, it follows immediately from (47) that $\tilde{x}^{*}=\hat{x}^{*}$ is an equilibrium of $(\tilde{\mathcal{N}},\tilde{u}^{*})$ . Since $\tilde{u}^{*}\in\mathcal{T}(\tilde{\mathcal{N}})$ , $\tilde{x}^{*}$ is a complex balanced equilibrium and thus the unique equilibrium in $(\tilde{x}^{*}+\tilde{\mathscr{S}})\cap\mathbb{R}^{n-q}_{>0}$ . Further, denote $V_{1}(\tilde{x},\tilde{u}^{*})=\sum_{i=1}^{n-q}(\tilde{x}_{i}(\ln\tilde{x}_{i}-\ln\tilde{x}_{i}^{*}-1)+\tilde{x}_{i}^{*})$ , then we have

\nabla_{\tilde{x}}^{\top}V_{1}(\tilde{x},\tilde{u}^{*})\sum_{j=1}^{\tilde{r}}\tilde{u}^{*}_{j}\tilde{x}^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\leq 0,

(48)

with the equality holding if and only if $\tilde{x}=\tilde{x}^{*}$ . Therefore, it follows

\begin{split}\nabla^{\top}_{\hat{x}}V_{0}(\hat{x},\hat{u}^{*})\sum_{j=1}^{\hat{r}}\hat{u}^{*}_{j}\hat{x}^{\hat{v}_{\cdot j}}\left(\hat{v}_{\cdot j}^{\prime}-\hat{v}_{\cdot j}\right)&=\nabla^{\top}_{\hat{x}}V_{1}(\hat{x},\hat{u}^{*})D_{1}\sum_{j=1}^{\hat{r}}\hat{u}^{*}_{j}\hat{x}^{\hat{v}_{\cdot j}}\left(\hat{v}_{\cdot j}^{\prime}-\hat{v}_{\cdot j}\right)\\ &=\nabla^{\top}_{\tilde{x}}V_{1}(\tilde{x},\tilde{u}^{*})\sum_{j=1}^{\tilde{r}}\tilde{u}^{*}_{j}\tilde{x}^{\tilde{v}_{\cdot j}}\left(\tilde{v}_{\cdot j}^{\prime}-\tilde{v}_{\cdot j}\right)\\ &\leq 0\end{split}

(49)

with the equality holding if and only if $\hat{x}=\tilde{x}^{*}=\hat{x}^{*}$ , and the condition (3) in lemma 13 holds.

Appendix B An example of a non-weakly reversible biological network with practical significance

Consider the p21-activated kinase 1 (PAK1) network

2X_{1}\overset{u_{1}}{\longrightarrow}X_{1}+X_{2},\quad X_{1}\overset{u_{2}}{\longleftarrow}X_{2}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\mathrm{\text{$u_{3}$}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\mathrm{\text{$u_{4}$}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\mathrm{\text{$u_{3}$}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\mathrm{\text{$u_{4}$}}\mkern 9.0mu$}}}}$}}{}X_{3},

where $X_{1},X_{2},X_{3}$ denote unphosphorylated PAK1, monophosphorylated PAK1, and diphosphorylated PAK1, respectively. This network plays a crucial role in the regulation of cell motility and morphology [24]. The dynamics of this time-varying MAS is

\left(\begin{array}[]{c}\dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\end{array}\right)=\left(\begin{array}[]{c}-u_{1}x_{1}^{2}+u_{2}x_{2}\\ u_{1}x_{1}^{2}-u_{2}x_{2}-u_{3}x_{2}+u_{4}x_{3}\\ u_{3}x_{2}-u_{4}x_{3}\end{array}\right).

(50)

This MAS was said to be linear conjugate to the following MAS $(\tilde{\mathcal{N}},\tilde{u})$

2X_{1}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\mathrm{\text{$\tilde{u}_{1}$}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\mathrm{\text{$\tilde{u}_{2}$}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\mathrm{\text{$\tilde{u}_{1}$}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\mathrm{\text{$\tilde{u}_{2}$}}\mkern 9.0mu$}}}}$}}{}X_{2}{}\mathrel{\hbox to0.0pt{\raisebox{0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 9.0mu{}\mathrm{\text{$\tilde{u}_{3}$}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu\hphantom{{}\mathrm{\text{$\tilde{u}_{4}$}}}\mkern 5.0mu$}}}}$}\hss}\raisebox{-0.94722pt}{$\mathrel{\mathop{\makebox[0.0pt]{\to}}\limits^{\mkern 5.0mu\hphantom{{}\mathrm{\text{$\tilde{u}_{3}$}}}\mkern 9.0mu}_{\makebox{\raisebox{3.76735pt}[0.0pt]{$\scriptstyle\mkern 9.0mu{}\mathrm{\text{$\tilde{u}_{4}$}}\mkern 9.0mu$}}}}$}}{}X_{3}

with dynamics

\left(\begin{array}[]{c}\dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\end{array}\right)=\left(\begin{array}[]{c}-2\tilde{u}_{1}x_{1}^{2}+2\tilde{u}_{2}x_{2}\\ \tilde{u}_{1}x_{1}^{2}-\tilde{u}_{2}x_{2}-\tilde{u}_{3}x_{2}+\tilde{u}_{4}x_{3}\\ \tilde{u}_{3}x_{2}-\tilde{u}_{4}x_{3}\end{array}\right),

(51)

linked by $D=\mathrm{diag}(2,1,1)$ . Clearly, the second MAS is weakly reversible, and has deficiency $\delta=3-1-2=0$ , so it follows $\mathcal{T}(\tilde{\mathcal{N}})=\mathbb{R}^{4}_{>0}$ . Further, for $z=(2,0,0)$ we have

\begin{split}D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)&=\left\{u_{1}(-2,1,0)^{\top}:u\in\mathbb{R}^{4}_{>0}\right\}\\ C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)&=\left\{4\tilde{u}_{1}(-2,1,0)^{\top}:\tilde{u}\in\mathcal{T}(\tilde{\mathcal{N}})=\mathbb{R}^{4}_{>0}\right\}\end{split}

and thus $D~C_{\mathcal{R},\mathbb{R}^{r}_{>0}}(z)\subseteq C_{\tilde{\mathcal{R}},\mathcal{T}(\tilde{\mathcal{N}})}(z)$ . Similarly, one can verify that the above inclusion relation holds for any $z\in\mathcal{C}_{react}$ . By proposition 25 we know that the set $\mathcal{U}$ defined in (22) satisfies $\mathcal{U}=\mathbb{R}^{4}_{>0}$ . In addition, by setting $u^{*}=(1,4,2,5)^{\top}$ , we obtain an equilibrium $x^{*}=(2,1,0.4)^{\top}$ . Then according to theorem 20, for any compact set $\mathbb{U}\subseteq\mathbb{R}^{4}_{>0}$ containing $u^{*}$ , the MAS (50), with input-value set $\mathbb{U}$ , is semiglobal ISS with respect to $(x^{*},u^{*})$ .

Also, we make numerical simulations for this system with two kinds of different inputs starting from the same initial point $x_{0}=(1.9,1.3,0.2)^{\top}$ , shown in fig. 3. As can be seen, for the first kind of input (converging to $u^{*}$ ), the state will converge to $x^{*}$ , while for the second kind of input (oscillating along $u^{*}$ ), the long-term behavior of the state is controlled by a function of $|u-u^{*}|$ .

References

[1] D. F. Anderson, G. Craciun, M. Gopalkrishnan, and C. Wiuf (2015) Lyapunov functions, stationary distributions, and non-equilibrium potential for reaction networks. Bulletin of Mathematical Biology 77 (9), pp. 1744–1767. External Links: Document Cited by: §1.
[2] D. F. Anderson, B. Joshi, and A. Deshpande (2021) On reaction network implementations of neural networks. Journal of the Royal Society Interface 18 (177). External Links: Document Cited by: §1, §5.
[3] D. F. Anderson and B. Joshi (2025) Chemical mass-action systems as analog computers: Implementing arithmetic computations at specified speed. Theoretical Computer Science 1025. External Links: Document Cited by: §1.
[4] D. Angeli, P. De Leenheer, and E. D. Sontag (2011) Persistence results for chemical reaction networks with time-dependent kinetics and no global conservation laws. SIAM Journal on Applied Mathematics 71 (1), pp. 128–146. External Links: Document Cited by: §1, §1, §2.3.
[5] B. Boros (2019) Existence of positive steady states for weakly reversible mass-action systems. SIAM Journal on Mathematical Analysis 51 (1), pp. 435–449. External Links: Document Cited by: §1.
[6] C. Chalk, N. Kornerup, W. Reeves, and D. Soloveichik (2019) Composable rate-independent computation in continuous chemical reaction networks. IEEE/ACM Transactions on Computational Biology and Bioinformatics 18 (1), pp. 250–260. External Links: Document Cited by: §1.
[7] M. Chaves and E. D. Sontag (2002) State-estimators for chemical reaction networks of Feinberg-Horn-Jackson zero deficiency type. European Journal of Control 8 (4), pp. 343–359. External Links: Document Cited by: §1, §2.3, §2.3, §3.1, §3.1, §3.2.
[8] M. Chaves (2005) Input-to-state stability of rate-controlled biochemical networks. SIAM Journal on Control and Optimization 44 (2), pp. 704–727. External Links: Document Cited by: §1, §2.3, §2.3, §2.3, §3.1, §3.2, Remark 17.
[9] Z. Chen, J. M. Linton, S. Xia, X. Fan, D. Yu, J. Wang, R. Zhu, and M. B. Elowitz (2024) A synthetic protein-level neural network in mammalian cells. Science 386 (6727), pp. 1243–1250. External Links: Document Cited by: §1.
[10] K. M. Cherry and L. Qian (2025) Supervised learning in DNA neural networks. Nature 645 (8081), pp. 639–647. External Links: Document Cited by: §1.
[11] G. Craciun, A. Deshpande, and J. Jin (2024) Weakly reversible deficiency one realizations of polynomial dynamical systems. Discrete and Continuous Dynamical Systems-B 29 (6). External Links: Document Cited by: §4.1.
[12] G. Craciun, J. Jin, and M. Sorea (2025) The structure of the toric locus of a reaction network. Nonlinearity 38 (1). External Links: Document Cited by: Appendix A, 1st item, §1, §2.3, §3.2, §3.2, Definition 15.
[13] G. Craciun, J. Jin, and P. Y. Yu (2020) An efficient characterization of complex-balanced, detailed-balanced, and weakly reversible systems. SIAM Journal on Applied Mathematics 80 (1), pp. 183–205. External Links: Document Cited by: §4.1.
[14] G. Craciun, J. Jin, and P. Y. Yu (2023) An algorithm for finding weakly reversible deficiency zero realizations of polynomial dynamical systems. SIAM Journal on Applied Mathematics 83 (4), pp. 1717–1737. External Links: Document Cited by: §4.1.
[15] G. Craciun, F. Nazarov, and C. Pantea (2013) Persistence and permanence of mass-action and power-law dynamical systems. SIAM Journal on Applied Mathematics 73 (1), pp. 305–329. External Links: Document Cited by: §5, Example 39.
[16] G. Craciun and C. Pantea (2008) Identifiability of chemical reaction networks. Journal of Mathematical Chemistry 44 (1), pp. 244–259. External Links: Document Cited by: §4.1.
[17] A. Deshpande (2023) Source-only realizations, weakly reversible deficiency one networks, and dynamical equivalence. SIAM Journal on Applied Dynamical Systems 22 (2), pp. 1502–1521. External Links: Document Cited by: §4.1.
[18] Y. Fan, X. Zhang, C. Gao, and D. Dochain (2025) Automatic implementation of neural networks through reaction networks-part I: circuit design and convergence analysis. IEEE Transactions on Automatic Control 70 (10), pp. 6356–6371. External Links: Document Cited by: §1, §5.
[19] Z. Fang and C. Gao (2019) Lyapunov function partial differential equations for chemical reaction networks: some special cases. SIAM Journal on Applied Dynamical Systems 18 (2), pp. 1163–1199. External Links: Document Cited by: §1.
[20] M. Feinberg (1987) Chemical reaction network structure and the stability of complex isothermal reactors—I. the deficiency zero and deficiency one theorems. Chemical Engineering Science 42 (10), pp. 2229–2268. External Links: Document Cited by: §1, §2.2, §2.3.
[21] M. Feinberg (2019) Foundations of chemical reaction network theory. Vol. 10, Springer. Cited by: §1, §2.1, §2.
[22] M. Gopalkrishnan, E. Miller, and A. Shiu (2014) A geometric approach to the global attractor conjecture. SIAM Journal on Applied Dynamical Systems 13 (2), pp. 758–797. External Links: Document Cited by: §5.
[23] L. Hoessly (2021) Stationary distributions via decomposition of stochastic reaction networks. Journal of Mathematical Biology 82 (7). External Links: Document Cited by: §1.
[24] H. Hong, J. Kim, M. Ali Al-Radhawi, E. D. Sontag, and J. K. Kim (2021) Derivation of stationary distributions of biochemical reaction networks via structure transformation. Communications Biology 4 (1). External Links: Document Cited by: Appendix B.
[25] F. Horn and R. Jackson (1972) General mass action kinetics. Archive for Rational Mechanics and Analysis 47 (2), pp. 81–116. External Links: Document Cited by: §1, §2.2, §2.3.
[26] R. Jiang, Y. Fan, D. Fan, C. Gao, and D. Dochain (2025) Input-to-state stability-based chemical reaction networks composition for molecular computations. arXiv preprint arXiv:2506.12056. Cited by: §2.3, §2.3.
[27] R. Jiang, C. Gao, and D. Dochain (2026) Structure-conditioned input-to-state stability for layer-by-layer molecular computations in parallel chemical reaction networks. Automatica 193. External Links: Document Cited by: §2.3, §2.3, §5, §5.
[28] M. D. Johnston, D. Siegel, and G. Szederkényi (2012) A linear programming approach to weak reversibility and linear conjugacy of chemical reaction networks. Journal of Mathematical Chemistry 50 (1), pp. 274–288. External Links: Document Cited by: §4.1.
[29] M. D. Johnston, D. Siegel, and G. Szederkényi (2013) Computing weakly reversible linearly conjugate chemical reaction networks with minimal deficiency. Mathematical Biosciences 241 (1), pp. 88–98. External Links: Document Cited by: §4.1.
[30] M. D. Johnston and D. Siegel (2011) Linear conjugacy of chemical reaction networks. Journal of Mathematical Chemistry 49 (7), pp. 1263–1282. External Links: Document Cited by: 2nd item, §2.3, §4.1, §4.1, §4, Definition 19, Example 23.
[31] M. Ke, Z. Fang, and C. Gao (2019) Complex balancing reconstructed to the asymptotic stability of mass-action chemical reaction networks with conservation laws. SIAM Journal on Applied Mathematics 79 (1), pp. 55–74. External Links: Document Cited by: Appendix A, 2nd item, §2.3, §4.2, §4.2, §4.2, §4, Lemma 29, Definition 30, Definition 33, Lemma 34, Example 37.
[32] P. M. Loriaux and A. Hoffmann (2013) A protein turnover signaling motif controls the stimulus-sensitivity of stress response pathways. PLoS Computational Biology 9 (2). External Links: Document Cited by: Example 18.
[33] C. D. Nicholls, K. G. McLure, M. A. Shields, and P. W. Lee (2002) Biogenesis of p53 involves cotranslational dimerization of monomers and posttranslational dimerization of dimers: implications on the dominant negative effect. Journal of Biological Chemistry 277 (15), pp. 12937–12945. External Links: Document Cited by: Example 18.
[34] E. D. Sontag (1989) Smooth stabilization implies coprime factorization. IEEE Transactions on Automatic Control 34 (4), pp. 435–443. External Links: Document Cited by: §1.
[35] G. Szederkényi, K. Hangos, and Z. Tuza (2012) Finding weakly reversible realizations of chemical reaction networks using optimization. MATCH Communications in Mathematical and in Computer Chemistry 67 (1), pp. 193–212. Cited by: §4.1.
[36] M. A. Vághy and G. Szederkényi (2023) Persistence and stability of generalized ribosome flow models with time-varying transition rates. Plos One 18 (7). External Links: Document Cited by: §2.3.
[37] M. Vasić, C. Chalk, A. Luchsinger, S. Khurshid, and D. Soloveichik (2022) Programming and training rate-independent chemical reaction networks. Proceedings of the National Academy of Sciences 119 (24). External Links: Document Cited by: §1.