On the stability of nucleic acid feedback controllers

TL;DR

This paper develops theoretical insights into the stability and dynamics of nucleic acid feedback controllers, revealing that nonlinear effects can destabilize systems even when linear designs suggest stability.

Contribution

It introduces foundational theoretical results on equilibria and stability of nucleic acid controllers, highlighting the nonlinear dynamics introduced by DNA strand displacement reactions.

Findings

Linear stability does not guarantee chemical network stability.

Nonlinear dynamics can cause instability in nucleic acid feedback systems.

Simulation confirms theoretical predictions of destabilization.

Abstract

Recent work has shown how chemical reaction network theory may be used to design dynamical systems that can be implemented biologically in nucleic acid-based chemistry. While this has allowed the construction of advanced open-loop circuitry based on cascaded DNA strand displacement (DSD) reactions, little progress has so far been made in developing the requisite theoretical machinery to inform the systematic design of feedback controllers in this context. Here, we develop a number of foundational theoretical results on the equilibria, stability, and dynamics of nucleic acid controllers. In particular, we show that the implementation of feedback controllers using DSD reactions introduces additional nonlinear dynamics, even in the case of purely linear designs, e.g. PI controllers. By decomposing the effects of these non-observable nonlinear dynamics, we show that, in general, the stability of the linear system design does not necessarily imply the stability of the underlying chemical network, which can be lost under experimental variability when feedback interconnections are introduced. We provide an in-depth theoretical analysis of an example illustrating this phenomenon, whereby the linear design does not capture the instability of the full nonlinear system implemented as a DSD reaction network, and we further confirm these results using VisualDSD, a bespoke software tool for simulating nucleic acid-based circuits. Our analysis highlight the many interesting and unique characteristics of this important new class of feedback control systems.

Figures

Tables

Table 1. Table 1: Nominal parameters for the example, and an asymmetrical parameterisation case which results in unstable dynamics.

Parameter	Nominal	Asymmetrical case
$k_1^{\pm }$	$0.001/$s	$0.00132/$s
$k_2^{\pm }$	$0.001/$s	$k_2^{+}=0.001320$/s, $k_2^{-}=0.000680$/s
${\gamma }_i^{\pm }$, $i=1,2,3$	$0.004/$s	${\gamma }_1^{\pm }={\gamma }_2^{\pm }=0.00528/$s, ${\gamma }_3^{\pm }=0.00272/$s
${\gamma }_i^{\pm }$, $i=6,7,8$	$0.008/$s	${\gamma }_6^{\pm }={\gamma }_7^{+}=0.01056/$s, ${\gamma }_7^{-}={\gamma }_8^{\pm }=0.00544/$s
${\gamma }_i^{\pm }$, $i=4,5$	$4\times {10}^{-6}/$s	${\gamma }_4^{+}=2.72\times {10}^{-6}/$s, ${\gamma }_4^{-}={\gamma }_5^{\pm }=5.28\times {10}^{-6}/$s
$k_P^{\pm }={\gamma }_4^{\pm }/{\gamma }_5^{\pm }$	$1$	$k_P^{+}=0.5152$, $k_P^{-}=1$
$k_0^{\pm }$	$0.00045/$s	$0.000594/$s
$\eta$	$5\times {10}^5/\text{M}/$s	$5\times {10}^5/\text{M}/$s

Table 2. Table 2: Poles with maximum real part, for the I/O and linearised dynamics, for the nominal and asymmetrical parametrisations.

Matrix $𝐌$	Poles corresponding to $\alpha \left\{𝐌\right\}$	Stability
${\overline{𝐑}}_{11}$	$-3.96\times {10}^{-6}$	${\overline{𝐑}}_{11}\in \mathscr{H}$
${\overline{𝐀}}_s$	$-3.96\times {10}^{-6}$	${\overline{𝐀}}_s\in \mathscr{H}$
${𝐑}_{11}$	$-5.23\times {10}^{-6}$	${𝐑}_{11}\in \mathscr{H}$
${𝐀}_s$	$+3.16\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}\pm i1.26\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$	${𝐀}_s\notin \mathscr{H}$

Table 3. Table 3: Assuming the I/O system is stable, we can state properties about CRN stability and the unforced equilibria 𝐱 ∗ superscript 𝐱 \mathbf{x}^{*} .

Parameters	Cascaded	With feedback
Nominal ${\overline{𝐑}}_{11}\in \mathscr{H}$	${𝐱}^{\ast }=0$ Unforced dynamics are GAS	Possible ${𝐱}^{\ast }>0$ Unforced dynamics are bounded
Asymmetrical ${𝐑}_{11}\in \mathscr{H}$	${𝐱}^{\ast }=0$ Unforced dynamics are GAS if additionally	Possible ${𝐱}^{\ast }>0$ CRN may be unstable

Equations

$$a_1{X}_1+a_2{X}_2\stackrel{\gamma }{\to }b{X}_3\Rightarrow {\dot{x}}_3=b\gamma x_1^{a_1}{x}_2^{a_2}$$

$$X_i\stackrel{\gamma }{\to }{X}_i+X_j,X_j\stackrel{\gamma }{\to }\varnothing ,X_i+X_j\stackrel{\eta }{\to }\varnothing$$

$$X_4^{\pm }\stackrel{k_1^{\pm }}{\to }{X}_4^{\pm }+X_5^{\pm },X_5^{\pm }\stackrel{k_2^{\pm }}{\to }\varnothing ,X_5^{+}+X_5^{-}\stackrel{\eta }{\to }\varnothing$$

$${\dot{x}}_1^{\pm }$$

$${\dot{x}}_2^{\pm }$$

$${\dot{x}}_3^{\pm }$$

$${\dot{x}}_4^{\pm }$$

$${\dot{x}}_5^{\pm }$$

$$\dot{\mathbf{p}}={\mathbf{A}}_p\mathbf{p}+{\mathbf{B}}_{p}r,y=p_5$$

$$\displaystyle\mathbf{A}_{p}=\left[\begin{array}[]{ccccc}-\gamma_{3}&0&0&0&-\gamma_{2}\\ \gamma_{4}&-\gamma_{5}&0&0&0\\ k_{0}&0&0&0&0\\ 0&\gamma_{6}&\gamma_{7}&-\gamma_{8}&0\\ 0&0&0&k_{1}&-k_{2}\end{array}\right],\,\mathbf{B}_{p}=\left[\begin{array}[]{c}\gamma_{1}\\ 0\\ 0\\ 0\\ 0\end{array}\right]$$

$$\dot{\mathbf{x}}=\left({\mathbf{A}}^{\text{boxbslash}}-\mathbf{D}\left\{\left|\mathbf{a}\right|\right\}\right)\mathbf{x}+\mathbf{Br}-\eta \left(\mathbf{Px}\right)\circ \mathbf{x}$$

$$\displaystyle\mathbf{x}=\left[\begin{array}[]{c|c}\left(\mathbf{x}^{+}\right)^{T}&\left(\mathbf{x}^{-}\right)^{T}\end{array}\right]^{T}=\left[\begin{array}[]{ccc|ccc}x_{1}^{+}&\ldots&x_{N}^{+}&x_{1}^{-}&\ldots&x_{N}^{-}\end{array}\right]^{T}$$

$$\displaystyle\mathbf{P}=\left[\begin{array}[]{cc}\mathbf{0}&\mathbf{I}\\ \mathbf{I}&\mathbf{0}\end{array}\right]\Rightarrow\left(\mathbf{Px}\right)\circ\mathbf{x}=\left[\begin{array}[]{c}\mathbf{x}^{+}\circ\mathbf{x}^{-}\\ \mathbf{x}^{+}\circ\mathbf{x}^{-}\end{array}\right]$$

$$\displaystyle\begin{array}[]{c}\dot{\mathbf{x}}^{+}=\mathbf{A}_{1}^{+}\mathbf{x}^{+}+\mathbf{A}_{2}^{-}\mathbf{x}^{-}+\mathbf{B}_{1}^{+}r^{+}-\eta\mathbf{x}^{+}\circ\mathbf{x}^{-}\\ \dot{\mathbf{x}}^{-}=\mathbf{A}_{2}^{+}\mathbf{x}^{+}+\mathbf{A}_{1}^{-}\mathbf{x}^{-}+\mathbf{B}_{1}^{-}r^{-}-\eta\mathbf{x}^{+}\circ\mathbf{x}^{-}\end{array}$$

$$\displaystyle\mathbf{A}=\left[\begin{array}[]{cc}\mathbf{A}_{1}^{+}&\mathbf{A}_{2}^{-}\\ \mathbf{A}_{2}^{+}&\mathbf{A}_{1}^{-}\end{array}\right],\>\mathbf{B}=\left[\begin{array}[]{cc}\mathbf{B}_{1}^{+}&\mathbf{0}\\ \mathbf{0}&\mathbf{B}_{1}^{-}\end{array}\right]$$

$${\mathbf{A}}_1^{\pm }={\left({\mathbf{A}}_1^{\pm }\right)}^{\text{boxbslash}}+\mathbf{D}\left\{{\mathbf{a}}_1^{\pm }\right\},{\mathbf{a}}_1^{\pm }\le 0,{\mathbf{A}}_2^{\pm }={\left({\mathbf{A}}_2^{\pm }\right)}^{\text{boxbslash}}$$

$$\displaystyle\left[\begin{array}[]{c}\mathbf{p}\\ \mathbf{q}\end{array}\right]=\left[\begin{array}[]{cc}\mathbf{I}&-\mathbf{I}\\ \hdashline\mathbf{I}&\;\;\>\mathbf{I}\end{array}\right]\mathbf{x}=\left[\begin{array}[]{c}\mathbf{W}_{p}\\ \hdashline\mathbf{W}_{q}\end{array}\right]\mathbf{x}=\mathbf{W}\mathbf{x}$$

$$\displaystyle\begin{aligned} \left[\begin{array}[]{c}\dot{\mathbf{p}}\\ \dot{\mathbf{q}}\end{array}\right]=&\left[\begin{array}[]{cc}\mathbf{R}_{11}&\mathbf{R}_{12}\\ \mathbf{R}_{21}&\mathbf{R}_{22}\end{array}\right]\left[\begin{array}[]{c}\mathbf{p}\\ \mathbf{q}\end{array}\right]+\left[\begin{array}[]{c}\mathbf{W}_{p}\\ \mathbf{W}_{q}\end{array}\right]\mathbf{B}\mathbf{r}\\ &-\frac{\eta}{2}\left[\begin{array}[]{c}0\\ \mathbf{q}\circ\mathbf{q}-\mathbf{p}\circ\mathbf{p}\end{array}\right]\end{aligned}$$

$$\begin{array}{rl}& {\mathbf{R}}_{22}=\frac{{\left({\mathbf{A}}_1^{+}+{\mathbf{A}}_1^{-}+{\mathbf{A}}_2^{+}+{\mathbf{A}}_2^{-}\right)}^{\text{boxbslash}}}{2}-\frac{\mathbf{D}\left\{\left|{\mathbf{a}}_1^{+}\right|+\left|{\mathbf{a}}_1^{-}\right|\right\}}{2}\\ & {\mathbf{R}}_{11}=\frac{{\left({\mathbf{A}}_1^{+}+{\mathbf{A}}_1^{-}-{\mathbf{A}}_2^{+}-{\mathbf{A}}_2^{-}\right)}^{\text{boxbslash}}}{2}-\frac{\mathbf{D}\left\{\left|{\mathbf{a}}_1^{+}\right|+\left|{\mathbf{a}}_1^{-}\right|\right\}}{2}\\ & {\mathbf{R}}_{12}=\frac{{\left({\mathbf{A}}_1^{+}-{\mathbf{A}}_1^{-}-{\mathbf{A}}_2^{+}+{\mathbf{A}}_2^{-}\right)}^{\text{boxbslash}}}{2}-\frac{\mathbf{D}\left\{\left|{\mathbf{a}}_1^{+}\right|-\left|{\mathbf{a}}_1^{-}\right|\right\}}{2}\\ & {\mathbf{R}}_{21}=\frac{{\left({\mathbf{A}}_1^{+}-{\mathbf{A}}_1^{-}+{\mathbf{A}}_2^{+}-{\mathbf{A}}_2^{-}\right)}^{\text{boxbslash}}}{2}-\frac{\mathbf{D}\left\{\left|{\mathbf{a}}_1^{+}\right|-\left|{\mathbf{a}}_1^{-}\right|\right\}}{2}\end{array}$$

$$\dot{\mathbf{p}}$$

$$\dot{\mathbf{q}}$$

$$k{q}_j^2+|m_{jj}|q_j-\sum _{i\ne j}{m}_{ji}{q}_i=0$$

$$\Rightarrow q_j=\frac{1}{2k}\left(-|m_{jj}|+\sqrt{m_{jj}^2+4k\sum _{i\ne j}{m}_{ji}{q}_i}\right)\ge 0$$

$$\dot{x}=-d_{1}x-{\mathbf{c}}_2\mathbf{y}+u,\dot{y}=-d_{2}y+c_{1}x$$

$$\Rightarrow \{\begin{array}{l}\begin{array}{r}{U}^{\pm }\stackrel{1}{\to }{U}^{\pm }+X^{\pm },X^{\pm }\stackrel{d_1}{\to }\varnothing ,X^{+}+X^{-}\stackrel{k}{\to }\varnothing \\ X^{\pm }\stackrel{c_1}{\to }{X}^{\pm }+Y^{\pm },Y^{\pm }\stackrel{d_2}{\to }\varnothing ,Y^{+}+Y^{-}\stackrel{k}{\to }\varnothing \\ Y^{\pm }\stackrel{c_2}{\to }{Y}^{\pm }+X^{\mp },U^{+}+U^{-}\stackrel{k}{\to }\varnothing \end{array}\end{array}$$

$$\displaystyle\Rightarrow\bar{\mathbf{R}}_{11}=\left[\begin{array}[]{cc}-d_{1}&-c_{2}\\ c_{1}&-d_{2}\end{array}\right],\quad\bar{\mathbf{R}}_{22}=\left[\begin{array}[]{cc}-d_{1}&c_{2}\\ c_{1}&-d_{2}\end{array}\right]$$

$$\begin{array}{rlr}{k}^3{q}_1^4+2k^2{d}_1{q}_1^3& +\left(d_1^2+c_2{d}_2\right)k{q}_1^2& \\ & +c_2\left(d_2{d}_1-c_2{c}_1\right)q_1=0& \end{array}$$

$${\dot{V}}_1\left\{\mathbf{q}\right\}$$

$$\dot{\mathbf{p}}$$

$$\dot{\mathbf{s}}$$

$$q_j^2-m_{jj}{q}_j-\sum _{i=1}^{j-1}{m}_{ji}{q}_i=0$$

$$\lambda \left\{{\mathbf{R}}_{11}\right\}:{\lambda }^2+\lambda \left(d_1+d_2\right)+d_1{d}_2+c_1{c}_2=0$$

$$\Rightarrow \lambda =-\frac{\left(d_1+d_2\right)}{2}\pm \frac{1}{2}\sqrt{{\left(d_1+d_2\right)}^2-4d_1{d}_2-4c_1{c}_2}$$

$$\lambda \left\{{\bar{\mathbf{R}}}_{22}\right\}:{\lambda }^2+\lambda \left(d_1+d_2\right)+d_1{d}_2-c_1{c}_2=0\Rightarrow$$

$$\lambda =-\frac{\left(d_1+d_2\right)}{2}\pm \frac{1}{2}\sqrt{{\left(d_1+d_2\right)}^2+4\left(d_1{d}_2-c_1{c}_2\right)}$$

$$\displaystyle\mathbf{Q}=\left[\begin{array}[]{cccc|cccc}1&0&\ldots&0&0&0&\dots&0\\ 0&0&\ldots&0&1&0&\ldots&0\\ \hline\cr 0&1&\dots&0&0&0&\ldots&0\\ 0&0&\ldots&0&0&1&\dots&0\\ \hline\cr\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ \hline\cr 0&0&\ldots&1&0&0&\ldots&0\\ 0&0&\ldots&0&0&0&\ldots&1\end{array}\right]$$

$$\displaystyle\mathbf{z}=\mathbf{Q}\mathbf{x}=\left[\begin{array}[]{cc|cc|c|cc}x_{1}^{+}&x_{1}^{-}&x_{2}^{+}&x_{2}^{-}&\dots&x_{N}^{+}&x_{N}^{-}\end{array}\right]^{T}$$

Full text

††thanks: Corresponding author

1 Introduction

Recent advances in synthetic biology have seen the incorporation of many control engineering design principles into the construction of biomolecular circuits [1, 2, 3, 4]. One of the current urgent needs of this aspect of synthetic biology is the development of bespoke feedback control theory that can be used to systematically design synthetic controllers for biomolecular processes. A promising direction for this work is to exploit chemical reaction network (CRN) theory, since CRN’s act as a “bridge” between mathematical designs based on ordinary differential equations (ODEs) and biological implementations in nucleic acid-based chemistry using DNA strand displacement (DSD) reactions [5, 6]. The capability of such circuits to operate in vivo and interface with endogenous cellular machinery has been demonstrated in mammalian cells, with some notable examples including engineered oligonucleotide AND gates responding to microRNA inputs [7], multi-input logic based on DNA circuitry interacting with native mRNA [8], and reliable strand displacement probes triggered by mRNA being transcribed into cells [9]. This makes circuits based on nucleic acids strong potential candidates for implementing many computing and control applications in synthetic biology.

The CRN to DNA design framework [5] assigns a formal species in the CRN to sets of DNA species, allowing the construction of circuits supported by a high level of automation using available syntax and software tools [10, 11, 12]. In the context of feedback control, however, a key challenge with employing CRNs is their inability to directly represent negative signals (since concentrations of chemical species are always positive). For example, CRNs generally can only compute a positive difference between two positive inputs, i.e. “one-sided” subtraction [13]. The use of the so-called dual-rail representation with nucleic-acids [14] circumvents this problem by representing each signal as the difference of concentrations of two different species. Although it doubles the number of required reactions, the dual rail representation enables the computation of rational linear functions as the steady state of a CRN [15], including two-sided subtraction. It provides an Internally Positive Representation (IPR), where a positive state-space system, together with input, state and output transformations, can realize arbitrary input/output dynamics [16]. Unimolecular reactions of catalysis and degradation, and bimolecular reactions of annihilation, can then be used to construct CRNs to approximately represent transfer functions [17], linear feedback systems [18, 12, 19], and nonlinear controllers [20].

In all these systems, annihilation reactions operating with very fast timescales are essential, in order to ensure that species concentrations remain within the bounds of experimental feasibility. However, these reactions result in a nonlinear IPR. As noted in [21], these annihilation reactions introduce additional internal nonlinear dynamics that are not observable in the represented input/output linear dynamics, but become important in the presence of inevitable experimental variability in the biomolecular implementations. Here, we formally characterise the effects of the nonlinear dynamics introduced through these annihilation reactions on the equilibria and stability of closed-loop nucleic acid systems. These results provide many useful insights that can guide the design and construction of these circuits in vitro and in vivo, and also highlight some of the associated technical challenges and limitations.

1.1 Notation and Preliminaries

We represent the elements of vectors and matrices $\mathbf{x=Mv}$ with $x_{j}=\left[\mathbf{Mv}\right]_{j}=\sum_{i}m_{ji}v_{i}$. $\mathbf{1}$ is a vector with elements $1$, and $\mathbf{I}$ is the identity matrix. The element-wise product is represented with $\mathbf{x}=\mathbf{v}\circ\mathbf{u}\Rightarrow x_{j}=v_{j}u_{j}$. For a vector $\mathbf{v}\geq 0$, $\|\mathbf{v}\|_{1}=\mathbf{1}^{T}\mathbf{v}$ and $\|\mathbf{v}\|_{2}^{2}=\mathbf{1}^{T}\left(\mathbf{v}\circ\mathbf{v}\right)$. In the system dynamics, for brevity, time dependency is implicit, i.e. $x_{j}\equiv x_{j}(t)$, $x_{j}^{*}\equiv x_{j}(\infty)$ represents steady state conditions, and $X_{j}(s)$ is the Laplace transform of $x_{j}$. $\rho\left\{\mathbf{M}\right\}$ denotes the set of the eigenvalues $\lambda_{i}$ of matrix $\mathbf{M}$. We represent the set of Hurwitz matrices with $\mathcal{H}$. Given the spectral abscissa $\alpha\left\{\mathbf{M}\right\}=\max_{i}\Re\left\{\lambda_{i}\right\}$, if $\mathbf{M}\in\mathcal{H}$, then $\alpha\left\{\mathbf{M}\right\}<0$. Given the set of lower triangular matrices $\mathcal{L}$, then for $\mathbf{M}\in\mathcal{L}$ we have that $m_{ji}=0,i>j$, and $\lambda_{i}\left\{\mathbf{M}\right\}=m_{ii}$. Given the set $\mathcal{I}$ of irreducible matrices [22], if $\mathbf{M}\in\mathcal{I}$, then there is no permutation such that $\mathbf{M}\in\mathcal{L}$. Also, if $\mathbf{M}\in\mathcal{L}$, then $\mathbf{M}\notin\mathcal{I}$.

$\mathbb{R}_{0}^{+}$ is the positive orthant, where all the coordinates of a vector $v_{j}\geq 0$. $\mathbf{M}\geq 0$ means all elements $m_{ji}\geq 0$, and $\mathbf{M}\in\mathbb{R}_{0}^{+}$. The operator $\mathbf{D}\left\{\mathbf{v}\right\}$ is defined as a square matrix where $d_{jj}=v_{j}$ and $d_{ji}=0,\,j\neq i$. If $\mathbf{m}$ is the diagonal of $\mathbf{M}$, the matrix of off-diagonal elements $\mathbf{M}^{\boxbslash}$ is defined as $\mathbf{M}^{\boxbslash}=\mathbf{M}-\mathbf{D}\left\{\mathbf{m}\right\}$. Defining $\mathcal{M}$ as the group of Metzler matrices, if $\mathbf{M}\in\mathcal{M}$, then $\mathbf{M}^{\boxbslash}\geq 0$ and if $\mathbf{M}\in\mathcal{M,H}$ then $\mathbf{m}<0$.

A CRN is composed of a set of reactions between chemical species $X_{j}$, at a rate $\gamma$. The dynamics of a chemical reaction can be approximated by ODEs using mass action kinetics [23], i.e.,

[TABLE]

The stoichiometric coefficients $a_{1}$, $a_{2}$ and $b$ indicate, respectively, the relative number of molecules consumed and produced during the reaction. Note that expressing the dynamics of (1) in their natural coordinates, the concentrations, results in a non-negative state vector $\mathbf{x}\in\mathbb{R}_{0}^{+}$. The use of chemical concentrations as state variables is therefore not suitable for circuits involving negative signals, such as feedback control loops. To circumvent the above problem, it is now standard practice to represent both positive and negative signals with a dual-rail representation, where each signal is the difference between two positive quantities. Consider, for example, the linear feedback control system depicted in Fig. 1. In the dual-rail representation, all signals in this system are split into contributions from two molecular concentrations $p_{j}=x_{j}^{+}-x_{j}^{-}$, where $x_{j}^{+}\geq 0$, $x_{j}^{-}\geq 0$, and $p_{j}\in\mathbb{R}$. Linear mathematical operators and transfer functions can be represented with elementary CRNs [18, 17] using unimolecular reactions of catalysis and degradation, and bimolecular reactions of annihilation, which are respectively given by

[TABLE]

Employing the dual rail representation entails duplicating all the catalysis and degradation reactions. In the following we compact the notation so that $X^{\pm}$ represents simultaneously both species $X^{+}$ and $X^{-}$, and $x^{\pm}$ their respective two concentrations $x^{+}$ and $x^{-}$. $Y^{\pm}\xrightarrow{\gamma^{\pm}}Y^{\pm}+X^{\mp}$ is an abbreviation for the two parallel reactions $Y^{+}\xrightarrow{\gamma^{+}}Y^{+}+X^{-}$ and $Y^{-}\xrightarrow{\gamma^{-}}Y^{-}+X^{+}$.

2 Representation of linear feedback control systems with chemical reaction networks

Each reaction in (2) has equivalent representation with DSD reactions, and the sets of these reactions can be systematically converted to implementable reactions based on nucleic acids (for example with [5]). Since the dual representation admits infinite combinations of $x_{j}^{+}$ and $x_{j}^{-}$ for the same difference $p_{j}=x_{j}^{+}-x_{j}^{-}$, in practice, the annihilation reaction in (2) is used to ensure one of the concentrations is kept close to zero, and $p_{j}\approx x^{+}_{j}$ or $p_{j}\approx-x_{j}^{-}$. A rate $\eta$ for the annihilation reactions that is on a much faster timescale than the dynamics of the system is used to keep the concentrations of all molecular species low (i.e. experimentally feasible) even in the presence of transients.

Assumption 1**.**

The nominal parameterisation and nominal implementation assume perfectly designed reaction rates in the absence of variability, and a symmetrical parametrisation where the reaction rates are the same in each dual reaction with $\gamma_{j}^{+}=\gamma_{j}^{-}=\gamma_{j}$.

Definition 1**.**

The Input-Output (I/O) dynamics are the response $Y(s)=G(s)R(s)$ from $r=\left(r^{+}-r^{-}\right)$ to $y=\left(y^{+}-y^{-}\right)$, with $r,y\in\mathbb{R}$ and $r^{\pm},y^{\pm}\in\mathbb{R}^{+}_{0}$.

For example, in Fig. 1 we represent the plant with the following set of chemical reactions

[TABLE]

Computing the I/O dynamics from Definition 1 under Assumption 1 ($k_{i}^{+}=k_{i}^{-}=k_{i}$), the nonlinear terms cancel out and the I/O dynamics of (3) represent a first-order linear system $Y(s)=k_{1}\left(s+k_{2}\right)^{-1}U(s)$, with $y=x_{5}^{+}-x_{5}^{-}$ and $u=x_{4}^{+}-x_{4}^{-}$. The use of bimolecular reactions results in an IPR of a linear system based on nonlinear internal positive dynamics, in contrast to IPRs based on linear positive dynamics [16].

The linear feedback system in Fig. 1 is represented by combining (3) with the CRNs of the linear operations of integration $\dot{p}_{3}=k_{0}p_{1}$, gain $p_{2}=k_{P}p_{1}$, and summation $p_{4}=p_{2}+p_{3}$. Despite the positivity of the concentrations, the dual rail representation allows us to represent the error $p_{1}=r-p_{5}$ with a CRN (see e.g. [18, 12, 19]).

Example 2.1**.**

The mass action kinetics of the linear feedback structure in Fig 1 results in the following dynamics

[TABLE]

Under Assumption 1, the ODEs for $\dot{p}_{j}=\dot{x}^{+}_{j}-\dot{x}^{-}_{j}$ and $r=r^{+}-r^{-}$ result in the I/O system

[TABLE]

Remark 2.2**.**

The I/O dynamics are an approximation to the original linear system in Fig. 1, since the representations of subtraction, gain, and sum are exact only at steady state. The impact of the respective transient dynamics (4a), (4b), and (4d) can be mitigated by increasing $\gamma_{j}$.

3 Dynamics of the chemical reaction network

We now define the class of systems analysed in this work, where we retain the natural non-negative coordinates, so that the states are the species concentrations $x_{j}^{\pm}$, and the input vector contains both positive and negative components for the reference $\mathbf{r}=\left[r^{+},r^{-}\right]^{T}$, $r^{\pm}\in\mathbb{R}_{0}^{+}$.

Assumption 2**.**

*Assume the dynamics we wish to represent result in stable I/O dynamics, and therefore $\mathbf{A}_{p}\in\mathcal{H}$ and $\mathbf{A}_{p}^{-1}$ exists. *

3.1 The dynamics in the natural coordinates are positive and nonlinear

Definition 2**.**

Defining the state $\mathbf{x}\in\mathbb{R}_{0}^{+}$ as the vector of species concentrations, the mass action kinetics of the constructed CRN result in

[TABLE]

where $\mathbf{A}=\mathbf{A}^{\boxbslash}-\mathbf{D}\left\{\left|\mathbf{a}\right|\right\}$ and

[TABLE]

Compared to (5), the model in (17) includes the components from the bimolecular reactions $-\eta\mathbf{x}^{+}\circ\mathbf{x}^{-}$. The unimolecular reactions depend linearly on the state with $\mathbf{Ax}$, where by construction the catalysis rates end up on the off-diagonal elements $\mathbf{A}^{\boxbslash}\geq 0$ and the degradation rates result in non-positive elements in the diagonal of $\mathbf{D}\left\{\mathbf{a}\right\}$ ($\mathbf{a}\leq 0$). $\mathbf{A}$ is Metzler, and we can decompose the dynamics into non-negative and non-positive contributions where $\mathbf{D}\left\{\mathbf{a}\right\}-\eta\left(\mathbf{Px}\right)\circ\mathbf{x}\leq 0$, and $\mathbf{A}^{\boxbslash}\mathbf{x}+\mathbf{Br}\geq 0$. With $\mathbf{B}\geq 0$, $\mathbf{r}\geq 0$, and $\mathbf{g}\left\{\mathbf{x}\right\}=-\eta\left(\mathbf{Px}\right)$, the following Lemma 3.3 shows that the nonlinear dynamics in their natural coordinates in (17) are non-negative.

Lemma 3.3**.**

For a vector function $\mathbf{g}\left\{\mathbf{x}\right\}$, if $\mathbf{M}\in\mathcal{M}$, $\mathbf{v}\geq 0$, and $\mathbf{x}\left(0\right)>0$, the dynamics $\dot{\mathbf{x}}=\mathbf{Mx}+\mathbf{x}\circ\mathbf{g}\left\{\mathbf{x}\right\}+\mathbf{v}$ are non-negative.

{pf}

For each component $\dot{x}_{j}=\left[\mathbf{Mx}\right]_{j}+x_{j}\left[\mathbf{g}\left\{\mathbf{x}\right\}\right]_{j}+v_{j}$. If $x_{j}=0$ and $\exists_{i\neq j}:x_{i}>0$, then $\dot{x}_{j}=\left[\mathbf{Mx}\right]_{j}+\mathbf{v}\geq 0$ and the trajectory remains in $\mathbb{R}^{+}_{0}$. ∎ Rewriting (17) according to the partition in (20) so that

[TABLE]

we have matrices $\mathbf{A}$ and $\mathbf{B}$ structured into

[TABLE]

$\mathbf{A}_{j}^{\pm}\in\mathcal{M}$, and from Definition 2 the degradation rates are in the diagonal of $\mathbf{A}_{1}^{\pm}$, while $\mathbf{A}_{2}^{\pm}$ contain only the catalysis reactions used to represent subtraction. Because the catalysis and degradation reactions are duplicated, both matrices $\mathbf{A}_{i}^{\pm}$ retain the same structure, but not necessarily the same parameterisation (similarly for the pair $\mathbf{B}_{1}^{\pm}$). Matrices $\mathbf{A}_{i}^{+}$ and $\mathbf{B}_{1}^{+}$ are populated with the reaction rates $\gamma_{j}^{+}$, and their counterparts $\mathbf{A}_{i}^{-}$ and $\mathbf{B}_{1}^{-}$ with $\gamma_{i}^{-}$.

3.2 The positive nonlinear dynamics are unobservable in the I/O dynamics of the linear representation

Definition 3**.**

The rotated coordinates $p_{j}=x_{j}^{+}-x_{j}^{-}\in\mathbb{R}$ and $q_{j}=x_{j}^{+}+x_{j}^{-}\in\mathbb{R}_{0}^{+}$ result from the similarity transformation $\mathbf{W}$, where

[TABLE]

We then have that $\mathbf{W}^{-1}=\frac{1}{2}\mathbf{W}^{T}$, $\mathbf{W}_{p}\left(\left(\mathbf{Px}\right)\circ\mathbf{x}\right)=0$ and $\mathbf{W}_{q}\left(\left(\mathbf{Px}\right)\circ\mathbf{x}\right)=2\eta\left(\mathbf{x}^{+}\circ\mathbf{x}^{-}\right)$. We can use $\mathbf{W}$ to split the dynamics into the I/O dynamics $\dot{\mathbf{p}}$ from Section 2 and the remaining nonlinear positive dynamics $\dot{\mathbf{q}}$, and infer their interconnections. The rotated dynamics are then given by

[TABLE]

Remark 3.4**.**

From the structures in Definition 2 and (29) (recall that $\mathbf{a}_{1}^{\pm}\leq 0$), we have that

[TABLE]

The diagonal of ${\mathbf{R}}_{22}$ is nonpositive, given by the average of the diagonals of $\mathbf{A}^{\pm}_{1}$. Also $\mathbf{A}_{j}^{\pm}\in\mathcal{M}\Rightarrow\mathbf{R_{22}}\in\mathcal{M}$.

Definition 4**.**

*Considering the condition of perfectly identical reaction rates from Assumption 1, we define the nominal matrices (represented with an upper bar), where we have that $\mathbf{A}_{1}^{\pm}=\bar{\mathbf{A}}_{1}$, $\mathbf{A}_{2}^{\pm}=\bar{\mathbf{A}}_{2}$, $\mathbf{B}_{1}^{\pm}=\bar{\mathbf{B}}_{1}$. *

Proposition 3.5**.**

For the nominal symmetrical parameterisation in Definition 4, the nonlinear dynamics are unobservable in the I/O system, due to the serial structure of the nominal rotated dynamics given by

[TABLE]

{pf}

Applying Definition 4 to the matrices in Remark 3.4, it follows immediately that $\bar{\mathbf{R}}_{12}=\bar{\mathbf{R}}_{21}=0$, $\bar{\mathbf{R}}_{11}=\bar{\mathbf{A}}_{1}-\bar{\mathbf{A}}_{2}$, $\bar{\mathbf{R}}_{22}=\bar{\mathbf{A}}_{1}+\bar{\mathbf{A}}_{2}$, and thus the serial structure of (38a-38b) (illustrated in Fig. 2) means that $\mathbf{p}$ evolves independently of $\mathbf{q}$, making $\mathbf{q}$ unobservable in any output of the I/O dynamics.∎

4 Equilibria of the chemical reaction network

We now compare the equilibria of the CRN with and without feedback, to analyse how feedback changes the fundamental properties of the system.

Definition 5**.**

We define a cascaded system as a set of DSD reactions without feedback, where the catalysis reactions do not depend directly or indirectly on the chemical species downstream.

Cascaded strand displacement reactions are well suited to systematically build large computational and logic gate circuitry [24]. The cascaded structure of the represented linear system results in a state matrix which can be permuted so that $\bar{\mathbf{R}}_{11}=\mathbf{A}_{p}\in\mathcal{L}$. Under Assumptions 1 and 2, and from Remark 3.4, we have $\bar{\mathbf{R}}_{11}\in\mathcal{L}\Rightarrow\bar{\mathbf{A}}_{1},\bar{\mathbf{A}}_{2}\in\mathcal{L}$, and $\bar{\mathbf{R}}_{11}\in\mathcal{L,H}\Rightarrow\bar{\mathbf{R}}_{22}\in\mathcal{L,H}$.

For example, representing the open loop of Fig. 1 without the reactions $X_{5}^{\pm}\xrightarrow{\gamma_{2}}X_{5}^{\pm}+X_{1}^{\mp}$ results in the cascade of serial and parallel reactions in Fig. 3a. In this particular case it also results in $\bar{\mathbf{A}}_{2}=0$, but in general, we can have $\bar{\mathbf{A}}_{2}\geq 0$ if there are subtractions in the cascaded I/O dynamics. Including feedback in the I/O dynamics leads to feedback within the network, and the cascaded structure is lost. The reactions $X_{5}^{\pm}\xrightarrow{\gamma_{2}}X_{5}^{\pm}+X_{1}^{\mp}$ connect the output to the input of the open loop cascade of reactions, and mass is transferred back into the input of the CRN.

Due to the triangular structure, the equilibrium of the unforced dynamics can be easily computed sequentially for each coordinate to show that there is a unique equilibrium at $\mathbf{q}=0$ (Lemma A.20 in Appendix). In the presence of feedback this is no longer possible since the states will depend on the output, and it follows that $\exists_{i>j}:\left[{\mathbf{A}}^{\pm}_{2}\right]_{ji}>0$, and $\bar{\mathbf{A}}_{2}\notin\mathcal{L}$. Consequently, the states involved in the closed loop become interdependent, and $\bar{\mathbf{R}}_{11}=\bar{\mathbf{A}}_{1}-\bar{\mathbf{A}}_{2}$ cannot be a lower triangular matrix.

Remark 4.6**.**

The interdependent evolution of all the states is reflected in the irreducibility [22] of the state matrix $\bar{\mathbf{R}}_{22}$. If $\bar{\mathbf{R}}_{22}\in\mathcal{I},\mathcal{M}$, for each coordinate $j$, $\exists_{i}:\left[\bar{\mathbf{R}}_{22}\right]_{ji}>0$. Therefore the trajectory of $q_{j}$ will always depend on another coordinate $q_{i}$, making the network irreducible.

Proposition 4.7**.**

Consider $\mathbf{M}\in\mathcal{I},\mathcal{M}$ such that $\mathbf{M}=\mathbf{M}^{\boxbslash}+\mathbf{D}\left\{\mathbf{m}\right\}$, $\mathbf{m}\leq 0$, and the dynamics $\dot{\mathbf{q}}=\mathbf{M}\mathbf{q}-k\mathbf{q}\circ\mathbf{q}$ with equilibrium $\mathbf{q}^{*}$. Then we have the following: i) $\exists_{j}q_{j}^{*}=0\Rightarrow q_{i\neq j}^{*}=0$; ii) the unforced dynamics may admit a second positive equilibrium $\mathbf{q}^{*}>0$, proportional to $k^{-1}$.

{pf}

From the equilibrium condition for each coordinate $j$ we take the non-negative roots

[TABLE]

i) If $\sum_{i\neq j}m_{ji}q_{i}=0$, then $q_{j}=0$, and we disregard the negative solution $q_{j}=-|m_{jj}|/k$. Since $\mathbf{M}\in\mathcal{I}$, for every coordinate $j$, $\exists_{l\neq j}:m_{jl}>0$, and $q_{j}=0\Leftrightarrow q_{l}=0$. We also have that for any $i\neq j:m_{ij}>0$, $q_{i}=0\Leftrightarrow q_{j}=0$. Hence, if $q_{j}=0\Rightarrow\forall_{i\neq j},\>q_{i}=0$, and we cannot have an equilibrium where only some of the states are at zero.

ii) If $\exists_{i\neq j}:m_{ji}>0$ and the coordinate $i$ is at a positive equilibrium $q_{i\neq j}^{*}>0$, then $\sum_{i\neq j}m_{ji}q_{i}^{*}>0$. The non-negative roots for each coordinate $j$ result from solving the system (40). Note that even if $m_{jj}=0$, then $q_{j}^{*}>0$.

Combining i) and ii), if $\mathbf{M}\in\mathcal{I}$, the system may have a positive equilibrium $\mathbf{q}^{*}>0$, which can be scaled down with $k$, since $\lim_{k\rightarrow\infty}q_{j}^{*}=0$.∎

Example 4.8**.**

Consider the CRN representation of a linear system with a single input $u$, which has negative feedback between its states $x$ and $y$ ($c_{2}>0$), resulting in

[TABLE]

Without feedback, so that $c_{2}=0$, the system simplifies to a reducible serial cascade where $\bar{\mathbf{R}}_{11}=\bar{\mathbf{R}}_{22}=\bar{\mathbf{A}}_{1}\in\mathcal{L}$, and the unforced dynamics $\dot{\mathbf{q}}=\bar{\mathbf{R}}_{22}\mathbf{q}-k\mathbf{q}\circ\mathbf{q}$ have a single non-negative equilibrium at $\mathbf{q}=0$.

With feedback, so that $c_{2}>0$, we can replace $q_{2}=c_{2}^{-1}\left(kq_{1}+d_{1}\right)q_{1}$ in the equilibrium conditions for $q_{1}$ and obtain the polynomial

[TABLE]

Using Descartes’ rule of signs, if $c_{2}>{d_{2}d_{1}}{c_{1}^{-1}}$, we have one positive root and the equilibrium $q_{1}^{*}>0$ exists.

Remark 4.9**.**

*Note that the use of $\bar{\mathbf{A}}_{2}$ to represent negative feedback in the I/O dynamics (5) with ${\mathbf{A}}_{p}=\bar{\mathbf{R}}_{11}=\bar{\mathbf{A}}_{1}-\bar{\mathbf{A}}_{2}$, results in positive feedback in the nonlinear dynamics in (38b) with $\bar{\mathbf{R}}_{22}=\bar{\mathbf{A}}_{1}+\bar{\mathbf{A}}_{2}$ *

In (47), $c_{2}$ impacts the spectral radius of $\bar{\mathbf{R}}_{11}$ and $\bar{\mathbf{R}}_{22}$ differently. From their characteristic polynomials, we have stable I/O dynamics ($\bar{\mathbf{R}}_{11}\in\mathcal{H}$) for any $c_{2}>0$, but for a sufficiently high gain $c_{2}>{d_{2}d_{1}}{c_{1}^{-1}}$, we get $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$. Not coincidentally, it is the same domain for which $\mathbf{q}^{*}>0$ exists. Further details in the Appendix.

Remark 4.10**.**

The existence of positive equilibrium conditions for linear feedback systems has direct consequences for the experimental construction of these circuits. Operating at an equilibrium corresponding to high concentrations aggravates leaky reactions, where undesired triggering of strand displacement leads to unwanted outputs in the absence of inputs. Furthermore, if ${\mathbf{q}}^{*}\geq 0$ with input $\mathbf{r}=0$, then the reactions persist even if the I/O dynamics are at rest $\mathbf{p}=0$, leading to unnecessary, irreversible, and costly consumption of fuel species. This is in direct contrast to cascaded DSD reactions, where without input to the I/O dynamics, the CRN is at equilibrium at $\mathbf{x}=0$, and no reactions occur.

5 Stability

We begin by proving the following lemma which is applicable to the unforced dynamics of (17) and (38b).

Lemma 5.11**.**

If $\mathbf{M}\in\mathcal{M},\mathcal{H}$, and $\mathbf{g}\left\{\mathbf{x}\right\}<0$ for $\mathbf{x}>0$, then the system $\dot{\mathbf{x}}=\mathbf{Mx}+\mathbf{x}\circ\mathbf{g}\left\{\mathbf{x}\right\}$ is globally asymptotically stable (GAS) at $\mathbf{x}=0$.

{pf}

From the stability of Metzler matrices [22], $\mathbf{M}\in\mathcal{M,H}\Rightarrow\exists_{\mathbf{d}>0}:\mathbf{M}^{T}\mathbf{D}\{\mathbf{d}\}+\mathbf{D}\{\mathbf{d}\}\mathbf{M}=-\mathbf{I}$. We take the Lyapunov function $V_{d}\left\{\mathbf{x}\right\}=\mathbf{x}^{T}\mathbf{D}\left\{\mathbf{d}\right\}\mathbf{x}>0$, and since $\mathbf{D\{d\}}\left(\mathbf{x}\circ\mathbf{g}\left\{\mathbf{x}\right\}\right)=\mathbf{d}\circ\mathbf{x}\circ\mathbf{g}\left\{\mathbf{x}\right\}<0$, $\forall_{\mathbf{x}>0}$, we have that $\dot{V}_{d}\left(\mathbf{x}\right)=-\mathbf{I}+2\mathbf{g}\left\{\mathbf{x}\right\}^{T}\left(\mathbf{d}\circ\mathbf{x}\circ\mathbf{x}\right)<0$∎ With $\mathbf{g}\left\{\mathbf{x}\right\}=-\mathbf{Px}$, Lemma 5.11 ensures that if the network of catalysis and degradation reactions is stable, $\mathbf{A}\in\mathcal{H}$, the bimolecular reactions cannot destabilise (17). A stable CRN with $\mathbf{A}\in\mathcal{H}$ can occur if the degradation of each species is faster than their overall production, and $\mathbf{A}$ has a dominant diagonal. However, this is not the general case. The dynamics without the bimolecular reactions result in the positive feedback loop between two positive systems of Fig. 3b. Since we cannot stabilise non-negative systems with positive gains $\mathbf{A}_{2}^{\pm}\geq 0$ [25], it is sufficient to have $\mathbf{A}_{1}^{\pm}\notin\mathcal{H}$ to give $\mathbf{A}\notin\mathcal{H}$. Even for the nominal symmetrical parameterisation, the representation has modes that are not present in the original linear system $\rho\left\{\bar{\mathbf{R}}\right\}=\rho\left\{\bar{\mathbf{R}}_{11}\right\}\cup\rho\left\{\bar{\mathbf{R}}_{22}\right\}$. While this is a problem for IPR with linear positive systems [16], the presence of the bimolecular reactions are sometimes sufficient for stabilisation, even if $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$.

5.1 The I/O dynamics determine the stability for the nominal symmetrical case

While at first glance it seems precarious to have unobservable nonlinear dynamics, for the nominal symmetrical case in Definition 4, it is possible to provide guarantees for stability and boundedness.

Proposition 5.12**.**

The cascaded systems from Definition 5 representing stable I/O dynamics, have GAS unforced nonlinear dynamics, for $\mathbf{x}>0$.

{pf}

From Remark 3.4, in cascaded systems $\bar{\mathbf{R}}_{11},\bar{\mathbf{R}}_{22}\in\mathcal{L}$, and $\rho\left\{\bar{\mathbf{R}}_{11}\right\}=\rho\left\{\bar{\mathbf{R}}_{22}\right\}$. If the I/O system is stable, then $\alpha\left\{\bar{\mathbf{R}}_{11}\right\}=\alpha\left\{\bar{\mathbf{R}}_{22}\right\}<0$ and Lemma 5.11 ensures $\mathbf{\dot{q}}=\mathbf{\bar{R}_{22}}\mathbf{q}-\frac{\eta}{2}\mathbf{q}\circ\mathbf{q}$ is GAS at $\mathbf{q}=0$.∎

Remark 5.13**.**

We can apply Proposition 5.12 to the representation of individual linear operations, which by themselves are cascaded reactions. It results directly that the CRNs for summation, gain, and subtraction by themselves, have GAS unforced dynamics, and are bounded for bounded inputs. More importantly, applying it to CRNs assembled from those linear operations in a cascade fashion, results in a single stable equilibrium for the complete circuit.

Recalling that with the introduction of feedback, we lose the cascaded structure and create an irreducible system, even for the representation of stable I/O linear dynamics ($\bar{\mathbf{R}}_{11}\in\mathcal{H}$), if feedback leads to $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$, then the following Lemma states that unforced trajectories diverge away from the origin due to a diverging mode of $\bar{\mathbf{R}}_{22}$.

Lemma 5.14**.**

For the dynamics $\dot{\mathbf{q}}={\mathbf{M}}\mathbf{q}-k\mathbf{q}\circ\mathbf{q}$ with ${\mathbf{M}}\in\mathcal{M,I}$ but ${\mathbf{M}}\notin\mathcal{H}$, the equilibrium at the origin $\mathbf{q}=0$ is unstable.

{pf}

From applying the Frobenius-Perron theorem to Metzler matrices [22, 26], $\mathbf{M}\in\mathcal{M,I}\Rightarrow\exists_{\mathbf{w_{F}>0}}:\mathbf{w}_{F}^{T}{\mathbf{M}}=\lambda_{F}\mathbf{w}_{F}^{T}$ and $\lambda_{F}=\alpha\left\{\mathbf{M}\right\}$. Defining the Lyapunov function $V_{F}\left(\mathbf{q}\right)=\mathbf{w}_{F}^{T}\mathbf{q}$, we have that $\mathbf{q}>0\Rightarrow V_{F}\left(\mathbf{q}\right)>0$ and $\dot{V}_{F}\left(\mathbf{q}\right)=\mathbf{w}_{F}^{T}\dot{\mathbf{q}}=\mathbf{w}_{F}^{T}\left(\mathbf{q}\circ\left(\lambda_{F}\mathbf{1}-k\mathbf{q}\right)\right)$. Since ${\mathbf{M}}\notin\mathcal{H}\Rightarrow\lambda_{F}>0$, hence $\forall_{j},\,q_{j}<\frac{\lambda_{F}}{k}$ gives that $\dot{V}_{F}\left\{\mathbf{q}\right\}>0$, and the system is divergent close to the origin. ∎ The IPR of a stable system using only linear positive systems is therefore not guaranteed to be stable [16]. However, for the nonlinear positive dynamics (38), we can still ensure boundedness with the following result.

Lemma 5.15**.**

For $\mathbf{M}\in\mathcal{M}$ , $\mathbf{q}\left(0\right)>0$, and a bounded input $\mathbf{v}\geq 0$, if $\mathbf{g}\left\{\mathbf{q}\right\}\leq-k\mathbf{q}$ then the non-negative trajectories of $\dot{\mathbf{q}}=\mathbf{M}\mathbf{q}+\mathbf{q}\circ\mathbf{g}\left\{\mathbf{q}\right\}+\mathbf{v}$ are bounded by $\|\mathbf{q}\|_{2}<k^{-1}\left(\sqrt{N}\|\mathbf{M}\|_{2}+\|\mathbf{v}\|_{1}\|\mathbf{q}\|_{2}^{-1}\right)$

{pf}

Lemma 3.3 guarantees that the trajectories are nonnegative for $\mathbf{q}\left(0\right)>0$. If $\mathbf{M}\in\mathcal{H}$ is Hurwitz, Lemma 5.11 guarantees that the system is asymptotically stable in $\mathbb{R}_{0}^{+}$ with equilibrium at $\mathbf{q}=0$. If $\mathbf{M}\notin\mathcal{H}$ is not Hurwitz, we can still show boundedness, using the linear Lyapunov function $V_{1}\left\{\mathbf{q}\right\}=\|\mathbf{q}\|_{1}=\sum_{j}q_{j}>0$, in the domain $\mathbf{q}>0$. We then have

[TABLE]

We can always find large enough values of $\mathbf{q}$ such that $\|\mathbf{q}\|_{2}>\frac{\sqrt{N}}{k}\|\mathbf{M}\|_{2}+\frac{1}{k}\frac{\|\mathbf{v}\|_{1}}{\|\mathbf{q}\|_{2}}$ where we have $\dot{V}_{1}\left\{\mathbf{q}\right\}<0$. ∎ Applying Lemma 5.15 with $\mathbf{g}\left\{\mathbf{q}\right\}=-\frac{\eta}{2}\mathbf{q}$ to the unforced dynamics in (38b) we have $\|\mathbf{q}\|_{2}<\eta^{-1}2\sqrt{N}\|\bar{\mathbf{R}}_{22}\|_{2}$ (see illustration in Fig. 10 in Appendix). In general, Lemma 5.15 is not applicable to the nonlinear dynamics (17), due to the matrix $\mathbf{P}$. Moreover, it relies on the assumption of a stable $\bar{\mathbf{R}}_{11}$.

Proposition 5.16**.**

Consider the nominal dynamics in (38a-38b), with the symmetrical parameterisation from Assumption 1. Under Assumption 2, the I/O dynamics (38a) are stable, and the concentrations in the complete CRN are bounded and can be scaled down with a faster annihilation reaction rate $\eta$.

{pf}

Assumption 2 ensures the trajectories of $\mathbf{p}$ are bounded. We can treat $\mathbf{p}$ as an additional input to the system (38b) and apply Lemma 5.15 with $\mathbf{v}=\mathbf{W}_{q}\bar{\mathbf{B}}\mathbf{r}+\frac{\eta}{2}\mathbf{p}\circ\mathbf{p}$. The unobserved dynamics are then bounded for bounded inputs $\mathbf{r},\mathbf{p}>0$, and are scaled down by increasing $\eta$. ∎ The same feedback responsible for a stable I/O linear dynamics can result in $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$ (see Remark 4.9). Designing feedback to ensure that $\bar{\mathbf{R}}_{11},\bar{\mathbf{R}}_{22}\in\mathcal{H}$ is impractical since it would put constraints on which I/O systems could be represented. It is one of the challenges of representing stable linear systems relying only on linear positive systems [16], where we would need $\bar{\mathbf{A}}\in\mathcal{H}$ for the IPR to be stable. Lemma 5.15 lifts this constraint, albeit at the cost of a positive equilibrium.

Remark 5.17**.**

With the introduction of feedback, the concentrations involved in the irreducible parts of the CRN will have positive equilibria, and $\exists_{j}\;q_{j}(t)>0$ even if $\mathbf{r}=0$ and the I/O dynamics are stable $\alpha\left\{\bar{\mathbf{R}}_{11}\right\}<0$. This result also explains why in experimental practice the annihilation rate $\eta$ is set as high as possible, to minimise the concentrations in the circuit during operation or at equilibrium.

Remark 5.18**.**

*In the presence of integrators $\bar{\mathbf{a}}_{1}\leq 0$, it is not possible to use positive feedback $\bar{\mathbf{A}}_{2}$ such that $\bar{\mathbf{A}}_{1}+\bar{\mathbf{A}}_{2}$ becomes Hurwitz [27, 25]. Starting from a marginally stable state matrix $\alpha\left\{\bar{\mathbf{A}}_{1}\right\}=0$, the introduction of feedback leads to $\alpha\left\{\bar{\mathbf{R}}_{22}\right\}\geq 0$.This raises an interesting tradeoff, when controllers that introduce integrators in the loop transfer function (for example in PI control) lead to a positive equilibrium, which is inconvenient for implementation. *

5.2 Local stability with asymmetrical parameterisation from experimental variability

The construction of the I/O dynamics in (5) assumes the symmetrical parameterisation in Definition 4. For a parametric analysis of the I/O system $\dot{\mathbf{p}}=\mathbf{A}_{p}\mathbf{p}+\mathbf{B}_{p}r$ Assumption 1 still holds. Hence, as long as the I/O linear dynamics are stable, Proposition 5.16 guarantees that the nonlinear dynamics are bounded.

Once we (realistically) allow that all the parameters in Example 2.1 can vary independently, we get an asymmetric parameterisation that deviates from Assumption 1. The dynamics for the I/O signals $p_{j}$ are still linear ($\mathbf{W}_{p}\left(\mathbf{Px}\circ\mathbf{x}\right)=0$), however, they depend on the nonlinear dynamics through the term $\mathbf{R}_{12}\mathbf{q}$ (absent in (5) and (38a))

[TABLE]

Remark 5.19**.**

With experimental variability, we lose the serial structure from (38b), and the I/O linear system and the underlying positive dynamics become interconnected (dashed connection in Fig. 2). A stable I/O dynamics $\mathbf{R}_{11}\in\mathcal{H}$ no longer provides guarantees of boundedness, since it ignores the feedback between the I/O linear dynamics and the underlying nonlinear dynamics. Therefore, we need to analyse the stability of the complete nonlinear dynamics of (17).

We investigate the stability of the nonlinear system using Lyapunov’s indirect method, and the eigenvalues of the linearisation at the equilibrium of the system. For an equilibrium $\mathbf{x}=\mathbf{x}^{*}$, $\mathbf{r}=0$, and $\mathbf{J}\left\{\mathbf{x}^{*}\right\}=-\mathbf{D}\left\{\mathbf{P}\mathbf{x}^{*}\right\}-\mathbf{D}\left\{\mathbf{x}^{*}\right\}\mathbf{P}$, the linearisation of (17) results in the following

[TABLE]

If $\alpha\left\{\mathbf{A}_{s}\right\}<0$ then the system is locally exponentially stable around the equilibrium [28]. The equilibrium $\mathbf{x}^{*}=0$ is stable if and only if $\mathbf{A}\in\mathcal{H}$, in agreement with Lemma 5.11. With the participation of $\mathbf{J}\left\{\mathbf{x}^{*}\right\}$, even if $\mathbf{A}$ is not Hurwitz, the linearisation can still be stable around the equilibrium $\mathbf{x}^{*}>0$, showing the stabilising role of the bimolecular reactions. It is also noteworthy that $\mathbf{W}_{p}\mathbf{J}\left\{\mathbf{x}^{*}\right\}=0$, hence $\mathbf{A}_{p}$ and the stability of the linear I/O dynamics does not depend on the equilibrium.

For the particular case of cascaded systems, as long as all species degrade with some non-zero rate, we show in Appendix that the CRN is stable.

6 Analysis of an example nucleic acid feedback control system

To illustrate the application of the above results, we now analyse the feedback system given in Example 2.1. We first consider the nominal parameterisation in Table 1, and analyse the dynamics in the natural coordinates $x^{\pm}_{j}$ in (17) and the I/O linear dynamics from (5). For simulation we assume that the reference signal is a sequence of steps, where only one of the concentrations $r^{+}>0$ or $r^{-}>0$ at any given time. The response with the nominal parameterisation is shown in Fig. 4 where $\mathbf{p}$ and $\mathbf{q}$ are recovered with (36). The output $y=p_{5}$ tracks successfully the reference $r$, while $\mathbf{q}\geq 0$ reveals the underlying dynamics. Since $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$ the origin is unstable (Lemma 5.14), and for $t>7\times 10^{4}\text{ s}$, when the reference returns to $r^{\pm}=0$, the state converges to a positive equilibrium $\bar{\mathbf{x}}^{+*}=\bar{\mathbf{x}}^{-*}>0$.

Table 2 shows that the nominal $\bar{\mathbf{R}}_{11}$ and the linearisation around the nominal equilibrium $\bar{\mathbf{A}}_{s}$ are Hurwitz. However, in reality, experimental variability in the reaction rates leads to asymmetric parameterisations, and the stability of I/O dynamics does not guarantee stability of the CRN. To account for realistic levels of experimental variability, we introduced an uncertainty of $\pm 33\%$ in the reaction rates, which includes the asymmetrical parameterisation shown in Table 1. Perturbing the unforced nonlinear dynamics for this case around its equilibrium ($\mathbf{r}=0$), results in the unstable response of Fig. 5.

The poles in Table 2 show that the linearisation with the asymmetrical parameterisation $\mathbf{A}_{s}$ captures the instability in a pair of conjugated poles on the right-hand plane, despite the stability of the I/O linear system $\mathbf{R}_{11}\in\mathcal{H}$. Indeed, integrating the rotated dynamics with a decoupled matrix $\mathbf{R}$ where we force $\mathbf{R}_{21}=\mathbf{R}_{12}=0$, we obtain the response of Fig. 6, where both ${\mathbf{p}}$ and ${\mathbf{q}}$ have bounded trajectories. This shows that the source of the instability of the complete nonlinear system is neither $\dot{\mathbf{p}}$ nor $\dot{\mathbf{q}}$ individually, and stability must be analysed for the complete interconnected dynamics.

7 Stability of the controller implementation with DSD reactions

It remains to verify whether the stability properties of Example 2.1 predicted from analysing the system CRNs are observed when the closed-loop system is implemented with nucleic acids. In a DSD reaction, a strand of DNA displaces another strand from its binding to a complementary strand, in a random thermodynamic process, which decreases the Gibbs free energy. The single-stranded overhangs, or toeholds, provide initial binding sites for incoming strands to initiate a toehold-mediated branch migration process that can result in strand displacement [10, 29]. Tuning the affinities of the toeholds, based on the base-pair affinities and the nucleotides sequences [29], allows the mapping of the desired reaction rates for the CRN into the DSD implementation.

Following [5], the chemical reactions result in bimolecular DSD reactions, which produce waste in the form of inactivated double strands of DNA which cannot participate in any reaction. Auxiliary fuel species are consumed irreversibly as fuel, and the reactions stop if these are not replenished (details in [5]). The fuel species are initialised at a high concentration $C_{max}=10^{4}\,\text{nM}$, to prevent their consumption from impacting the dynamics significantly.

The DNA strand displacement reactions are simulated using VisualDSD, a rapid-prototyping tool that allows precise analysis of computational devices implemented using DNA strand displacement reactions [11]. The translation of the CRN system follows the construction proposed in [5], with a two-domain programmming structure [30]. The model was parameterised with the unstable parameterisation of Table 1, applying the correspondence between the reaction rates in the CRN and the DSD implementation: $q_{max}=2\eta$, $q_{ki}^{\pm}=2k_{i}^{\pm}/C_{max}$, $i\in\{0,1,2\}$, and $q_{i}^{\pm}=2\gamma_{i}^{\pm}/C_{max}$, $i\in\{1,\dots,8\}$.

The behaviour of the nominal symmetrical parameterisation is first verified in Fig. 7, where $p_{5}$ tracks the step inputs of $r$. After $6\times 10^{6}$ s, the system converges to the positive equilibrium, and Fig. 8 shows the concentrations of the auxiliary species involved in the annihilation reactions remain around $C_{max}$ but are still depleted when $r^{\pm}=p_{j}=0$. Fig. 9 shows that the parameterisation which destabilises the CRN also destabilises the DSD implementation, emphasising the practical relevance of the stability results.

8 Conclusions

Several recent works have applied the dual-rail representation of CRN’s to obtain linear I/O models of synthetic feedback control systems, but have not explicitly considered the potential impact of the underlying nonlinear annihilation reactions in their analysis. This new class of IPR derived from CRNs relies on internally nonlinear positive dynamics. We decomposed the dynamics of the CRN’s involved in a typical linear controller design, and highlighted the effects of the non-observable and nonlinear dynamics - in particular, we showed that the stability of these I/O models does not imply the stability of the underlying chemical network. Under inevitable experimental variability, stability can be affected by the looped interconnection between the nonlinear dynamics arising from biochemical implementation and the linear I/O dynamics resulting from the controller designs. We presented an example of this phenomenon, where the I/O linear system does not capture the instability of the full nonlinear system, and verified this result via simulation of the DSD network that would be implemented experimentally. Our results confirm that the stability of nucleic acid-based controllers must be analysed using the linearisation of the complete nonlinear system, and provide a rigorous theoretical approach for conducting such an analysis.

Acknowledgements

DGB acknowledges funding from the University of Warwick, the EPSRC/BBSRC Centre for Doctoral Training in Synthetic Biology via grant EP/L016494/1 and the BBSRC/EPSRC Warwick Integrative Synthetic Biology Centre via grant BB/M017982/1.

Appendix A The representation of cascaded dynamics has a single equilibrium

For cascaded systems $\bar{\mathbf{R}}_{22}\in\mathcal{L}$, and we can apply Lemma A.20 to conclude that the unforced equilibrium of the cascaded dynamics $\dot{\mathbf{q}}=\bar{\mathbf{R}}_{22}\mathbf{q}-\frac{\eta}{2}\mathbf{q}\circ\mathbf{q}$ is unique at $\mathbf{q}=0$.

Lemma A.20**.**

If $\mathbf{M}\in\mathcal{H}$ is also a lower triangular matrix $\mathbf{M}\in\mathcal{L}$, then $\dot{\mathbf{q}}=\mathbf{Mq}-k\mathbf{q}\circ\mathbf{q}$ has a single equilibrium $\mathbf{q}=0$.

{pf}

Given $\mathcal{M}\in\mathcal{L}$, the solution for $q_{j}$ depends only on $q_{i},\,i\leq j$, with

[TABLE]

For $j=1$, the solutions are $q_{1}\in\{0,m_{11}\}$. $\mathcal{M}\in\mathcal{H},\mathcal{L}\Rightarrow m_{11}\leq 0$ and the only non-negative solution is $q_{1}=0$. Solving sequentially for the remaining coordinates $j=2,3,\ldots,N$, knowing that $q_{i<j}=0$, we have that $\sum_{i=1}^{j-1}m_{ji}q_{i}=0\Rightarrow q_{j}\in\{0,m_{jj}\}$. Since $m_{jj}\leq 0$, the non-negative solution is always $q_{j}=0$.∎

Appendix B Stability analysis and positive equilibrium conditions for Example 4.8

Following on from Remark 4.9, we see how $c_{2}$ impacts differently the spectral radius of $\bar{\mathbf{R}}_{11}$ and $\bar{\mathbf{R}}_{22}$. The characteristic polynomial for the I/O dynamics

[TABLE]

shows that the I/O system is stable for any $c_{2}>0$. On the other hand

[TABLE]

and a gain $c_{2}>{d_{2}d_{1}}{c_{1}^{-1}}$ which stabilises the linear I/O dynamics leads to $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$. Furthermore, the domain for which $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$ and $\mathbf{q}^{*}>0$ exists is the same: $c_{2}>{d_{2}d_{1}}{c_{1}^{-1}}$.

Appendix C Representation of the stability bounds

Fig. 10 illustrates the stability results from Lemmas 5.14 and 5.15.

Appendix D Stability of the CRN representation for a cascaded system, under parameter variability

We can easily state a stability condition for the representation of a cascaded system, even if experimental variability results in an asymmetrical parameterisation.

Proposition D.21**.**

Take the representation of a stable cascaded system $\dot{\mathbf{p}}=\bar{\mathbf{R}}_{11}\mathbf{p},\bar{\mathbf{R}}_{11}\in\mathcal{L},\mathcal{H}$. For an asymmetrical parameterisation (without Assumption 1), if $\mathbf{a}_{1}^{\pm}<0$, the unforced dynamics $\dot{\mathbf{x}}=\mathbf{Ax}-\eta\mathbf{x}\circ\left(\mathbf{Px}\right)$ are GAS for $\mathbf{x}=0$.

{pf}

Given a cascaded I/O dynamics, then we can permutate the state $\mathbf{p}$ so that $\bar{\mathbf{R}}_{11}\in\mathcal{L}$, resulting also $\bar{\mathbf{A}}_{1},\bar{\mathbf{A}}_{2}\in\mathcal{L}$. In the presence of variability, $\mathbf{A_{1}^{\pm}}$ have the same structure as $\bar{\mathbf{A}}_{1}$ but with different parameterisations, resulting $\mathbf{A_{1}^{\pm}}\in\mathcal{L}$. In the same way, $\mathbf{A}_{2}^{\pm}\in\mathcal{L}$. Now take the permutation matrix $\mathbf{Q}$

[TABLE]

such that

[TABLE]

The dynamics of the permuted state result

[TABLE]

where $\mathbf{z}>0\Rightarrow\mathbf{g}\left\{\mathbf{z}\right\}<0$ and

[TABLE]

The structures of $\mathbf{L}_{jj}$ are determined by the structure of $\mathbf{A}_{1}^{\pm}$, and $\mathbf{L}_{21}\geq 0$ contains the cross terms which result in subtractions in the I/O dynamics (elements in $\mathbf{A}_{2}^{\pm}$). Since $\mathbf{A_{1}^{\pm}},\mathbf{A_{2}^{\pm}}\in\mathcal{L}$, $\mathbf{L}$ results triangular, Moreover, $\mathbf{L}=\mathbf{L}^{\boxbslash}+\mathbf{D}\left\{\mathbf{l}\right\}$ where \mathbf{l}=\mathbf{Q}\left[\begin{array}[]{cc}\mathbf{a}_{1}^{+T}&\mathbf{a}_{1}^{-T}\end{array}\right].

It results directly that $\rho\left\{\mathbf{A}\right\}=\rho\left\{\mathbf{L}\right\}=\rho\left\{\mathbf{A}_{1}^{+}\right\}\cup\rho\left\{\mathbf{A}_{1}^{-}\right\}$, and $\mathbf{a}_{1}^{\pm}<0\Leftrightarrow\mathbf{L}\in\mathcal{H}\Leftrightarrow\mathbf{A}\in\mathcal{H}$. If the represented cascaded linear I/O dynamics are stable, then $\mathbf{A}\in\mathcal{H}$. Moreover, even with uncertainty, as long as the degradation rates remain positive $\mathbf{A}\in\mathcal{H}$. Since $\mathbf{A}\in\mathcal{H}$, we can invoke Lemma 5.11 to establish $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}-\eta\left(\mathbf{Px}\right)\circ\mathbf{x}$ is GAS around $\mathbf{x}=0$. ∎

In the presence of variability we can have mismatching rates, but as long as all species degrade with some non-zero rate, the unforced dynamics of the cascaded system will have a single stable nonnegative equilibrium. Without input, the CRNs will converge to rest at $\mathbf{x}=0$. Table 3 summarises the derived properties, depending on the structure of the DSD network (cascaded versus with feedback).