Robust and structural ergodicity analysis of stochastic biomolecular networks involving synthetic antithetic integral controllers

TL;DR

This paper develops exact algebraic and polynomial methods to verify ergodicity and output controllability in stochastic biomolecular networks with synthetic antithetic integral controllers, accounting for uncertainties.

Contribution

It introduces robust, exact theoretical techniques for analyzing ergodicity and controllability in stochastic reaction networks with uncertainties, extending previous linear programming approaches.

Findings

Exact algebraic methods verify ergodicity under uncertainty.

Polynomial techniques confirm output controllability.

Illustrative examples demonstrate the methods' effectiveness.

Abstract

Ergodicity and output controllability have been shown to be fundamental concepts for the analysis and synthetic design of closed-loop stochastic reaction networks, as exemplified by the use of antithetic integral feedback controllers. In [Gupta, Briat & Khammash, PLoS Comput. Biol., 2014], some ergodicity and output controllability conditions for unimolecular and certain classes of bimolecular reaction networks were obtained and formulated through linear programs. To account for context dependence, these conditions were later extended in [Briat & Khammash, CDC, 2016] to reaction networks with uncertain rate parameters using simple and tractable, yet potentially conservative, methods. Here we develop some exact theoretical methods for verifying, in a robust setting, the original ergodicity and output controllability conditions based on algebraic and polynomial techniques. Some examples are given for illustration.

Equations

$${\mathcal{R}}_k:\sum _{i=1}^d{\zeta }_{k,i}^l{\mathbf{X}}_{\mathbf{i}}\stackrel{{\rho }_k}{\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\to }\sum _{i=1}^d{\zeta }_{k,i}^r{\mathbf{X}}_{\mathbf{i}},k=1,\dots ,K$$

$$\lambda (x)=\left[\begin{array}{c}{w}_0({\rho }_0)\\ W({\rho }_u)x\\ Y({\rho }_b,x)\end{array}\right]$$

$$A({\rho }_u):=S_{u}W({\rho }_u)\text{and}{b}_0({\rho }_0):=S_0{w}_0({\rho }_0).$$

$$\varnothing \stackrel{\mu }{\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\to }{\mathbf{Z}}_1,\varnothing \stackrel{\theta X_{\ell }}{\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\to }{\mathbf{Z}}_2,{\mathbf{Z}}_1+{\mathbf{Z}}_2\stackrel{\eta }{\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\to }\varnothing ,\varnothing \stackrel{k{Z}_1}{\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\to }{\mathbf{X}}_1$$

$$v^{T}A<0,w^{T}A+e_{\ell }^T=0\text{and}\frac{\mu }{\theta }>\frac{v^T{b}_0}{c{e}_{\ell }^{T}v}$$

$$S_u=\left[\begin{array}{ccc}{S}_{dg}& S_{ct}& S_{cv}\end{array}\right]$$

$${\rho }_{\bullet }\in {\mathcal{P}}_{\bullet }:=[{\rho }_{\bullet }^{-},{\rho }_{\bullet }^{+}],0\le {\rho }_{\bullet }^{-}\le {\rho }_{\bullet }^{+}<\infty$$

$$A^{+}({\rho }_{cv}):=A({\rho }_{dg}^{-},{\rho }_{ct}^{+},{\rho }_{cv})$$

$$A({\rho }_{dg},{\rho }_{ct},{\rho }_{cv})\le A^{+}({\rho }_{cv})=A({\rho }_{dg}^{-},{\rho }_{ct}^{+},{\rho }_{cv})$$

$$A^{+}({\rho }_{cv}):=A({\rho }_{dg}^{-},{\rho }_{ct}^{+},{\rho }_{cv})$$

$$\begin{array}[]{rcl}v(\rho_{cv})^{T}&=&-(-1)^{d}\mathds{1}_{d}^{T}\det(A^{+}(\rho_{cv}))A^{+}(\rho_{cv})^{-1}\\ &=&(-1)^{d+1}\mathds{1}^{T}_{d}\textnormal{Adj}(A^{+}(\rho_{cv}))>0\end{array}$$

$$v({\rho }_u)>0,v({\rho }_u{)}^T{S}_b=0\text{and}v({\rho }_u{)}^{T}A({\rho }_u)<0$$

$$\tilde{v}({\rho }_{cv}{)}^T{S}_b^{\perp }>0\text{and}\tilde{v}({\rho }_{cv}{)}^T{S}_b^{\perp }{A}^{+}({\rho }_{cv})<0$$

$$\begin{array}[]{rcl}v(\rho_{cv})^{T}A(\rho_{u})&=&v(\rho_{cv})^{T}(A^{+}(\rho_{cv})-\Delta(\rho_{dg},\rho_{ct}))\\ &\leq&v(\rho_{cv})^{T}A^{+}(\rho_{cv})<0\end{array}$$

$$A({\rho }_{dg},{\rho }_{cv},0)\le A(\theta \text{mathds}1,{\rho }_{cv},0)$$

$$\det (A({\rho }_u))=\det (A_{\rho })\det (I-D({\rho }_{ct})W_{ct}{A}_{\rho }^{-1}{S}_{ct})$$

$$\det (I-D({\rho }_{ct})M)>0$$

$$A^{22}({\rho }_{dg},{\rho }_{cv}):={\left[\begin{array}{c}0\\ I\end{array}\right]}^{T}A({\rho }_{dg},{\rho }_{cv},0)\left[\begin{array}{c}0\\ I\end{array}\right]$$

$$A=\left[\begin{array}{ccc}-{\gamma }_s& 0& k_{rs}\\ 0& -({\gamma }_i+k_{ir})& 0\\ 0& k_{ir}& -({\gamma }_r+k_{rs})\end{array}\right],S_b=\left[\begin{array}{c}-1\\ 1\\ 0\end{array}\right]$$

$${\tilde{v}}^T{S}_b^{\perp }A<0\iff {\tilde{v}}^T\left[\begin{array}{cc}-({\gamma }_i+k_{ir})& k_{rs}\\ k_{ir}& -({\gamma }_r+k_{rs})\end{array}\right]<0.$$

$${\tilde{v}}^T\left[\begin{array}{cc}-({\gamma }_i+k_{ir})& k_{rs}\\ k_{ir}& -({\gamma }_r+k_{rs})\end{array}\right]=\left[\begin{array}{cc}-{\gamma }_i& -{\gamma }_r\end{array}\right]<0$$

$$A=\left[\begin{array}{ccccc}-{\delta }_{M_A}& 0& 0& 0& 0\\ {\beta }_A& -{\delta }_A& 0& 0& 0\\ 0& 0& -{\delta }_{M_R}& 0& 0\\ 0& 0& {\beta }_R& -{\delta }_R& {\delta }_A\\ 0& 0& 0& 0& -{\delta }_A\end{array}\right],S_b=\left[\begin{array}{c}0\\ -1\\ 0\\ -1\\ 1\end{array}\right]$$

$${\tilde{v}}^T{S}_b^{\perp }A<0\iff {\tilde{v}}^T\left[\begin{array}{cccc}-{\delta }_{M_A}& 0& 0& 0\\ {\beta }_A& -{\delta }_A& 0& 0\\ 0& 0& -{\delta }_{M_R}& 0\\ 0& 0& {\beta }_R& -{\delta }_R\end{array}\right]<0$$

$$W_{ct}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]\text{and}{S}_{ct}={\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right]}^T$$

$$A=\left[\begin{array}{ccc}-({\gamma }_1+{\alpha }_1)& 0& k_1\\ k_2& -({\gamma }_2+{\alpha }_2)& 0\\ 0& k_3& -k_1\end{array}\right].$$

$$A_{\text{mathds}1}=\left[\begin{array}{ccc}-2& 0& 1\\ 1& -2& 0\\ 0& 1& -1\end{array}\right]$$

$$A_{\text{mathds}1}=\left[\begin{array}{ccc}-1& 0& 1\\ 0& -1& 0\\ 0& 0& -1\end{array}\right]\text{and}{A}_{\text{mathds}1}^{-1}=\left[\begin{array}{ccc}-1& 0& -1\\ 0& -1& 0\\ 0& 0& -1\end{array}\right]$$

$$W_{ct}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]\text{and}{S}_{ct}={\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\end{array}\right]}^T$$

$$W_{ct}{A}_{\text{mathds}1}^{-1}{S}_{ct}=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$$

$$A^{+}(k_1)=\left[\begin{array}{ccc}-{\gamma }_1^{-}& 0& k_1\\ k_2^{+}& -{\gamma }_2^{-}& 0\\ 0& k_3^{+}& -k_1\end{array}\right].$$

Full text

1 Introduction

The main objective of synthetic biology is the rational and systematic design of biological networks that can achieve de-novo functions such as the heterologous production of a metabolite of interest [1]. Besides the obvious necessity of developing experimental methodologies allowing for the reliable implementation of synthetic networks, tailored theoretical and computational tools for their design, their analysis and their simulation also need to developed. Indeed, theoretical tools that could predict certain properties (e.g. a stable/oscillatory/switching behavior, controllable trajectories, etc.) of a synthetic biological network from an associated model formulated, for instance, in terms of a reaction network [2, 3, 4], could pave the way to the development of iterative procedures for the systematic design of efficient synthetic biological networks. Such an approach would allow for a faster design procedure than those involving fastidious experimental steps, and would give insights on how to adapt the current design in order to improve a certain design criterion. This way, synthetic biology would become conceptually much closer to existing theoretically-driven engineering disciplines, such as control engineering.

However, while such methods are well-developed for deterministic models (i.e. deterministic reaction networks), they still lag behind in the stochastic setting. This lack of tools is quite problematic since it is now well-known that stochastic reaction networks [4] are versatile modeling tools that can capture the inherent stochastic behavior of living cells [5, 6] and can exhibit several interesting properties that are absent for their deterministic counterparts [7, 8, 9, 10]. Under the well-mixed assumption, it is known [11, 12] that such random dynamics can be well represented as a continuous-time jump Markov process evolving on the $d$-dimensional nonnegative integer lattice where $d$ is the number of distinct molecular species involved in the network. Sufficient conditions for checking the ergodicity of open unimolecular and bimolecular stochastic reaction networks has have been proposed in [13] and formulated in terms of linear programs. The concept of ergodicity is of fundamental importance as it can serve as a basis for the development of a control theory for biological systems. Indeed, verifying the ergodicity of a control system, consisting for instance of an endogenous biomolecular network controlled by a synthetic controller (see Fig. 1A), would prove that the closed-loop system is well-behaved (e.g. ergodic with bounded first- and second-order moments) and that the designed control system achieves its goal (e.g. set-point tracking and perfect adaptation properties). This procedure is analogous to that of checking the stability of a closed-loop system in the deterministic setting; see e.g. [14]. Additionally, designing synthetic circuits achieving a given function that are provably ergodic could allow for the rational design of synthetic networks that can exploit noise in their function. A recent example is that of the antithetic integral feedback controller proposed in [9] (see also Fig. 1B) that has been shown to induce an ergodic closed-loop network when some conditions on the endogenous network to be controlled are met.

A major limitation of the ergodicity conditions obtained in [13, 9] is that they only apply to networks with fixed and known rate parameters – an assumption that is rarely met in practice as the rate parameters are usually poorly known and context dependent. This has motivated the consideration of networks with uncertain rate parameters in [15] using two different approaches. The first one is quantitative and considers networks having an interval matrix as characteristic matrix [16]. It was notably shown that checking the ergodicity and the output controllability of those networks reduces to checking the Hurwitz stability of a single matrix and the output controllability of a single positive linear system. It was also shown that these conditions exactly write as a simple linear program having the same complexity as the program associated with the nominal case; i.e. in the case of constant and fixed rate parameters. The second approach is qualitative and is based on the theory of sign-matrices [17, 18] which has been extensively studied and considered for the qualitative analysis of dynamical systems. Sign-matrices have also been considered in the context of reaction networks, albeit much more sporadically; see e.g. [19, 20, 21, 15, 22]. In this case, again, the conditions obtained in [21, 15] can be stated as a very simple linear program that can be shown to be equivalent to some graph theoretical conditions.

However, these approaches can be very conservative when the entries of the characteristic matrix of the network are not independent – a situation that appears when conversion reactions are involved in the networks. In order to solve this problem, the actual parameter dependence needs to be exactly captured and many approaches exist to attack this problem such as, to cite a few, $\mu$-analysis [23, 24], small-gain methods [24, 25], eigenvalue and perturbation methods [26, 27], interval matrices [16, 15], sign-matrices [19, 20, 21, 15, 22], Lyapunov methods [28, 29, 30, 31], etc.

As some of the methods just cited above can be conservative or may yield too complex conditions, we propose to develop an approach that is tailored to our problem by exploiting its inherent properties. Several conditions for the robust ergodicity of unimolecular and biomolecular networks are first obtained in terms of a sign switching property for the determinant of the upper-bound of the characteristic matrix of the network. This condition also alternatively formulates as the existence of a positive vector depending polynomially on the uncertain parameters and satisfying certain inequality conditions. The complexity of the problem is notably reduced by exploiting the Metzler structure111A matrix is Metzler if its off-diagonal elements are nonnegative. of the matrices involved and through the use of various algebraic results such as the Perron-Frobenius theorem. The structural ergodicity of unimolecular networks is also considered and shown to reduce to the analysis of constant matrices when some certain realistic assumptions are met. It is notably shown in the examples that this latter result can be applied to bimolecular networks in some situations. The examples also illustrate that the proposed approach can be used to establish the robust or structural ergodicity of reaction networks for which the methods proposed in [15] fail.

Outline. Preliminaries on reaction networks, ergodicity analysis and antithetic integral control are given in Section 2. Section 3 is devoted to the robust ergodicity analysis of unimolecular and bimolecular reactions networks while the problem of establishing the structural ergodicity of unimolecular reaction networks is addressed in Section 4. Examples are finally treated in Section 5.

Notations. The standard basis for $\mathbb{R}^{d}$ is denoted by $\{e_{i}\}_{i=1}^{d}$. The sets of integers, nonnegative integers, nonnegative real numbers and positive real numbers are denoted by $\mathbb{Z}$, $\mathbb{Z}_{\geq 0}$, $\mathbb{R}_{\geq 0}$ and $\mathbb{R}_{>0}$, respectively. The $d$-dimensional vector of ones is denoted by $\mathds{1}_{d}$ (the index will be dropped when the dimension is obvious). For vectors and matrices, the inequality signs $\leq$ and $<$ act componentwise. Finally, the vector or matrix obtained by stacking the elements $x_{1},\ldots,x_{d}$ is denoted by $\operatorname*{col}_{i=1}^{d}(x_{i})$ or $\operatorname*{col}(x_{1},\ldots,x_{d})$. The diagonal operator $\operatorname*{diag}(\cdot)$ is defined analogously. The spectral radius of a matrix $M\in\mathbb{R}^{n\times n}$ is defined as $\varrho(M)=\max\{|\lambda|:\det(\lambda I-M)=0\}$.

2 Preliminaries on reaction networks

2.1 Reaction networks

We consider here a reaction network with $d$ molecular species $\boldsymbol{X_{1}},\ldots,\boldsymbol{X_{d}}$ that interacts through $K$ reaction channels $\mathcal{R}_{1},\ldots,\mathcal{R}_{K}$ defined as

[TABLE]

where $\rho_{k}\in\mathbb{R}_{>0}$ is the reaction rate parameter and $\zeta_{k,i}^{l},\zeta_{k,i}^{r}\in\mathbb{Z}_{\geq 0}$. Each reaction is additionally described by a stoichiometric vector and a propensity function. The stoichiometric vector of reaction $\mathcal{R}_{k}$ is given by $\zeta_{k}:=\zeta_{k}^{r}-\zeta_{k}^{l}\in\mathbb{Z}^{d}$ where $\zeta_{k}^{r}=\operatorname*{col}(\zeta_{k,1}^{r},\ldots,\zeta_{k,d}^{r})$ and $\zeta_{k}^{l}=\operatorname*{col}(\zeta_{k,1}^{l},\ldots,\zeta_{k,d}^{l})$. In this regard, when the reaction $\mathcal{R}_{k}$ fires, the state jumps from $x$ to $x+\zeta_{k}$. We define the stoichiometry matrix $S\in\mathbb{Z}^{d\times K}$ as $S:=\begin{bmatrix}\zeta_{1}\ldots\zeta_{K}\end{bmatrix}$. When the kinetics is mass-action, the propensity function of reaction $\mathcal{R}_{k}$ is given by $\textstyle\lambda_{k}(x)=\rho_{k}\prod_{i=1}^{d}\frac{x_{i}!}{(x_{i}-\zeta_{k,i}^{l})!}$ and is such that $\lambda_{k}(x)=0$ if $x\in\mathbb{Z}_{\geq 0}^{d}$ and $x+\zeta_{k}\notin\mathbb{Z}_{\geq 0}^{d}$. We denote this reaction network by $(\boldsymbol{X},\mathcal{R})$. Under the well-mixed assumption, this network can be described by a continuous-time Markov process $(X_{1}(t),\ldots,X_{d}(t))_{t\geq 0}$ with state-space $\mathbb{Z}_{\geq 0}^{d}$; see e.g. [11].

2.2 Ergodicity of unimolecular and bimolecular reaction networks

Let us assume here that the network $(\boldsymbol{X},\mathcal{R})$ is at most bimolecular and that the reaction rates are all independent of each other. In such a case, the propensity functions are polynomials of at most degree 2 and we can write the propensity vector as

[TABLE]

where $w_{0}(\rho_{0})\in\mathbb{R}^{n_{0}}_{\geq 0}$, $W(\rho_{u})x\in\mathbb{R}^{n_{u}}_{\geq 0}$ and $Y(\rho_{b},x)\in\mathbb{R}^{n_{b}}_{\geq 0}$ are the propensity vectors associated the zeroth-, first- and second-order reactions, respectively. Their respective rate parameters are also given by $\rho_{0}$, $\rho_{u}$ and $\rho_{b}$, and according to this structure, the stoichiometric matrix is decomposed as $S=:\begin{bmatrix}S_{0}&S_{u}&S_{b}\end{bmatrix}$. Before stating the main results of the section, we need to introduce the following terminology:

Definition 1

The characteristic matrix $A(\rho_{u})$ and the offset vector $b_{0}(\rho)$ of a bimolecular reaction network $(\boldsymbol{X},\mathcal{R})$ are defined as

[TABLE]

A particularity is that the matrix $A(\rho_{u})$ is Metzler (i.e. all the off-diagonal elements are nonnegative) for all $\rho_{u}\geq 0$. This property plays an essential role in the derivation of the results of [9] and will also be essential for the derivation of the main results of this paper. It is also important to define the property of ergodicity:

Definition 2** ([32])**

The Markov process associated with the reaction network $(\boldsymbol{X},\mathcal{R})$ is said to be ergodic if its probability distribution globally converges to a unique stationary distribution. It is exponentially ergodic if the convergence to the unique stationary distribution is exponential.

We then have the following result:

Theorem 3** ([13])**

Let us consider an irreducible222Computationally tractable conditions for checking the irreducibility of reaction networks are provided in [33]. bimolecular reaction network $(\boldsymbol{X},\mathcal{R})$ with fixed rate parameters; i.e. $A=A(\rho_{u})$ and $b_{0}=b_{0}(\rho_{0})$. Assume that there exists a vector $v\in\mathbb{R}^{d}_{>0}$ such that $v^{T}S_{b}=0$ and $v^{T}A<0$. Then, the reaction network $(\boldsymbol{X},\mathcal{R})$ is exponentially ergodic and all the moments are bounded and converging.

We also have the following immediate corollary pertaining on unimolecular reaction networks:

Corollary 4

Let us consider an irreducible unimolecular reaction network $(\boldsymbol{X},\mathcal{R})$ with fixed rate parameters; i.e. $A=A(\rho_{u})$ and $b_{0}=b_{0}(\rho_{0})$. Assume that there exists a vector $v\in\mathbb{R}^{d}_{>0}$ such that $v^{T}A<0$. Then, the reaction network $(\boldsymbol{X},\mathcal{R})$ is exponentially ergodic and all the moments are bounded and converging.

2.3 Antithetic integral control of unimolecular networks

Antithetic integral control has been first proposed in [9] for solving the perfect adaptation problem in stochastic reaction networks. The underlying idea is to augment the open-loop network $(\boldsymbol{X},\mathcal{R})$ with an additional set of species and reactions (the controller). The usual set-up is that this controller network acts on the production rate of the molecular species $\boldsymbol{X_{1}}$ (the actuated species) in order to steer the mean value of the controlled species $\boldsymbol{X_{\ell}}$, $\ell\in\{1,\ldots,d\}$, to a desired set-point (the reference). To the regulation problem, it is often sought to have a controller that can ensure perfect adaptation for the controlled species. As proved in [9], the antithetic integral control motif $(\boldsymbol{Z},\mathcal{R}^{c})$ defined with

[TABLE]

solves this control problem with the set-point being equal to $\mu/\theta$. Above, $\boldsymbol{Z_{1}}$ and $\boldsymbol{Z_{2}}$ are the controller species. The four controller parameters $\mu,\theta,\eta,k>0$ are assumed to be freely assignable to any desired value. The first reaction is the reference reaction as it encodes part of the reference value $\mu/\theta$ as its own rate. The second one is the measurement reaction that produces the species $\boldsymbol{Z_{2}}$ at a rate proportional to the current population of the controlled species $\boldsymbol{X_{\ell}}$. The third reaction is the comparison reaction as it compares the populations of the controller species and annihilates one molecule of each when these populations are both positive. Finally, the fourth reaction is the actuation reaction that produces the actuated species $\boldsymbol{X_{1}}$ at a rate proportional to the controller species $\boldsymbol{Z_{1}}$.

The following fundamental result states conditions under which a unimolecular reaction network can be controlled using an antithetic integral controller:

Theorem 5** ([9])**

Suppose that the open-loop reaction network $(\boldsymbol{X},\mathcal{R})$ is unimolecular and that the state-space of the closed-loop reaction network $((\boldsymbol{X},\boldsymbol{Z}),\mathcal{R}\cup\mathcal{R}^{c})$ is irreducible. Let us assume that $\rho_{0}$ and $\rho_{u}$ are fixed and known (i.e. $A=A(\rho_{u})$ and $b_{0}=b_{0}(\rho_{0})$) and assume, further, that there exist vectors $v\in\mathbb{R}_{>0}^{d}$, $w\in\mathbb{R}_{\geq 0}^{d}$, $w_{1}>0$, such that

[TABLE]

where $c>0$ verifies $v^{T}(A+cI)\leq 0$.

Then, for any values for the controller rate parameters $\eta,k>0$, (i) the closed-loop network is ergodic, (ii) $\mathbb{E}[X_{\ell}(t)]\to\mu/\theta$ as $t\to\infty$ and (iii) $\mathbb{E}[X(t)X(t)^{T}]$ is bounded over time.

We can see that the conditions above consist of the combination of an ergodicity condition (i.e. $v^{T}A<0$) and an output controllability condition for Hurwitz stable matrices $A$ (i.e. $w^{T}A+e_{\ell}^{T}=0$ with $w_{1}>0$), which are fully consistent with the considered control problem. Note, however, that unlike in the deterministic case, the above result proves that the closed-loop network cannot be unstable if the conditions on the open-loop network are met; i.e. have trajectories that grow unboundedly with time. This is illustrated in more details in the supplemental material of [9].

As the robust/structural output controllability problem has been completely solved in [15], we will only focus on checking the ergodicity condition in the rest of the paper.

3 Robust ergodicity of reaction networks

3.1 Preliminaries

The following lemma will be useful in proving the main results of this section:

Lemma 6

Let us consider a parameter-dependent Metzler matrix $M(\theta)\in\mathbb{R}^{d\times d}$, $\theta\in\Theta\subset\mathbb{R}^{N}_{\geq 0}$, where $\Theta$ is compact and connected. Then, the following statements are equivalent:

1. (a)

The matrix $M(\theta)$ is Hurwitz stable for all $\theta\in\Theta$. 2. (b)

The coefficients of the characteristic polynomial of $M(\theta)$ are positive $\Theta$. 3. (c)

The following conditions hold:

1. (c1)

there exists a $\theta^{*}\in\Theta$ such that $M(\theta^{*})$ is Hurwitz stable, and 2. (c2)

for all $\theta\in\Theta$ we have that $(-1)^{d}\det(M(\theta))>0$.

*Proof : * The proof of the equivalence between (a) and (b) follows, for instance, from [34] and is omitted. It is also immediate to prove that (b) implies (c) since if $M(\theta)$ is Hurwitz stable for all $\theta\in\Theta$ then (c1) holds and the constant term of the characteristic polynomial of $M(\theta)$ is positive on $\theta\in\Theta$. Using now the fact that that constant term is equal to $(-1)^{d}\det(M(\theta))$ yields the result.

To prove that (c) implies (a), we use the contraposition. Hence, let us assume that there exists at least a $\theta_{u}\in\Theta$ for which the matrix $M(\theta_{u})$ is not Hurwitz stable. If such a $\theta_{u}$ can be arbitrarily chosen in $\Theta$, then this implies the negation of statement (c1) (i.e. for all $\theta^{*}\in\Theta$ the matrix $M(\theta^{*})$ is not Hurwitz stable) and the first part of the implication is proved.

Let us consider now the case where there exists some $\theta_{s}\in\Theta$ such that $M(\theta_{s})$ is Hurwitz stable. Let us then choose a $\theta_{u}$ and a $\theta_{s}$ such that $M(\theta_{u})$ is not Hurwitz stable and $M(\theta_{s})$ is. Since $\Theta$ is connected, then there exists a path $\mathscr{P}\subset\Theta$ from $\theta_{s}$ and $\theta_{u}$. From Perron-Frobenius theorem, the dominant eigenvalue, denoted by $\lambda_{PF}(\cdot)$, is real and hence, we have that $\lambda_{PF}(M(\theta_{s}))<0$ and $\lambda_{PF}(M(\theta_{u}))\geq 0$. Hence, from the continuity of eigenvalues then there exists a $\theta_{c}\in\mathscr{P}$ such that $\lambda_{PF}(M(\theta_{c}))=0$, which then implies that $(-1)^{d}\det(M(\theta_{c}))=0$, or equivalently, that the negation of (c2) holds. This concludes the proof. $\diamondsuit$

Before stating the next main result of this section, let us assume that $S_{u}$ in Definition 1 has the following form

[TABLE]

where $S_{dg}\in\mathbb{R}^{d\times n_{dg}}$ is a matrix with nonpositive columns, $S_{ct}\in\mathbb{R}^{d\times n_{ct}}$ is a matrix with nonnegative columns and $S_{cv}\in\mathbb{R}^{d\times n_{cv}}$ is a matrix with columns containing exactly one negative entry and at least one positive entry. Also, decompose accordingly $\rho_{u}$ as $\rho_{u}=:\operatorname*{col}(\rho_{dg},\rho_{ct},\rho_{cv}\}$ and define

[TABLE]

where $\bullet\in\{dg,ct,cv\}$ and let $\mathcal{P}_{u}:=\mathcal{P}_{dg}\times\mathcal{P}_{ct}\times\mathcal{P}_{cv}$.

In this regard, we can alternatively rewrite the matrix $A(\rho_{u})$ as $A(\rho_{dg},\rho_{ct},\rho_{cv})$. We then have the following result:

Lemma 7

The following statements are equivalent:

1. (a)

The matrix $A(\rho_{u})$ is Hurwitz stable for all $\rho_{u}\in\mathcal{P}_{u}$. 2. (b)

The matrix

[TABLE]

is Hurwitz stable for all $\rho_{cv}\in\mathcal{P}_{cv}$.

*Proof : * The proof that (a) implies (b) is immediate. To prove that (b) implies (a), first note that we have

[TABLE]

since for all $(\rho_{dg},\rho_{ct},\rho_{cv})\in\mathcal{P}_{u}$. Using the fact that for two Metzler matrices $B_{1},B_{2}$, the inequality $B_{1}\leq B_{2}$ implies $\lambda_{PF}(B_{1})\leq\lambda_{PF}(B_{2})$ [35], then we can conclude that $A(\rho_{dg}^{-},\rho_{ct}^{+},\rho_{cv})$ is Hurwitz stable for all $\rho_{cv}\in\mathcal{P}_{cv}$ if and only if the matrix $A(\rho_{dg},\rho_{ct},\rho_{cv})$ is Hurwitz stable for all $(\rho_{dg},\rho_{ct},\rho_{cv})\in\mathcal{P}_{u}$. This completes the proof. $\diamondsuit$

3.2 Unimolecular networks

The following theorem states the main result on the robust ergodicity of unimolecular reaction networks:

Theorem 8

Let $A(\rho_{u})\in\mathbb{R}^{d\times d}$ be the characteristic matrix of some unimolecular network and $\rho_{u}\in\mathcal{P}_{u}$. Then, the following statements are equivalent:

1. (a)

The matrix $A(\rho_{u})$ is Hurwitz stable for all $\rho_{u}\in\mathcal{P}_{u}$. 2. (b)

The matrix

[TABLE]

is Hurwitz stable for all $\rho_{cv}\in\mathcal{P}_{cv}$. 3. (c)

There exists a $\rho_{cv}^{s}\in\mathcal{P}_{cv}$ such that the matrix $A^{+}(\rho_{cv}^{s})$ is Hurwitz stable and the polynomial

$(-1)^{d}\det(A^{+}(\rho_{cv}))$* is positive for all $\rho_{cv}\in\mathcal{P}_{cv}$.* 4. (d)

There exists a polynomial vector-valued function ${v:\mathcal{P}_{cv}\mapsto\mathbb{R}_{>0}^{d}}$ of degree at most $d-1$ such that

${v(\rho_{cv})^{T}A^{+}(\rho_{cv})<0}$* for all $\rho_{cv}\in\mathcal{P}_{cv}$.*

*Proof : * The equivalence between the statement (a), (b) and (c) directly follows from Lemma 6 and Lemma 7. To prove the equivalence between the statements (b) and (d), first remark that (b) is equivalent to the fact that for any $q(\rho_{cv})>0$ on $\mathcal{P}_{cv}$, there exists a unique parameterized vector $v(\rho_{cv})\in\mathbb{R}^{d}$ such that $v(\rho_{cv})>0$ and $v(\rho_{cv})^{T}A^{+}(\rho_{cv})=-q(\rho_{cv})^{T}$ for all $\rho_{cv}\in\mathcal{P}_{cv}$. Choosing $q(\rho_{cv})=(-1)^{d}\mathds{1}_{n}\det(A^{+}(\rho_{cv}))$, we get that such a $v(\rho_{cv})$ is given by

[TABLE]

for all $\rho_{cv}\in\mathcal{P}_{cv}$. Since the matrix $A^{+}(\rho_{cv})$ is affine in $\rho_{cv}$, then the adjugate matrix $\textnormal{Adj}(A^{+}(\rho_{cv})$ contains entries of at most degree $d-1$ and the conclusion follows. $\diamondsuit$

Checking the condition (c) amounts to solving two problems. The first one is is concerned with the construction of a stabilizer $\rho_{cv}\in\mathcal{P}_{cv}$ for the matrix $A^{+}(\rho_{cv})$ whereas the second one is about checking whether a polynomial is positive on a compact set. The first problem can be easily solved by checking whether $A^{+}(\rho_{cv})$ is Hurwitz stable for some randomly chosen point in $\mathcal{P}_{cv}$. For the second one, optimization-based methods can be used such as those based on the Handelman’s Theorem combined with linear programming [36, 37] or Putinar’s Positivstellensatz combined with semidefinite programming [38, 39]. Note also that the degree $d-1$ is a worst case degree and that, in fact, polynomials of lower degree will in general be enough for proving the Hurwitz stability of the matrix $A^{+}(\rho_{cv})$ for all $\rho_{cv}\in\mathcal{P}_{cv}$. For instance, the matrices $A(\rho)$ and $A^{+}(\rho_{cv})$ are very sparse in general due to the particular structure of biochemical reaction networks. The sparsity property is not considered here but could be exploited to refine the necessary degree for the polynomial vector $v(\rho_{cv})$.

In is important to stress here that Theorem 8 can only be considered when the rate parameters are time-invariant (i.e. constant deterministic or drawn from a distribution). When they are time-varying (e.g. time-varying stationary random variables), a possible workaround relies on the use of a constant vector $v$ as formulated below:

Proposition 9** (Constant $v$)**

Let $A(\rho_{u})\in\mathbb{R}^{d\times d}$ be the characteristic matrix of some unimolecular network and $\rho_{u}\in\mathcal{P}_{u}$. Then, the following statements are equivalent:

1. (a)

There exists a vector $v\in\mathbb{R}^{d}_{>0}$ such that ${v^{T}A^{+}(\rho_{cv})<0}$ holds for all $\rho_{cv}\in\mathcal{P}_{cv}$. 2. (b)

There exists a vector $v\in\mathbb{R}^{d}_{>0}$ such that $v^{T}A^{+}(\theta)<0$ holds for all $\theta\in\textnormal{vert}(\mathcal{P}_{cv})$ where $\textnormal{vert}(\mathcal{P}_{cv})$ denotes the set of vertices of the set $\mathcal{P}_{cv}$.

*Proof : * The proof exploits the affine, hence convex, structure of the matrix $A^{+}(\rho_{cv})$. Using this property, it is indeed immediate to show that the inequality ${v^{T}A^{+}(\rho_{cv})<0}$ holds for all $\rho_{cv}\in\mathcal{P}_{cv}$ if and only if $v^{T}A^{+}(\theta)<0$ holds for all $\theta\in\textnormal{vert}(\mathcal{P}_{cv})$ (see e.g. [28] for a similar arguments in the context of quadratic Lyapunov functions). $\diamondsuit$

The above result is connected to the existence of a linear copositive Lyapunov function for a linear positive switched system with matrices in the family $\{A(\theta):\theta\in\textnormal{vert}(\mathcal{P}_{cv})\}$ for which many characterizations exist; see e.g. [40, 41].

3.3 Bimolecular networks

In the case of bimolecular networks, we have the following result:

Proposition 10

Let $A(\rho_{u})\in\mathbb{R}^{d\times d}$ be the characteristic matrix of some bimolecular network and $\rho_{u}\in\mathcal{P}_{u}$. Then, the following statements are equivalent:

1. (a)

There exists a polynomial vector-valued function ${v:\mathcal{P}_{u}\mapsto\mathbb{R}_{>0}^{d}}$ such that

[TABLE]

*for all $\rho_{u}\in\mathcal{P}_{u}$. * 2. (b)

There exists a polynomial vector-valued function ${\tilde{v}:\mathcal{P}_{cv}\mapsto\mathbb{R}^{d-n_{b}}}$ such that

[TABLE]

for all $\rho_{cv}\in\mathcal{P}_{cv}$ and where $n_{b}:=\operatorname*{rank}(S_{b})$ and $S_{b}^{\bot}S_{b}=0$, $S_{b}^{\bot}$ full-rank.

*Proof : * It is immediate to see that (a) implies (b). To prove the converse, first note that we have that $v(\rho_{cv})=(S_{b}^{\bot})^{T}\tilde{v}(\rho_{cv})$ verifies $v(\rho_{cv})^{T}S_{b}=0$ and $v(\rho_{cv})>0$ for all $\rho_{cv}\in\mathcal{P}_{cv}$. This proves the equality and the first inequality in (11). Observe now that for any $\rho_{u}\in\mathcal{P}_{u}$, there exists a nonnegative matrix $\Delta(\rho_{dg},\rho_{ct})\in\mathbb{R}^{d\times d}_{\geq 0}$ such that $A(\rho_{u})=A^{+}(\rho_{cv})-\Delta(\rho_{dg},\rho_{ct})$. Hence, we have that

[TABLE]

which proves the result. $\diamondsuit$

As in the unimolecular case, we have been able to reduce the number of parameters by using an upper-bound on the characteristic matrix. It is also interesting to note that the condition $\tilde{v}(\rho_{cv})^{T}S_{b}^{\bot}A^{+}(\rho_{cv})<0$ can be sometimes brought back to a problem of the form $\tilde{v}(\rho_{cv})^{T}M(\rho_{cv})<0$ for some square, and often Metzler, matrix $M(\rho_{cv})$ which can be dealt in the same way as in the unimolecular case.

The following result can be used when the parameters are time-varying and is the bimolecular analogue of Proposition 9:

Proposition 11** (Constant $v$)**

Let $A(\rho_{u})\in\mathbb{R}^{d\times d}$ be the characteristic matrix of some bimolecular network and $\rho_{u}\in\mathcal{P}_{u}$. Then, the following statements are equivalent:

1. (a)

There exists a vector $v\in\mathbb{R}^{d}_{>0}$ such that $v^{T}S_{b}=0$ and $v^{T}A^{+}(\rho_{cv})<0$ hold for all $\rho_{cv}\in\mathcal{P}_{cv}$. 2. (b)

There exists a vector $v\in\mathbb{R}^{d}_{>0}$ such that $v^{T}S_{b}=0$ and $v^{T}A^{+}(\theta)<0$ hold for all $\theta\in\textnormal{vert}(\mathcal{P}_{cv})$.

4 Structural ergodicity of unimolecular reaction networks

We are interested in this section in the structural stability of the characteristic matrix of given unimolecular network. Hence, we have in this case $\mathcal{P}_{\bullet}:=\mathbb{R}^{n_{\bullet}}_{>0}$, $\bullet\in\{u,dg,ct,cv\}$ where $n_{\bullet}$ is the dimension of the vector $\rho_{\bullet}$.

4.1 A preliminary result

Lemma 12

Let $A(\rho_{u})\in\mathbb{R}^{d\times d}$ be the characteristic matrix of some unimolecular network and $\rho_{u}\in\mathcal{P}_{u}$. Then, the following statements are equivalent:

1. (a)

For all $\rho_{dg}\in\mathcal{P}_{dg}$ and a $\rho_{cv}\in\mathcal{P}_{cv}$, the matrix $A(\rho_{dg},\rho_{cv},0)$ is Hurwitz stable. 2. (b)

The matrix $A(\mathds{1},\rho_{cv},0)$ is Hurwitz stable for all $\rho_{cv}\in\mathcal{P}_{cv}$.

*Proof : * The proof that (a) implies (b) is immediate. To prove the reverse implication, we use contraposition and we assume that there exist a $\rho_{dg}\in\mathcal{P}_{dg}$ and a $\rho_{cv}\in\mathcal{P}_{cv}$ such that $A(\rho_{dg},\rho_{cv},0)$ is not Hurwitz stable. Then, we clearly have that

[TABLE]

where $\theta=\min(\rho_{dg})$ and hence $A(\theta\mathds{1},\rho_{cv},0)$ is not Hurwitz stable. Since $A(\theta\mathds{1},\rho_{cv},0)$ is affine in $\theta$ and $\rho_{cv}$, then we have that $\theta A(\mathds{1},\rho_{cv}/\theta,0)$ and since $\theta$ is independent of $\rho_{cv}$, then we get that the matrix $A(\mathds{1},\tilde{\rho}_{cv},0)$ is not Hurwitz stable for some $\tilde{\rho}_{cv}\in\mathcal{P}_{cv}$. The proof is complete. $\diamondsuit$

4.2 Main result

Theorem 13

Let $A(\rho_{u})\in\mathbb{R}^{d\times d}$ be the characteristic matrix of some unimolecular network and $\rho_{u}\in\mathbb{R}^{n_{u}}_{>0}$. Then, the following statements are equivalent:

1. (a)

The matrix $A(\rho_{u})$ is Hurwitz stable for all $\rho_{u}\in\mathbb{R}^{n_{u}}_{>0}$. 2. (b)

*There exists a polynomial vector $v(\rho_{u})\in\mathbb{R}^{d}$ of degree at most $d-1$ such that $v(\rho_{u})>0$ and $v(\rho_{u})^{T}A(\rho_{u})<0$ for all $\rho_{u}\in\mathbb{R}^{n_{u}}_{>0}$. * 3. (c)

There exists a $\rho_{u}^{s}\in\mathbb{R}^{n_{u}}_{>0}$ such that the matrix $A(\rho_{u}^{s})$ is Hurwitz stable and the polynomial $(-1)^{d}\det(A(\rho_{u}))$ is positive for all $\rho_{u}\in\mathbb{R}^{n_{u}}_{>0}$. 4. (d)

For all $\rho_{dg}\in\mathbb{R}^{n_{dg}}_{>0}$ and a $\rho_{cv}\in\mathbb{R}^{n_{cv}}_{>0}$, the matrix $A_{\rho}:=A(\rho_{dg},\rho_{cv},0)$ is Hurwitz stable and we have that $\varrho(W_{ct}A_{\rho}^{-1}S_{ct})=0$. 5. (e)

*The matrix $A_{n}(\rho_{cv}):=A(\mathds{1},\rho_{cv},0)$ is Hurwitz stable for all $\rho_{cv}\in\mathbb{R}^{n_{cv}}_{>0}$ and $\varrho(W_{ct}A_{n}(\rho_{cv})^{-1}S_{ct})=0$ for all $\rho_{cv}\in\mathbb{R}^{n_{cv}}_{>0}$. *

Moreover, when each column of $S_{cv}$ contains exactly two nonzero entries, one being equal to $-1$ and one being equal to 1, then the above statements are also equivalent to

1. (f)

The matrix $A_{\mathds{1}}:=A(\mathds{1},\mathds{1},0)$ is Hurwitz stable and $\varrho(W_{ct}A_{\mathds{1}}^{-1}S_{ct})=0$.

*Proof : * The equivalence between the three first statements has been proved in Theorem 8. Let us prove now that (c) implies (d). Assuming that (c) holds, we get that the existence of a $\rho_{u}^{s}=\operatorname*{col}(\rho_{dg}^{s},\rho_{cv}^{s},\rho_{ct}^{s})$ such that the matrix $A(\rho_{u}^{s})$ is Hurwitz stable immediately implies that the matrix $A_{\rho}=A(\rho_{dg},\rho_{cv},0)$ is Hurwitz stable since we have that $A_{\rho}\leq A(\rho_{u})$ and, therefore $\lambda_{PF}(A_{\rho})\leq\lambda_{PF}(A(\rho_{u}))<0$. Using now the determinant formula, we have that

[TABLE]

where $D(\rho_{ct}):=\operatorname*{diag}(\rho_{ct})$ and $W_{ct}$ is defined such that ${\operatorname*{diag}(\rho_{ct})W_{ct}x}$ is the vector of propensity functions associated with the catalytic reactions. Since $A_{\rho}$ is Hurwitz stable then the determinant has fixed sign and is positive if $d$ is even, negative otherwise. Hence, this implies that

[TABLE]

for all $\rho_{ct}\in\mathbb{R}^{n_{ct}}_{>0}$ where $M:=-W_{ct}A_{\rho}^{-1}S_{ct}$. Since the matrices $W_{ct},S_{ct}$ are nonnegative, the diagonal entries of $D(\rho_{ct})$ are positive and $A_{\rho}^{-1}$ is nonpositive (since $A_{\rho}$ is Metzler and Hurwitz stable), then $M$ is nonnegative. Then, by the Perron-Frobenius theorem, we have that $\lambda_{PF}(M)=\varrho(M)$. Clearly, the fact that (15) holds for all $\rho_{ct}\in\mathbb{R}^{n_{ct}}_{>0}$ implies that $\lambda_{PF}(M)=\varrho(M)=0$ since, otherwise, there would exist a $D(\rho_{ct})$ such that $\det(I-D(\rho_{ct})M)=0$. This completes the argument.

The converse (i.e. (d) implies (c)) can be proven by noticing that if $A_{\rho}$ is Hurwitz stable, then $A_{\rho}+\epsilon S_{ct}W_{ct}$ remains Hurwitz stable for some sufficiently small $\epsilon>0$. This proves the existence of a $\rho_{u}^{s}\in\mathbb{R}^{d}_{>0}$ such that the matrix $A(\rho_{u}^{s})$. Using the determinant formula, it is immediate to see that the second statement implies the determinant condition of statement (c).

The equivalence between the statements (d) and (e) comes from Lemma 12 and the fact that the sign-pattern of the inverse of a Hurwitz stable Metzler matrix is uniquely defined by its sign-pattern.

Let us now focus on the equivalence between the statements (d) and (f) under the assumption that each column of $S_{cv}$ contains exactly one entry equal to $-1$ and one equal to 1. Assume w.l.o.g that $S_{dg}=\operatorname*{col}(-I_{n_{dg}},0)$. Then, we have that $\mathds{1}_{d}^{T}A(\rho_{dg},\rho_{cv},0)=\begin{bmatrix}-\rho_{dg}^{T}&0\end{bmatrix}$. Hence, the function $V(z)=\mathds{1}_{d}^{T}z$ is a weak Lyapunov function for the linear positive system $\dot{z}=A(\rho_{dg},\rho_{cv},0)z$. Invoking LaSalle’s invariance principle, we get that the matrix is Hurwitz stable if the matrix

[TABLE]

is Hurwitz stable for all $(\rho_{dg},\rho_{cv})\in\mathbb{R}^{n_{dg}}_{>0}\times\mathbb{R}^{n_{cv}}_{>0}$. Note that this is a necessary condition for the matrix $A(\rho_{dg},\rho_{cv},0)$ to be Hurwitz stable for all rate parameters values. Hence, this means that the stability of the matrix $A_{\rho}$ is equivalent to the Hurwitz stability of $A_{\mathds{1}}:=A(\mathds{1},\mathds{1},0)$. Finally, since $A^{22}(\rho_{dg},\rho_{cv})$ is Hurwitz stable, then we have that $\mathds{1}^{T}A^{22}(\mathds{1},\mathds{1})<0$. The proof is complete. $\diamondsuit$

5 Examples

5.1 Example 1 - SIR model

Let us consider the open stochastic SIR model considered in [13] described by the matrices

[TABLE]

where all the parameters are positive. The first, second and third states are associated with susceptible, infectious and removed people, respectively. The reaction network consists of 6 reactions: 5 unimolecular reactions including 3 degradation reactions (with rates $\gamma_{s},\gamma_{i}$ and $\gamma_{r}$) and 2 conversion reactions (with rates $k_{ir}$ and $k_{rs}$), and 1 bimolecular reaction that implements the contamination reaction that converts one susceptible person to an infectious one.

The constraint $v^{T}S_{b}=0$ enforces that $v=\tilde{v}^{T}S_{b}^{\bot}$, $\tilde{v}>0$, where $S_{b}^{\bot}=\begin{bmatrix}1&1&0\\ 0&0&1\end{bmatrix}$. This leads to

[TABLE]

Since the entries are not independent, the use of sign-matrices or interval matrices are conservative. However, if we use Theorem 13, then we can just substitute the parameters by 1 and observe that the resulting matrix is Hurwitz stable to prove the structural stability of the matrix. Alternatively, we can take $\tilde{v}=\mathds{1}$ and obtain

[TABLE]

from which the same result follows.

5.2 Example 2 - Circadian Clock

We consider the circadian clock-model of [7] which is described by the matrices

[TABLE]

where all the parameters are positive. Even if initially the network consists of 9 species and 16 reactions, the problem can be reduced to the above problem involving 5 species (the activator protein $A$ and its corresponding mRNA $M_{A}$, the activator protein $R$ and its corresponding mRNA $M_{R}$ and a dimer $C$ composed of an activator protein $A$ and a repressor protein $R$). There are 8 reactions: 7 unimolecular reactions including 4 degradation reactions (with rates $\delta_{M_{A}},\delta_{A},\delta_{M_{R}}$ and $\delta_{R}$), 2 catalytic reactions (with rates $\beta_{A}$ and $\beta_{R}$) and one conversion reaction with rate $\delta_{A}$ (with rate identical to the degradation rate of the activator protein $A$), and one bimolecular reaction corresponding to the binding reaction of $A$ and $R$.

As in the previous example, the condition reduces to

[TABLE]

where $\tilde{v}>0$. Clearly, we have four degradation reactions and two catalytic ones. Using the last statement of Theorem 13 we get that the matrix $A_{\mathds{1}}=-I$. We also have in this case that

[TABLE]

and, hence, $W_{ct}A_{\mathds{1}}^{-1}S_{ct}=0$. Hence, the system is structurally stable. Alternatively, the triangular structure of the matrix would also lead to the same conclusion.

5.3 Example 3 - Toy model

Let us consider here the following toy network where

[TABLE]

Assume that $\alpha_{1}=k_{2}$ and $\alpha_{2}=k_{3}$. In such a case, the network consists of 5 unimolecular reactions: 2 degradation reactions with rates $\gamma_{1}$ and $\gamma_{2}$ and 3 conversion reactions with rates $k_{1}$, $k_{2}$ and $k_{3}$. Then, we get that

[TABLE]

is Hurwitz stable and hence that the matrix is structurally stable. However, if we assume now that $\alpha_{1}=\alpha_{2}=0$, then the network still consists of 5 unimolecular reactions but now we have 2 degradation reactions with rates $\gamma_{1}$ and $\gamma_{2}$, 1 conversion reaction with rate $k_{1}$ and 2 catalytic reactions with rates $k_{2}$ and $k_{3}$. In such a case, we get

[TABLE]

where $A_{\mathds{1}}$ is Hurwitz stable. We have in this case that

[TABLE]

and hence

[TABLE]

which has a spectral radius equal to 1. Hence, the matrix is not structurally stable. Define now

[TABLE]

Using a perturbation argument, we can prove that the 0-eigenvalue of $A^{+}(0)$ locally bifurcates to the open left half-plane for some sufficiently small $k_{1}>0$ if and only if $k_{2}^{+}k_{3}^{+}-\gamma_{1}^{-}\gamma_{2}^{-}<0$. Hence, there exists a $k_{1}>0$ such that $A^{+}(k_{1})$ is Hurwitz stable if and only if $k_{2}^{+}k_{3}^{+}-\gamma_{1}^{-}\gamma_{2}^{-}<0$. Noting now that

[TABLE]

and, hence, the determinant never switches sign, which proves that the matrix $A^{+}(k_{1})$ is structurally stable.

6 Conclusion

Several conditions have been proposed for checking the robust and the structural ergodicity of unimolecular and a certain class of bimolecular reaction networks. When certain conditions are met, the conditions can be verified using convex programming tools and/or through the computation of simple algebraic quantities. These results complement those of [15] by providing a more accurate set of conditions. Yet, improvements are still possible at the level of Theorem 8 and Theorem 13 in order to avoid the use of complex numerical tools or to remove the need for simplifying assumptions.