A Theoretical Framework for Implementable Nucleic Acids Feedback Systems

Abstract

Chemical reaction networks can be utilised as basic components for nucleic acid feedback control systems’ design for Synthetic Biology application. DNA hybridisation and programmed strand-displacement reactions are effective primitives for implementation. However, the experimental validation and scale-up of nucleic acid control systems are still considerably falling behind their theoretical designs. To aid with the progress heading into experimental implementations, we provide here chemical reaction networks that represent two fundamental classes of linear controllers: integral and static negative state feedback. We reduced the complexity of the networks by finding designs with fewer reactions and chemical species, to take account of the limits of current experimental capabilities and mitigate issues pertaining to crosstalk and leakage, along with toehold sequence design. The supplied control circuits are quintessential candidates for the first experimental validations of nucleic acid controllers, since they have a number of parameters, species, and reactions small enough for viable experimentation with current technical capabilities, but still represent challenging feedback control systems. They are also well suited to further theoretical analysis to verify results on the stability, performance, and robustness of this important new class of control systems.

1. Introduction

With the recent advancement of cost-effective approaches for assembling oligonucleotides, synthetic DNA has become an appealing preference for constructing circuits using hybridisation and DNA strand displacement (DSD) reactions. Programmable and biologically compatible synthetic circuits that utilise DSD reactions have been shown to function effectively in living organisms (in vivo), as reported in a study by [1]. Therefore, these DSD networks are a promising candidate for building synthetic feedback controllers in biomolecular settings [2], to drive DNA-based molecular devices, or for the realisation of biochemical computing circuits with emerging opportunities and challenges also for applications in vivo [3].

Systematic methods to compose and program DNA circuits have been proposed [4] that enable the translation of chemical programs to reactions through the use of nucleic acids, where the binding affinities and the reaction rates can, to a certain degree, be controlled [5,6]. Feedback control systems can be assembled with elementary chemical reaction networks (CRNs) that provide an abstract layer for designing mathematical operators [7], prior to being translated to biomolecular applications. As a result, the networks of fundamental reactions such as catalysis, degradation, and annihilation can then be mapped to their equivalent DSD reactions [8], where the DNA sequences effectively program the biochemical interactions and each of the strand displacement reaction rates [5,9].

Nonetheless, at first glance, the concentrations of biomolecular species are seemingly ill-suited to represent signals in the theory of feedback control, given that they are limited to positive quantities. For instance, in a typical reference tracking control problem, the error between the reference and the plant output can take either a negative or a positive value, but then, this cannot be represented by a positive concentration. An alternative technique that works around this limitation is to employ the so-called dual-rail representation, in which pairs of concentrations are used to represent system state variables [10]. Specifically, a signal $x\left(t\right)$ is described by the difference between two concentrations $x\left(t\right)=x^{+}\left(t\right)-x^{-}\left(t\right)$, which corresponds to two chemical species $X^{+}$ and $X^{-}$. By making such an abstraction, this approach enables the representation of gains and negative signals using positive quantities that are essential for error computation in a negative feedback control system. The dual-rail representation and implementations for elementary catalysis, degradation and annihilation reactions provide a systematic framework to convert control theory into biochemistry, which can be implemented with synthetic DNA oligonucleotides, where the representation of transfer functions [11], linear feedback systems [12,13,14,15,16,17], and nonlinear controllers [18,19,20] can be realised through the assembly of DSD networks.

While open-loop cascades for analogue or logical computations implemented with DSD networks have been constructed [21], there remains a considerable gap between theory and experiment when it comes to dynamical circuits that implement negative feedback. Moreover, even the most basic types of feedback controllers have yet to be implemented and reported. This is an important limitation to further progress, because the proposed dual-rail representation requires annihilation reactions, which are essential to ensure that species concentrations remain within the bounds of experimental feasibility, but lead to nonlinear dynamics due to the bimolecular nature of annihilation reactions.

Ideally, these additional nonlinear dynamics are not observable in the represented input to output linear dynamics, but still determine some important properties of the CRN, which need careful consideration for experimental implementation, such as the level of concentrations at which the circuit operates and unforced positive equilibria with persistent consumption of chemical species [22,23]. Additionally, the implementation of the DSD reactions introduce leakage and spurious effects [24], which prompted study on leakage analysis and ways to mitigate it [25]. Furthermore, the variability and granular design of hybridisation rates [5] introduce error and mismatches in the reaction rates. All these aforementioned issues show that we are lacking an experimental validation of control systems based on DNA hybridisation, to help us understand how nonlinearities and uncertainty in the implementation impact the stability and performance of these control circuits.

Due to several technical challenges, which include cross-talk, designing toehold sequences, leakage, etc., there is a substantial lag in the construction, validation, and scale-up of nucleic acid control systems compared to the theoretical advancements’ counterpart [5,24]. While the nucleotide sequence can be predicted from the toehold affinities [5], the kinetics can be affected by undesirable binding scenarios, resulting in the modification of the reaction rate constants and the dynamics (e.g., in oscillating circuits [24] or the seesaw gate circuit [26]). The activation of these undesirable reactions that yields the leakage of the outputs in the absence of the input also poses another major experimental hurdles. A potential approach in addressing this issue is the use of clamps, which could impede the spurious hybridisations [24,27], or compartmentalisation and localisation strategies [28,29,30], which isolate strand complexes, which could cause leak reactions. Additionally, the severity of leakage increases when concentrations are high; thus, it is a common practice to ensure low concentrations of the reacting species [31]. This leakage at high concentrations along with the restrictions on hybridisation rates [5] impose upper bounds on the speed at which these circuits can operate. A technique that can be used to mitigate this is localisation, whereby the adjacent gates are placed for interaction without diffusion at faster rates [26]. Lastly, despite measures taken to prevent these spurious reactions, it is still essential that the sequence of reactions, which are competing for common reactants, are well-managed through the means of timescale separation or compartmentalisation [32]. With larger and more complex circuits, sequence similarity aggravates all these undesired effects, with a stronger need for toehold-mediated strand displacement circuits, which are robust to molecular noise [33].

All the experimental challenges described above can be addressed through the designs of circuits that have a lesser amount of reactions. With a lesser amount of reactions and species, this reduces the number of designs required for the template strands. This concurrently reduces the demand for the characterisation and tuning of the kinetic rates, lessens the sources of error, as well as curtails the possibility of leakage and undesired interactions. Ensuring fewer reactions and species are of pertinence, especially in using the dual-rail representation, which necessitates a duplication of nearly all of the reactions [19].

A key challenge from a theoretical perspective is, therefore, to design circuits of reduced complexity, using the least amount of chemical reactions, at the same time still able to represent a negative feedback control problem of interest in an accurate manner. All these are done to ensure that, with the presently available technical capacity, the likelihood of implementing a successful experiment can be maximised. Here, we suggest chemical representations for two fundamental types of feedback control systems in order to achieve the above objective. Our first type of circuit is a reference-tracking problem employing integral control for a stable first-order system, with a single tuning parameter. The second type of circuit is a combined reference tracking and stabilisation problem of a marginally stable plant, i.e., the classic double-integrator system employing the static state feedback control of the two integrated states with two tunable gains. The dual-rail framework is used in both of these circuits, and their parametrisations can be tuned to satisfy control requirements such as steady-state tracking and transient dynamics. The supplied circuits of reduced complexity are ideal candidates for the first experimental validation and robustness assessment of nucleic acid controllers implemented via the dual-rail representation.

An early version of this study can be found in [34], where we established the fundamental theories of chemical reaction network that are required to realise the reference tracking control problem. Here, we extended the work in [34] by applying the theories to design two types of controllers and investigated the impact of the nonlinearities and uncertainties on the experimental implementation, in particular the effect of leakage, which was not considered previously.

The remainder of this article is organised as follows. Section 2 introduces the definitions and methods used in this work. The formulations and constructions with DNA hybridisation reactions are detailed for control systems with integral action in Section 3 and with state feedback in Section 4. Section 5 discusses the impact of nonlinearities and uncertainty in performance and robustness, highlighting the properties of interest for the experimental investigation. Section 6 summarises the main conclusions.

2. Definitions and Preliminaries

2.1. Chemical Reaction Networks

The theory of CRNs has its roots in chemistry [35], whereby the same principle is applied in a similar manner to biomolecules to model and analyse biochemical reaction networks. With their extensive computational capabilities [36,37], CRNs provide a convenient representation for implementing elementary arithmetic operations [7] or the computation of polynomials [38], using any chemical system with mass action kinetics. They also provide an appropriate level of abstraction for designing complex circuits [39] and integrate the different elements necessary to build linear feedback control systems [12,13].

Definition 1.

Given a set of M reactions between N chemical species $X_j$, a chemical reaction network is represented by

$$\begin{array}{c}\hfill \sum _{j=1}^N{a}_{jm}{X}_j\stackrel{{\gamma }_m}{\to }\sum _{j=1}^N{b}_{jm}{X}_j,\left(m=1,\dots ,M\right)\end{array}$$

(1)

where the stoichiometric coefficients $a_{jm}$ and $b_{jm}$ determine, for each reaction m, the quantities of the reactants consumed on the left and the products output on the right, at a rate ${\gamma }_m$. When, due to degradation or sequestration, a product is removed or becomes nonreactive, it is omitted or replaced by the symbol Ø.

When the same species is utilised in multiple reactions, it creates links between different species that can be expressed as a graph to describe the network. There are various sources that provide the definitions of the models and properties pertaining to CRNs. These include Petri nets and monotone systems [40,41]; deficiency theorems and structural analysis work [42,43,44,45]; and biochemical networks’ application [46].

Different principles have been used to model the time evolution of the concentration of species [35]. In this study, the main principle that we used is the Mass-Action-Kinetics (MAK)-based deterministic models, which have been traditionally used to model biochemical processes [47,48,49], as well as used in the implementation of polynomial Ordinary Differential Equations (ODEs) and Turing universal CRNs [36,50]. Assuming the MAK, a set of ODE can be used to express a deterministic model [51,52,53] that is appropriate for the analysis of control and systems theory, despite being nonlinear, with the reaction fluxes depending on the product of the reatanct’s concentrations.

Definition 2.

For a deterministic model based on the MAK [35], we have that, for each species $X_j$ in (1), the dynamics of its concentration $x_j$ is given by

$$\begin{array}{c}\hfill {\dot{x}}_j=\sum _{m=1}^M\left({\gamma }_m\left(b_{jm}-a_{jm}\right)\prod _{j=1}^N{x}_j^{a_{jm}}\right)\end{array}$$

(2)

Example 1.

Take the reaction between the reactant species $X_1$ and $X_2$, which produces species $X_3$, represented by

$$\begin{array}{c}\hfill a_1{X}_1+a_2{X}_2\stackrel{\gamma }{\to }b{X}_3\end{array}$$

(3)

where the reactants on the left are converted into the product on the right at a rate γ, according to the stoichiometric coefficients $a_1$, $a_2$, and b. The MAK model results in the ODEs:

$${\dot{x}}_1=-a_1\gamma x_1^{a_1}{x}_2^{a_2}$$

(4)

$${\dot{x}}_2=-a_2\gamma x_1^{a_1}{x}_2^{a_2}$$

(5)

$${\dot{x}}_3=+b\gamma x_1^{a_1}{x}_2^{a_2}$$

(6)

where the stoichiometric coefficients $a_1$, $a_2$, and b indicate, respectively, the relative number of molecules consumed and produced during the reaction.

Let ${\mathbb{R}}_0^{+}$ denote the positive orthant, where $v_j\in {\mathbb{R}}_0^{+}\Rightarrow v_j\ge 0$; the dynamics of (4) to (6) are expressed in their natural coordinates, and the concentrations $x_j\in {\mathbb{R}}_0^{+}$ are always non-negative. The system is non-negative (${\forall }_j,x_j\left(0\right)\ge 0\Rightarrow x_j\left(t\right)\ge 0$) if the sets of reactions are modelled properly using realistic non-negative initial conditions.

2.2. Dual-Rail Representation

In the context of feedback control, a key challenge in employing CRNs is their inability to directly represent negative signals, since the concentrations of chemical species are always positive. In a linear negative feedback system, we cannot restrict ourselves to considering only a positive difference between two positive inputs, i.e., “one-sided” subtraction [54]. The limitation of the positivity of CRNs can be circumvented using the dual-rail representation [10,12,55], where this enables the positive and, in particular, the negative signals that are vital for error signals’ computation in a feedback control to be represented using molecular concentrations.

Definition 3.

Consider two chemical species $X_j^{+}$ and $X_j^{-}$ and respective concentrations $x_j^{+}\ge 0$ and $x_j^{-}\ge 0$. A dual-rail signal $p_j\in \mathbb{R}$ is represented by $p_j=x_j^{+}-x_j^{-}$, with dynamics given by ${\dot{p}}_j={\dot{x}}_j^{+}-{\dot{x}}_j^{-}$.

Here, we have signals that are generated by subtracting two positive quantities, which enables us to obtain positive or negative results when subtracting signals instead of concentrations.

Definition 4.

The Input–Output (I/O) system is the response $Y\left(s\right)=G\left(s\right)U\left(s\right)$, from an input $u=\left(u^{+}-u^{-}\right)$ to an output $y=\left(y^{+}-y^{-}\right)$. The states are also dual-rail $p_j=x_j^{+}-x_j^{-}$, where $u,y,p_j\in \mathbb{R}$ and $u^{\pm },y^{\pm },x_j^{\pm }\in {\mathbb{R}}_0^{+}$.

Combining the MAK and the dual-rail representation, we can model linear systems with both positive and negative signals, using CRNs composed of only the three types of reactions:

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

(7)

$$X_i\stackrel{\gamma }{\to }Ø$$

(8)

$$X_i+X_j\stackrel{\eta }{\to }Ø$$

(9)

where we have unimolecular reactions for catalysis (7) and degradation (8) and a bimolecular reaction in annihilation (9).

Example 2

(First-order system). To chemically represent the first-order system with a negative gain given by the transfer function:

$$\begin{array}{c}\hfill Y\left(s\right)=-\frac{k_1}{s+k_2}U\left(s\right)\end{array}$$

(10)

with $u,y\in \mathbb{R}$, we take the pairs of chemical species $\left\{U^{\pm },Y^{\pm }\right\}$ and the CRN from Figure 1 with

$$U^{\pm }\stackrel{k_1}{\to }{U}^{\pm }+Y^{\mp }$$

(11)

$$Y^{\pm }\stackrel{k_2}{\to }Ø$$

(12)

$$Y^{+}+Y^{-}\stackrel{\eta }{\to }Ø$$

(13)

From the MAK, we derive the ODEs for the concentrations:

$${\dot{y}}^{+}=-k_2{y}^{+}+k_1{u}^{-}-\eta y^{+}{y}^{-}$$

(14)

$${\dot{y}}^{-}=-k_2{y}^{-}+k_1{u}^{+}-\eta y^{+}{y}^{-}$$

(15)

From Definition 4, the linear I/O system results:

$$\begin{array}{ccc}\hfill {\dot{y}}^{+}-{\dot{y}}^{-}& =& -k_2\left(y^{+}-y^{-}\right)-k_1\left(u^{+}-u^{-}\right)\iff \dot{y}=-k_{2}y-k_{1}u\hfill \end{array}$$

(16)

The computation of a two-sided subtraction with the steady-state of a CRN or the representation of negative gains is enabled by the dual-rail representation [55]. Although this increases the number of required reactions (by duplicating the unimolecular reactions), it does provide a systematic framework with which to represent linear and negative feedback systems [11,12,14].

2.3. DNA Strand Displacement Reactions

Sequences of four types of nucleotides, i.e., Adenine (A), Guanine (G), Cytosine (C), and Thymine (T), are formed by chains of nucleotides in nucleic acids. The formation of hydrogen bonds between nucleotide pairs A-T and C-G results in a double-stranded DNA with a helical shape. As depicted in Figure 2, this DNA structure is composed of two strands of DNA that run antiparallel to each other. The bonding between the complementary pairs of A-T and C-G is stabilised by the hydrophilic backbone and the hydrophobic nature of the bases. As illustrated in Figure 2b, this results in the formation of double-stranded DNA through an enzyme-free hybridisation reaction between two antiparallel complementary strands of DNA [56].

The dynamics of simple chemical reactions such as catalysis, degradation, and annihilation can be readily translated into nucleic-acid-based reactions with comparable behaviours [8,9]. This approach enables technology for implementing circuits designed using CRNs, where the DNA sequence of the nucleotides within the strands can efficaciously be used to program the biochemical circuitry to realise both analogue and digital computations [4,5,26,54]. Each of these reactions has equivalent representations with DSD reactions, where formal species and reactions in the CRN are mapped to sets of DNA species. Various suggested architectures [8,9,24,32] permit the creation of circuits that are aided by a high-level of automation through the utilisation of existing syntax and software tools [4,13,57].

The DSD reaction uses DNA molecules as signal species, and they consist of binding domains and single-stranded overhanging and exposed toeholds. The function of these toeholds is to initiate and adjust the reaction rates of strand hybridisation [4,5], as depicted in Figure 2c. The hybridisation process leads to bimolecular reactions with unitary stoichiometric coefficients. Using the following DSD reaction of $A+B\to C+D$ as an illustration (Figure 3), the toehold 1 of the approaching strand A is hybridised to its complementary toehold $1^{*}$ in B, which leads to the commencement of branch migration. This then displaces Domain 2 and releases the two output species of C and D. From the thermodynamic point of view, the unbinding of toehold 1 is less favoured compared to the displacement of the incumbent domain 2, resulting in C being completely displaced. With the two output species C and D having no overhanging complementary toeholds for further hybridisation, the displacement mechanism is deemed not reversible, assuming leakage reaction rates are much slower than the displacement mediated by the toehold reactions. The reverse reaction, where C displaces the strand A from D, can be up to orders of magnitude slower [58] and be considered effectively unreactive.

With the assumption of an irreversible reaction coupled with the output strand C being able to take part in other reactions, the cascade of several reactions ensues. It is often a challenge to describe the hybridisation reaction with good accuracy. Nonetheless, there are several predictive models available that utilise the nucleotide sequences to provide a good approximation of the reaction rate [5]. Using Visual DSD [57], which is a specifically designed computer-aided design software for strand displacement reactions’ analysis, we conducted a simulation to analyse and verify the DSD circuity by programming the species and their affinities in the software. We ran the simulation in “default” mode, which is the mode that considers finite binding/unbinding rates. The affinities between the toeholds, which facilitate the initial hybridisation between the strands, are determined based on the reaction rates in the CRNs.

3. Integral Feedback Control System

We now propose a CRN with reduced complexity to represent the integral feedback control system in Figure 4. This is then followed by implementing it in Visual DSD using strand displacement reactions.

Example 3

(Reference tracking problem using integral control). Define the first-order plant and the control law:

$$\begin{array}{ccc}\hfill Y\left(s\right)& =& \frac{b}{s+a}V\left(s\right)\Rightarrow sY\left(s\right)=bV\left(s\right)-aY\left(s\right)\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill V\left(s\right)& =& \frac{k_i}{s}\left(R\left(s\right)-Y\left(s\right)\right)\hfill \end{array}$$

(18)

where $V\left(s\right)$ is the integral control action to have the output $Y\left(s\right)$ track a reference $R\left(s\right)$, according to Figure 4.

The transfer function for the closed-loop dynamics results:

$$\begin{array}{c}\hfill Y\left(s\right)=\frac{b{k}_i}{s^2+as+b{k}_i}R\left(s\right)\end{array}$$

(19)

where the integral control ensures steady-state tracking with $Y\left(0\right)=R\left(0\right)$. The transient dynamics are defined by the roots of the characteristic polynomial:

$$\begin{array}{c}\hfill \lambda =\frac{1}{2}\left(-a\pm \sqrt{a^2-4b{k}_i}\right)\end{array}$$

(20)

with a natural frequency ${\omega }_n=\sqrt{b{k}_i}$ and damping coefficient $\xi =\frac{a}{2\sqrt{b{k}_i}}$.

3.1. Representation with Chemical Reactions

The closed-loop response represented by a network of catalysis, degradation, and annihilation reactions is shown in Figure 5, where

$$V^{+}\stackrel{b}{\to }{V}^{+}+Y^{+},V^{-}\stackrel{b}{\to }{V}^{-}+Y^{-}$$

(21)

$$Y^{+}\stackrel{a}{\to }Ø,Y^{-}\stackrel{a}{\to }Ø$$

(22)

$$R^{+}\stackrel{k_i}{\to }{R}^{+}+V^{+},R^{-}\stackrel{k_i}{\to }{R}^{-}+V^{-}$$

(23)

$$Y^{+}\stackrel{k_i}{\to }{Y}^{+}+V^{-},Y^{-}\stackrel{k_i}{\to }{Y}^{-}+V^{+}$$

(24)

$$Y^{+}+Y^{-}\stackrel{\eta }{\to }Ø,V^{+}+V^{-}\stackrel{\eta }{\to }Ø$$

(25)

The reactions (21) and (22) describe the plant, where the reaction rates a and b, respectively, set the gain and the degradation reactions, with a stable pole at $s=-a$.

The additional dynamics of the chemical reactions used to solve the subtraction and compute the control error as found in [12] or [14] are removed. Instead of using separate CRNs for subtraction and integration, both operations are combined to reduce the number of reactions and simplify the circuit. Thus, Equations (23) and (24) describe the subtraction of the contributing reference and the output to the integral action through the crossing of the contributions from $Y^{\pm }$ to $V^{\mp }$. The control gain is applied by the reaction rate $k_i$ in the corresponding catalysis reactions. The annihilation reactions in (25) ensure low (i.e., experimentally feasible) concentrations of the chemical species [12,23].

The respective MAK governing the species concentrations can be written as

$${\dot{v}}^{+}=k_i{r}^{+}+k_i{y}^{-}-\eta v^{+}{v}^{-}$$

(26)

$${\dot{v}}^{-}=k_i{r}^{-}+k_i{y}^{+}-\eta v^{+}{v}^{-}$$

(27)

$${\dot{y}}^{+}=b{v}^{+}-a{y}^{+}-\eta y^{+}{y}^{-}$$

(28)

$${\dot{y}}^{-}=b{v}^{-}-a{y}^{-}-\eta y^{+}{y}^{-}$$

(29)

This gives us then that, for the dual-rail variables $y=y^{+}-y^{-}$, $v=v^{+}-v^{-}$, and $r=r^{+}-r^{-}$, the dynamics of the I/O are the outcome of the expressions $\dot{v}={\dot{v}}^{+}-{\dot{v}}^{-}$ and $\dot{y}={\dot{y}}^{+}-{\dot{y}}^{-}$, and this yields

$$\dot{v}=k_i\left(r^{+}-r^{-}\right)+k_i\left(-y^{+}+y^{-}\right)=k_i\left(r-y\right)$$

(30)

$$\dot{y}=bv-ay$$

(31)

The swapped contributions from $Y^{\pm }$ to $V^{\mp }$ still result in positive gains in the dynamics of the concentrations in (26) to (29), but result in a negative gain on the output y when computing the dynamics in (30) to (31).

The dynamics of the input–output given in (30) to (31) are linear, corresponding to the transfer function of the closed-loop system in (19). To further simplify the construction, we disregarded any annihilation reaction among the reference species $R^{+}+R^{-}\stackrel{\eta }{\to }Ø$ with the assumption that we have low and constant input concentrations during the whole circuit operation, and the dynamics in (30) to (31) rely solely on the difference of r instead of the $R^{+}$ and $R^{-}$ concentration levels.

3.1.1. Dual-Rail CRN as a Nonlinear Internal Positive Representation

The dual-rail representation uses a positive system (the MAK of the CRN) to represent a non-positive system. This approach of having an Internally Positive Representation (IPR) [59] lifts the constraint of positivity imposed by the implementation of the physical system using chemical species. The representation uses input and output transformations that translate the input and output of the system into and from positive states of the dynamics, respectively, to realise arbitrary I/O dynamics. In our case, the internal positive dynamics are represented by CRNs (Figure 6) and the positive variables are the natural coordinates of the MAK given by positive concentrations.

When representing linear I/O systems with a linear IPR, we need to ensure that the system is stable from input to output, but also internally stable. This results in limitations when building the internal linear process [59]. However, the dual-rail representation constructed with the elementary reactions of catalysis, degradation, and annihilation reactions results in a nonlinear MAK. The I/O system in (30) to (31) is linear because the dynamics of the differences of concentrations have the nonlinear terms of the MAK in (26) to (29) cancelling out. However, internally, we have a process that is not only positive, but also nonlinear and not observable in the I/O system.

Definition 5

(Rotated coordinates [23]). Define the transformation of coordinates $\overline{v}=v^{+}-v^{-}$ and $\widehat{v}=v^{+}+v^{-}$, where we can recover the concentrations with $v^{\pm }=\frac{1}{2}\left(\widehat{v}\pm \overline{v}\right)$.

In the new coordinates, we can clearly express the relationship between the I/O system and the internal positive nonlinear dynamics [23]. Since $\eta y^{+}{y}^{-}=\frac{\eta }{4}\left({\widehat{y}}^2-{\overline{y}}^2\right)$, we can express the MAK of Example 3 in the new coordinates, where we obtain

$$\dot{\overline{v}}=k_i\left(\overline{r}-\overline{y}\right)$$

(32)

$$\dot{\overline{y}}=b\overline{v}-a\overline{y}$$

(33)

$$\dot{\widehat{v}}=k_i\left(\widehat{r}+\widehat{y}\right)-\frac{\eta }{2}{\widehat{v}}^2+\frac{\eta }{2}{\overline{v}}^2$$

(34)

$$\dot{\widehat{y}}=b\widehat{v}+a\widehat{y}-\frac{\eta }{2}{\widehat{y}}^2+\frac{\eta }{2}{\overline{y}}^2$$

(35)

The dynamics of $\dot{\overline{v}}$ and $\dot{\overline{y}}$ correspond to the I/O system and $\dot{\widehat{v}}$ and $\dot{\widehat{y}}$ are the remaining dynamics in the dual-rail representation, which are nonlinear and depend on the trajectories of the I/O system. It is noteworthy that the I/O system does not rely on $\dot{\widehat{v}}$ or $\dot{\widehat{y}}$, and these underlying nonlinear dynamics are, by design, not observable in the represented linear system.

Assumption 1.

It was assumed that the designed I/O system results in a stable closed feedback system.

Assuming that the controller is designed to ensure stability for the feedback system, we have that the I/O system is a bounded and stable input into the internal nonlinear dynamics. Assumption 1 together with previous theoretical results from [23] ensure that, although nonlinear and partially unobservable, the dynamics in (32) to (35) remain stable.

3.2. Representation with Strand Displacement Reactions

In order for CRNs to be translated into nucleic acid chemistry, the DSD reactions are set with comparable dynamics to each of the three elementary reactions, coupled with toehold mediation to tune their reaction rates. Following [9], both catalysis and degradation reactions can be implemented using the Join–Fork templates. A signal strand can be released through the interaction of a $Join$ double-stranded complex with either one or a few single-stranded inputs. This then leads to a sequence of DSD reactions in the $Fork$ template, where one or more output strands are released. The auxiliary templates, namely the $Join$ and $Fork$ templates, as well as the auxiliary two-domain strands, need to be present in high concentrations and supplied at the start of the reaction.

Figure 7 shows an example of sets of DSD reactions equivalent to the unimolecular reaction of catalysis $R\stackrel{}{\to }R+V$, where the interactions with the intermediary strands and auxiliary templates are mediated utilising toeholds. The sets of reactions are

$$R+\mathit{JoinVR}\underset{k_{bnd}}{\overset{c_i\times k_t}{{\displaystyle \rightleftharpoons }}}JoinVR\_1+\mathit{aux}\_\mathit{hrtp}$$

(36)

$$JoinVR\_1+\mathit{aux}\_\mathit{tphpr}\underset{k_{bnd}}{\overset{k_{bnd}}{{\displaystyle \rightleftharpoons }}}JoinVR\_2+sig\_hprtq$$

(37)

$$sig\_hprtq+\mathit{ForkVR}\underset{k_{bnd}}{\overset{k_{bnd}}{{\displaystyle \rightleftharpoons }}}ForkVR\_1+\mathit{aux}\_\mathit{tphpr}$$

(38)

$$ForkVR\_1+\mathit{aux}\_\mathit{hrtp}\underset{k_t}{\overset{k_{bnd}}{{\displaystyle \rightleftharpoons }}}ForkVR\_2+R$$

(39)

$$ForkVR\_2+\mathit{aux}\_\mathit{hvtr}\underset{k_t}{\overset{k_t}{\rightleftharpoons }}ForkVR\_3+V$$

(40)

$$ForkVR\_3+\mathit{aux}\_\mathit{hitv}\stackrel{k_t}{\to }Ø$$

(41)

The species highlighted in bold in (36) to (41) are the template complexes and auxiliary single-stranded species, which are made available at a high concentration $C_{max}$, to avoid their irreversible consumption from significantly impacting the dynamics.

As shown in Figure 7, every signal species is a single-stranded DNA that consists of toehold (<tr>) and binding (<hr>) domains. The hybridisation to the multi-stranded complex $JoinVR$ is initiated by the toehold domain, which activates a series of processes, leading to an intermediary strand $sig\_hprtq$ being released. The presence of this new intermediary strand activates another series of strand displacements, but this time involving the $ForkVR$ complex. This results in the signal species V being released and a strand of R being returned (following the stoichiometry of $R\stackrel{}{\to }R+V$).

Remark 1.

The reversible reactions in (37) and (38) assume the same forward and backward binding rate $k_{bnd}$. Implicit in this assumption is that the sequences of the toeholds <tp> and <tp*> can be designed to have the same maximum binding rate as the toeholds <tq> and <tq*>. The same assumption is made in (40) for <tr>, <tr*>, <tv>, and <tv*>.

As depicted in Figure 8, the DSD reactions for the degradation reaction are less complicated given that the $JoinY$ complex only requires capturing the signal species Y in an irreversible manner, with

$$\begin{array}{c}\hfill \mathit{JoinY}+Y\stackrel{k_t\times c_a}{\to }JoinY\_1+dY\end{array}$$

(42)

Considering that there are no exposed toeholds in the product of the double-strand, this product becomes inactive and participates no further in the reactions.

To implement the annihilation reactions, we propose the cooperative hybridisation from [61] and applied in [32]. An auxiliary species $AnnVV$ at high concentrations quickly hybridises with both input signals to generate waste species, and in order for the cascade to become irreversible, the two input species need to be present. In the ensuing cooperative hybridisation mechanism, a single template per signal is then used to set the annihilation reactions $V+V^{\prime }\stackrel{}{\to }Ø$, where

$$\begin{array}{c}\hfill V+V^{\prime }+AnnVV\stackrel{c_{VV}\times k_t}{\to }Ø\end{array}$$

(43)

As shown in Figure 9, there is a set of DSD reactions that includes intermediary complexes together with reversible displacement reactions, resulting in the following CRNs:

$$\mathit{AnnVV}+V\underset{k_{ubnd}}{\overset{c_{VV}\times k_t}{{\displaystyle \rightleftharpoons }}}I\mathrm{v}$$

(44)

$$I\mathrm{v}+V^{\prime }\stackrel{c_{VV}\times k_t}{\to }W\mathrm{v}+W{\mathrm{v}}^{\prime }$$

(45)

$$\mathit{AnnVV}+V^{\prime }\underset{k_{ubnd}}{\overset{c_{VV}\times k_t}{{\displaystyle \rightleftharpoons }}}I{\mathrm{v}}^{\prime }$$

(46)

$$I{\mathrm{v}}^{\prime }+V\stackrel{c_{VV}\times k_t}{\to }W\mathrm{v}+W{\mathrm{v}}^{\prime }$$

(47)

Both V and $V^{\prime }$, when present, facilitate the production of two double-stranded waste complexes $W\mathrm{v}$ and $W{\mathrm{v}}^{\prime }$ from the second irreversible reaction.

The nucleotide sequences in the toeholds, which initiate strand hybridisation, allow for the programmability of the DSD reactions. It is well known that the predictable hybridisation kinetics is defined by the affinities between the base pairs [5] despite several factors that could likewise impact the effective reaction rate constants [24]. To facilitate our programming and analysis using Visual DSD, we adjusted the rates of the DSD reactions through the assignment of complementarity degrees between toeholds [13,57] based on the assumption that the binding affinities are weakened by the possibility of toehold design and nucleotide sequences mismatch. As illustrated in Figure 7, the complex $JoinVR$ has the toehold <tr*ci> that possesses the complementarity degree to the signal toehold <tr> of $0<c_i\le 1$. This slows down the reaction mediated by <tr> and <tr*ci> to $k_t\times c_i$, where $k_t$ is the maximum binding rate of two complementary toeholds <tr> and <tr*>.

Assumption 2.

The auxiliary strands were assumed to be present at high concentrations and to remain constant or approximately close to their initial value during the operation of the circuit (or indefinitely through replenishment).

Assumption 2 allowed us to approximate the bimolecular hybridisation reactions with the unimolecular reactions in (7) and (8). For example, the implementation in (42) can be approximated assuming $\left[JoinY\right]\left(t\right)\approx \left[JoinY\right]\left(0\right)=C_{max}$, resulting in the approximating unimolecular reaction:

$$\begin{array}{c}\hfill Y\stackrel{C_{max}\times c_a\times k_t}{\to }JoinY\_1+dY\end{array}$$

(48)

where the constant and large presence of the auxiliary species $JoinY$ is translated directly into the unimolecular reaction rate [8].

3.3. Initialisation and Modelling of the Integral Control System

The six catalysis and two degradation reactions from the CRN in (21) to (25) with the exception of one annihilation reaction were used in implementing the DSD reactions. In the analysis using Visual DSD, the reaction $Y^{+}+Y^{-}\stackrel{}{\to }Ø$ can be removed to further reduce the implementation given that both the $Y^{+}$ and $Y^{-}$ concentrations are sufficiently low.

The templates shown in Figure 7, Figure 8 and Figure 9 were employed to set the DSD reactions. A single $Fork$ template was used for both the outputs of the catalysis reaction, where through cooperative hybridisation, the identical sequestering template is shared between the two annihilated species. A high concentration of 15 double-stranded and 20 single-stranded DNA templates and auxiliary species, respectively, was used for circuit initialisation. The domains for signal species V, V′, Y, Y′, and R, R′ that, respectively, represent $V^{+}$, $V^{-}$, $Y^{+}$, $Y^{-}$, and $R^{+}$, $R^{-}$ are shown in Figure 10. The same figure also shows the provided auxiliary templates with the details of their toeholds and domains that were constructed for hybridisation and strand displacement.

In the Visual DSD simulation, we set $C_{max}={10}^4$ nM (the consumption of the auxiliary strands are not reversible and replenished). We also set the maximum binding and unbinding rate of the toehold to ${10}^{-3}{\left(\mathrm{nMs}\right)}^{-1}$ and $0.1{\mathrm{s}}^{-1}$, respectively.

Table 1. The integral control circuit’s parametrisation for Visual DSD simulation.

Variable	Description	Values	Units
$C_{max}$	Template and auxiliary species initial concentrations	${10}^4$	nM
$k_{bnd}$	Auxiliary strands’ (<tp>,<tq>) toehold maximum binding rate	${10}^{-3}$	${\left(\mathrm{nMs}\right)}^{-1}$
$k_{ubnd}$	Unbinding rate	$0.1$	${\left(\mathrm{s}\right)}^{-1}$
$k_t$	Signal species’ (<tr>,<ty>,<tv>) toehold maximum binding rate	${10}^{-4}$	${\left(\mathrm{nMs}\right)}^{-1}$
$c_a$	Degradation of $Y^{\pm }$ toehold complementarity degree	$2.5\times {10}^{-3}$	−
$c_b$	Plant catalysis reaction toehold complementarity degree	${10}^{-3}$	−
$c_i$	Feedback catalysis implementing integral control toehold complementarity degree	$5\times {10}^{-2}$	−
$c_{VV}$	Annihilation reaction toehold complementarity degree	$0.25$	−

provides the details of all the parameters.

From Figure 11, we note that there was an agreement between the output $y=y^{+}-y^{-}$ reference tracking behaviour, as well as the concentration evolution in the DSD reactions and the CRN concentrations obtained by the MAK, as in (26) to (29). Despite the DSD reactions circuit exhibiting a desirable tracking behaviour, the dynamics of the transient were damped further. One probable reason for this is the additional auxiliary species and bimolecular reactions that were present.

4. State Feedback Control System

The closed-loop dynamics of the static state feedback can be modified by the controller utilising only the gains on the plant state and adding no dynamics to the open-loop system. The plant considered here was the classic double-integrator that represents the simplest second-order system. Thus, there was an extra state for the feedback apart from the output feedback.

Compared to the previous example, the control of this plant was more difficult as it was a marginally stable system due to the presence of two poles at the origin. Hence, the closed-loop system not only had to achieve the reference tracking capability, but to stabilise the open-loop system as well.

Example 4

(Double-integrator with linear state feedback). The state space description of the process to be controlled is

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right]& =\left[\begin{array}{cc}0& 0\\ q& 0\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}q\\ 0\end{array}\right]r\end{array}$$

(49)

where each integration has a gain q (Figure 12), and reference tracking was achieved with the control law:

$$\begin{array}{c}\hfill v=r-k_{1}x-k_{2}y\end{array}$$

(50)

There are two parameters, $k_1$ and $k_2$, that can be used to adjust the closed-loop state space system dynamics:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right]& =\left[\begin{array}{cc}-q{k}_1& -q{k}_2\\ q& 0\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}q\\ 0\end{array}\right]r\end{array}$$

(51)

The closed-loop frequency response results in a second-order system with the transfer function:

$$\begin{array}{c}\hfill Y\left(s\right)=\frac{q^2}{s^2+q{k}_{1}s+q^2{k}_2}R\left(s\right)\end{array}$$

(52)

where the poles describing the transient response are given by

$$\begin{array}{c}\hfill \lambda =\frac{q}{2}\left(-k_1\pm \sqrt{k_1^2-4k_2}\right)\end{array}$$

(53)

There are only three parameters in the closed-loop system, i.e., two controller gains and one plant gain. From (52), there are q that define the system timescale, $1/k_2$, which is the static gain, ${\omega }_n=q\sqrt{k_2}$ that denote frequency, and $\xi =\frac{k_1}{2\sqrt{k_2}}$ that represent the damping coefficient. The parameter $k_2$ is required to be set to unity for achieving steady-state reference tracking. This implies that any implementation-related error or deviation in this parameter will be visible in the steady-state error. It follows also that, for an overdamped response, $\xi >1\Rightarrow k_1>2$.

4.1. Representation with Chemical Reactions

Similar to the previous example, further simplification of the CRN can be achieved through the combination of the integration of the first state with the sum of the feedback contributions and reference. The additional reactions to represent the sum operation suggested in the following studies [12,13,14] can be avoided through this simplification.

By considering the dual-rail representation, the eight catalysis and two annihilation reactions resulting from the CRN are given by

$$R^{+}\stackrel{q}{\to }{R}^{+}+X^{+},R^{-}\stackrel{q}{\to }{R}^{-}+X^{-}$$

(54)

$$X^{+}\stackrel{q}{\to }{X}^{+}+Y^{+},X^{-}\stackrel{q}{\to }{X}^{-}+Y^{-}$$

(55)

$$X^{+}\stackrel{q{k}_1}{\to }{X}^{+}+X^{-},X^{-}\stackrel{q{k}_1}{\to }{X}^{-}+X^{+}$$

(56)

$$Y^{+}\stackrel{q{k}_2}{\to }{Y}^{+}+X^{-},Y^{-}\stackrel{q{k}_2}{\to }{Y}^{-}+X^{+}$$

(57)

$$X^{+}+X^{-}\stackrel{\eta }{\to }Ø,Y^{+}+Y^{-}\stackrel{\eta }{\to }Ø$$

(58)

The plant double-integrator is represented by (54) and (55), which depicts the chain of two catalysis reactions. The negative gains meanwhile are represented by the reactions in (56) and (57), which can be realised through the exchange of contributions between the dual-rail species, as shown in Figure 13a. Lastly, to ensure that the concentrations are maintained at feasible levels, the annihilation reactions are implemented in (58).

The MAK for the chemical network results in

$${\dot{x}}^{+}=q\left(r^{+}+k_2{y}^{-}+k_1{x}^{-}\right)-\eta x^{+}{x}^{-}$$

(59)

$${\dot{x}}^{-}=q\left(r^{-}+k_2{y}^{+}+k_1{x}^{+}\right)-\eta x^{+}{x}^{-}$$

(60)

$${\dot{y}}^{+}=q{x}^{+}-\eta y^{+}{y}^{-}$$

(61)

$${\dot{y}}^{-}=q{x}^{-}-\eta y^{+}{y}^{-}$$

(62)

From the reversed contributions in (59) and (60), the negative signs appear in the gains in the I/O dynamics of $\dot{x}={\dot{x}}^{+}-{\dot{x}}^{-}$ and $\dot{y}={\dot{y}}^{+}-{\dot{y}}^{-}$, given by

$$\dot{x}=qr-q{k}_{2}y-q{k}_{1}x$$

(63)

$$\dot{y}=qx$$

(64)

In (63) to (64), we recover the linear closed-loop dynamics.

Having the catalysis reactions $X^{\pm }\stackrel{q{k}_1}{\to }{X}^{\pm }+X^{\mp }$ crossing the production of each dual-species $X^{\mp }$ together with a fast annihilation $X^{+}+X^{-}\stackrel{\eta }{\to }Ø$ corresponding to the alternative catalytic degradation proposed in [13], with a self-repressing gain analogous to having $X^{\pm }\stackrel{q{k}_1}{\to }Ø$ (for fast $\eta \gg q{k}_1$), hence we can alternatively replace the catalysis in (56) with degradation reactions. Interestingly, the absence of degradation reactions in the produced species to ensure bounded concentrations and the marginal stability of the stoichiometry in the catalysis reactions can be related to the marginal stability of the integrators in the plant. In this construction, the stabilising feedback gain on x can be directly related to degradation reactions on $x^{\pm }$, which introduce a stable pole in the state.

This results in the CRN from Figure 13b with six catalysis, two degradation, and two annihilation reactions given by

$$R^{+}\stackrel{q}{\to }{R}^{+}+X^{+},R^{-}\stackrel{q}{\to }{R}^{-}+X^{-}$$

(65)

$$X^{+}\stackrel{q}{\to }{X}^{+}+Y^{+},X^{-}\stackrel{q}{\to }{X}^{-}+Y^{-}$$

(66)

$$X^{+}\stackrel{q{k}_1}{\to }Ø,X^{-}\stackrel{q{k}_1}{\to }Ø$$

(67)

$$Y^{+}\stackrel{q{k}_2}{\to }{Y}^{+}+X^{-},Y^{-}\stackrel{q{k}_2}{\to }{Y}^{-}+X^{+}$$

(68)

$$X^{+}+X^{-}\stackrel{\eta }{\to }Ø,Y^{+}+Y^{-}\stackrel{\eta }{\to }Ø$$

(69)

Here, we have a different MAK, i.e.,

$$\begin{array}{cc}\hfill {\dot{x}}^{+}& =q\left(r^{+}+k_2{y}^{-}-k_1{x}^{+}\right)-\eta x^{+}{x}^{-}\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill {\dot{x}}^{-}& =q\left(r^{-}+k_2{y}^{+}-k_1{x}^{-}\right)-\eta x^{+}{x}^{-}\hfill \end{array}$$

(71)

albeit that $\dot{x}$ has the same I/O dynamics as in (63). In Figure 14, the comparison of the linear design with the CRN representation is shown. There is a good agreement of the trajectories between the linear control design and the dual-signals obtained from the CRN I/O dynamics, which follows the predefined reference tracking behaviour.

We must keep in mind, however, that the dynamics of the CRN are still a nonlinear system. Writing the MAK of (65) to (69) in the new coordinates from Definition 5, we obtain

$$\dot{\overline{x}}=q\left(\overline{r}-k_2\overline{y}-k_1\overline{x}\right)$$

(72)

$$\dot{\overline{y}}=q\overline{x}$$

(73)

$$\dot{\widehat{x}}=q\left(\widehat{r}+k_2\widehat{y}+k_1\widehat{x}\right)-\frac{\eta }{2}{\widehat{x}}^2+\frac{\eta }{2}{\overline{x}}^2$$

(74)

$$\dot{\widehat{y}}=q\widehat{x}-\frac{\eta }{2}{\widehat{y}}^2+\frac{\eta }{2}{\overline{y}}^2$$

(75)

Again, we have the static state feedback in the I/O system, with nonlinear dynamics unobservable in the representation of the linear feedback system. Although the negative state feedback gains in $\dot{\overline{x}}$ become positive state feedback gains in the dynamics of $\dot{\widehat{x}}$, the dynamics in (72) to (75) can be shown to be bounded [23].

4.2. Construction with Strand Displacement Reactions

For the DSD representation, the catalysis, degradation, and annihilation reactions are constructed again according to Section 3.2, applying the architectures from Figure 7, Figure 8 and Figure 9, to obtain DSD reactions equivalent to the CRN in (65) to (69).

Following the above discussion, the construction in (65) to (69) is adhered to for further simplification to the circuit. Here, we replace the catalysis reactions that are responsible for realising the state feedback with the degradation reactions given in (67) that utilise fewer species and a simpler template complex shown in Figure 8. Nevertheless, solely depending on catalysis (also annihilation) reactions for circuit building may have its own benefit. When proposing a catalytic degradation scheme, the authors in [13] advocated that spatial localisation of the catalysis reaction enables degradation at faster rates.

Moreover, through the Visual DSD simulation and analysis, we observed that the circuit was still functional even with the annihilation reaction $X^{+}+X^{-}\stackrel{\eta }{\to }Ø$ omitted considering their concentrations stayed low. This is yet another plausible simplification, depending on the experimental setup.

These aforementioned simplifications resulted in the implementation of state feedback having the same complexity level as the integral control problem (i.e., six catalysis, two degradation, and one annihilation reactions), despite the former needing to control a more complex second-order system with marginal stability with more degrees of freedom. In the Visual DSD simulation, there were, respectively, 15 and 20 double-stranded complexes and auxiliary single-stranded species that were initialised at $C_{max}={10}^4$ nM (the consumption was not reversible and replenished). Furthermore, the maximum toehold binding and unbinding rates were set to ${10}^{-3}{\left(\mathrm{nMs}\right)}^{-1}$ and $0.1{\mathrm{s}}^{-1}$, respectively. For the details of all the parameters and the supplied auxiliary strands, see

Table 2. The state feedback control circuit’s parametrisation for the Visual DSD simulation.

Variable	Description	Vales	Units
$C_{max}$	Template and auxiliary species’ initial concentrations	${10}^4$	nM
$k_{bnd}$	Auxiliary strands’ (<tp>,<tq>) toehold maximum binding rate	${10}^{-3}$	${\left(\mathrm{nMs}\right)}^{-1}$
$k_{ubnd}$	Unbinding rate	$0.1$	${\left(\mathrm{s}\right)}^{-1}$
$k_t$	Signal species’ (<tr>,<ty>,<tx>) toehold maximum binding rate	${10}^{-4}$	${\left(\mathrm{nMs}\right)}^{-1}$
$c_q$	Plant reaction toehold complementarity degree	$8\times {10}^{-3}$	−
$c_1$	Degradation of $X^{\pm }$ toehold complementarity degree	$8\times {10}^{-3}$	−
$c_2$	Feedback catalysis of $Y^{\pm }$ toehold complementarity degree	$8\times {10}^{-3}$	−
$c_{YY}$	Annihilation reaction toehold complementarity degree	1	−

and Figure 15, respectively.

In Figure 14, the reference-tracking behaviour of the DSD circuits is shown. We note that there was good agreement between the CRN I/O dynamics with the linear design. The DSD reactions on the other hand had slower transient dynamics. When we compared the concentrations of the CRN with the Visual DSD simulation counterpart (Figure 16), we also observed slower dynamics and lower equilibrium for $Y^{\pm }$ and $X^{\pm }$, respectively, suggesting state x was subjected to higher damping.

Nevertheless, the dynamics were similar to the desired design of the state feedback. This observation provided us with another strong indication that we can omit the annihilation reaction for $X^{\pm }$ given its limiting degradation.

5. Impacts of Nonlinearities and Uncertainty on Experimental Implementations

The goal of the proposed circuits was to represent two linear feedback systems with chemical reactions, with a number of species and reactions, which can be implemented with currently available capabilities.

Synthetic DNA and strand displacement reactions have been shown to be applicable for biocomputation and implementation of dynamical systems [9,24] and scalable to dozens of distinct molecules [26,32]. Distributed setups can increase even further the number of reactions and species, where the management of the chemical interactions is defined by spatial separation or co-localisation [30], rather than the specificity and orthogonality of the sequences of the programmed toeholds [62]. The constructions proposed in this work relied on systems of DSD reactions of a smaller scale and were based on already implemented schemes such as cooperative hybridisation [32,61].

Despite the obvious potential and capabilities of DNA-based chemistry, from the existing demonstrations of circuits based on DSD reactions, none dealt with the representation of feedback control. The biocomputation of linear feedback systems with nucleic acids remains to be validated experimentally, and beyond the basic validation of the designs, there are additional properties that deserve investigation with an experimental setup.

5.1. Positive Equilibrium

It was pointed out in Section 3.1.1 that the linear I/O systems rely on an internal positive and nonlinear representation, albeit unobservable and stable. For a linear IPR, the stability of the IPR implicates a single equilibrium at the origin [59]. In the case of the dual-rail representation applied here and, specifically in the representation of feedback with integral action, previous work from the authors showed that the positive nonlinear dynamics of the CRN realisation can admit an unforced positive equilibrium, which defines the levels of concentration at which the circuit operates and leads to persistent consumption of chemical species [23].

Take again Example 3, assuming that $k_i$ was designed for a stable closed-loop system. Then, $\overline{v}$ and $\overline{y}$ in (32) to (35) have stable trajectories, and for $r^{\pm }=0$ (no input concentrations), the equilibrium of the I/O system at ${\overline{y}}^{*}=0$ and ${\overline{v}}^{*}=0$ is globally asymptotically stable. For the internal dynamics with $r^{\pm }=0$, the unforced equilibrium conditions result in

$$\begin{array}{cc}\hfill 0& =k_i\widehat{y}-\frac{\eta }{2}{\widehat{v}}^2\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill 0& =b\widehat{v}+a\widehat{y}-\frac{\eta }{2}{\widehat{y}}^2\hfill \end{array}$$

(77)

Solving for $\widehat{y}$, we can write

$$\widehat{v}=\frac{\eta }{2b}{\widehat{y}}^2-\frac{a}{b}\widehat{y}\Rightarrow 0=k_i\widehat{y}-\frac{\eta }{2}{\left(\frac{\eta }{2b}{\widehat{y}}^2-\frac{a}{b}\widehat{y}\right)}^2\Rightarrow \left(\frac{{\eta }^3}{8b^2}{\widehat{y}}^3-\frac{a{\eta }^2}{2b^2}{\widehat{y}}^2+\frac{\eta a^2}{2b^2}\widehat{y}-k_i\right)\widehat{y}=0$$

(78)

Besides the trivial solution $\widehat{y}=0$, we have from Descartes’ rules of signs at least another equilibrium $\widehat{y}={\widehat{y}}^{*}>0$, and

$$\widehat{y}=0\Rightarrow \widehat{v}=0$$

(79)

$${\widehat{y}}^{*}>0\Rightarrow {\widehat{v}}^{*}=+\sqrt{\frac{2k_i}{\eta }{\widehat{y}}^{*}}>0$$

(80)

The reasons for the positive unforced concentrations have been theoretically predicted [23] and can be quickly understood for Examples 3 and 4. Writing the unforced dynamics ($\overline{v}=\overline{y}=r^{\pm }=0$) of the internal dynamics from (32) to (35), we have linear and nonlinear components with

$$\begin{array}{c}\hfill \left[\begin{array}{c}\dot{\widehat{v}}\\ \dot{\widehat{y}}\end{array}\right]=\left[\begin{array}{cc}0& k_i\\ b& -a\end{array}\right]\left[\begin{array}{c}\widehat{v}\\ \widehat{y}\end{array}\right]-\frac{\eta }{2}\left[\begin{array}{c}{\widehat{v}}^2\\ {\widehat{y}}^2\end{array}\right]\end{array}$$

(81)

The nonlinear quadratic contribution is always negative and stabilising. However, for the linear term, we have the eigenvalues of the matrix as the solutions of $s\left(s+a\right)-k_{i}b=0$ (which are different from the characteristic polynomial of the I/O system) and

$$\begin{array}{cc}\hfill 0& =s^2+as-k_{i}b\Rightarrow \lambda =\frac{1}{2}\left(-a\pm \sqrt{a^2+4k_{i}b}\right)\end{array}$$

(82)

Due to the positivity of the parameters, we have $k_i>0\Rightarrow \lambda =\frac{1}{2}\left(-a+\sqrt{a^2+4k_{i}b}\right)>0$. This means the matrix in the internal positive dynamics has modes that push the unforced dynamics away from the origin. The positive eigenvalue is, in fact, the Frobenius eigenvalue, and it can be used to show that the origin is an unstable unforced equilibrium of the MAK [23]. We have the same issue present in Example 4, where, from the dynamics in (72) to (75), we can write the unforced response of the positive dynamics with

$$\begin{array}{c}\hfill \left[\begin{array}{c}\dot{\widehat{x}}\\ \dot{\widehat{y}}\end{array}\right]=q\left[\begin{array}{cc}{k}_1& k_2\\ 1& 0\end{array}\right]\left[\begin{array}{c}\widehat{x}\\ \widehat{y}\end{array}\right]-\frac{\eta }{2}\left[\begin{array}{c}{\widehat{x}}^2\\ {\widehat{y}}^2\end{array}\right]\end{array}$$

(83)

and the eigenvalues of the matrix in the linear term has an unstable positive eigenvector with the associated eigenvalue ${\forall }_{k_1,k_2>0},\lambda =\frac{q}{2}\left(k_1+\sqrt{k_1^2+4k_2}\right)>0$.

The nonlinear quadratic terms ensure the system is bounded, and although for cascaded systems of CRNs without feedback, the annihilation reactions are a practical mechanism to keep the concentrations low, for the representation of feedback, including such reactions may become essential to keep the concentrations bounded [23].

5.2. Positive Concentrations and Persistent Consumption

Note that the levels of the operating concentrations and the existence of an unforced equilibrium have direct experimental consequences, such as the continuous depletion of auxiliary strands in the implementation with bimolecular strand displacement reactions.

The histories in Figure 17 show the consumptions of the auxiliary species from the simulation in Visual DSD of the state feedback control in Example 4. The consumed auxiliary species are converted into the outputs of the DSD cascades, which includes the signal species, but also the unreactive waste species that accumulate as the auxiliary species are consumed. Figure 18 shows the histories of the concentrations of the accumulated waste species that result from the DSD implementation in Example 4, where the species $Wy$ and $W{y}^{\prime }$ result from the implementation of $Y^{+}+Y^{-}\stackrel{}{\to }Ø$, while $JoinX\_1$ and $Join{X}^{\prime }\_1$ result from sequestering the signal species in $X^{\pm }\stackrel{}{\to }Ø$. At the end of the DSD cascades for the implementations of the catalysis reactions, there is a waste species: $ForkXR\_4$ and $Fork{X}^{\prime }{R}^{\prime }\_4$ result from $X^{\pm }\stackrel{}{\to }{X}^{\pm }+R^{\pm }$; $ForkYX\_4$ and $Fork{Y}^{\prime }{X}^{\prime }\_4$ result from $X^{\pm }\stackrel{}{\to }{X}^{\pm }+Y^{\pm }$; $ForkX{Y}^{\prime }\_4$ and $Fork{X}^{\prime }Y\_4$ result from $Y^{\pm }\stackrel{}{\to }{Y}^{\pm }+X^{\mp }$.

The simulation shows the consequences of a positive equilibrium, with the continuous consumption of auxiliary species and the production of waste, even when the output is reaching the steady-state, in Figure 14. The consumption of the auxiliary species should be slow to ensure the approximation of large and constant concentrations of the auxiliary species, during the complete duration of the operation of the DSD circuit. This calls for some care in the parametrisation and experimental setup to ensure the concentrations of auxiliary species remain large during the duration of operation of the circuit or the mechanisms for the replenishment of the auxiliary species and removal of waste species.

The properties of dual-rail implementations of feedback controllers predicted by the above theory have yet to be experimentally tested. The systems proposed in this paper are simple, but they still constitute powerful examples, wherein these theoretical results [22,23] can be tested. The role of annihilation reactions in the implementation of feedback systems, to bound the signal species and to slow down the consumption of auxiliary species and the production of waste species, deserves immediate experimental investigation.

5.3. Variability in the Reaction Rates

The underlying assumption of the dual-rail representation of I/O systems is that it relies on perfectly tuned reaction rates. In the derivation of the I/O dynamics in (30) to (31), it is assumed in the CRN in (21) to (25) that the pairs of equations have exactly the same reaction rates.

However, such an assumption is hindered by the design of affinities based on the complementarity of toeholds, which suffer from granularity and variability [5]. For example, the implementation of the catalysis assumes the designs of different toehold sequences with the same maximum hybridisation (see Remark 1), which we can hope to be similar, but not exactly the same. The differences will change the equilibrium conditions of the reversible reactions and introduce error and uncertainty into the effective rate of output release.

Moreover, besides the hybridisation rates, there is the approximation of large and constant concentrations of the auxiliary species from Assumption 2. Although, theoretically, the concentration of individual auxiliary species can be used to fine-tune the speed of the equivalent unimolecular reaction rates [11], experimental variations will always introduce errors and variability. In Figure 7, we have for catalysis that $k_i\propto C_{max}\times c_i\times k_t$, which is the limiting reaction in (36) to (41), and in Figure 8, for degradation, we have $a\propto C_{max}\times c_a\times k_t$. Hence, the accuracy of the initial concentrations $C_{max}$ impacts the resulting reaction speed.

Consider Example 3, where, despite the system’s simplicity, the physical parameters in the implementation of the degradation reactions impact directly on the dynamic response of the represented I/O system. If we look at the eigenvalues of the closed-loop represented in (19) given by

$$\begin{array}{c}\hfill \lambda =-\frac{a}{2}\pm \frac{1}{2}\sqrt{a^2-4k_{i}b}\end{array}$$

(84)

we have that, for a high gain $k_i>\frac{a^2}{4b}$, it results in $\Re \left\{\lambda \right\}=-\frac{a}{2}$. The real part of the pole depends directly on the implementation of the degradation reactions, and the variability in the designs of the toeholds $c_a\times k_t$, as well as the impact of the variations of $\left[JoinY\right]\left(t\right)$ during operation, can be observed in the response time.

With a realisation with strand displacement reactions, the parametrisation of the CRN will necessarily depend on several experimental effects, which introduce variability and uncertainty in the reaction rates. The known issues of the dual-rail representation with the mismatching of the reaction rates [11,12] were formally addressed by the author’s theoretical analysis and the characterisation of the asymmetric parametrisation of the dual-rail representations in [22,23].

In particular, the authors showed in [23] that, under parametric variability, the representation of a stable linear model does not necessarily result in a CRN with stable dynamics and the stability of the realisation with DSD reactions. The dual-rail representation relies on the assumption of perfectly symmetrical pairs of chemical reactions, which is broken when considering variability in the parametrisation of the reaction rates. Asymmetric parametrisation of the pairs in the CRNs causes feedback between the input-to-output dynamics and the internal nonlinear positive dynamics of the kinetics, which can lead to unstable feedback within the network, and stability must be analysed for the complete CRN of the dual-rail representation.

However, the emergence of the behaviours predicted by theoretical and numerical analysis and simulation in Visual DSD [22,23] still need validation with actual variability from experimental realisation with DSD reactions.

5.4. Robustness

One of the main questions for experimental testing are how the physical parameters impact the reaction rates of the CRN, not only for nominal performance, but also for robustness.

For example, in practice, experimental variability leads to mismatches not only on the representation of the plants with $a^{+}\ne a^{-}$, $b^{+}\ne b^{-}$, and $q^{+}\ne q^{-}$, but also in the representation of the control laws with $k_i^{+}\ne k_i^{-}$,${\left(q{k}_1\right)}^{+}\ne {\left(q{k}_1\right)}^{-}$, and ${\left(q{k}_2\right)}^{+}\ne {\left(q{k}_2\right)}^{-}$.

The reasons are clearer in the rotated dynamics. Take again Example 3, but now assuming independent reaction rates for all reactions. The MAK then results:

$${\dot{v}}^{+}=k_i^{+}{r}^{+}+k_i^{-}{y}^{-}-\eta v^{+}{v}^{-}$$

(85)

$${\dot{y}}^{+}=b^{+}{v}^{+}-a^{+}{y}^{+}-\eta y^{+}{y}^{-}$$

(86)

$${\dot{v}}^{-}=k_i^{-}{r}^{-}+k_i^{+}{y}^{+}-\eta v^{+}{v}^{-}$$

(87)

$${\dot{y}}^{-}=b^{-}{v}^{-}-a^{-}{y}^{-}-\eta y^{+}{y}^{-}$$

(88)

Expressing the dynamics in the coordinates from Definition 5, we obtain

$$\begin{array}{cc}\hfill \dot{\overline{v}}& =k_i^{+}{r}^{+}-k_i^{-}{r}^{-}+k_i^{-}{y}^{-}-k_i^{+}{y}^{+}\hfill \end{array}$$

(89)

$$\begin{array}{cc}\hfill \dot{\overline{y}}& =b^{+}{v}^{+}-b^{-}{v}^{-}-a^{+}{y}^{+}+a^{-}{y}^{-}\hfill \end{array}$$

(90)

$$\begin{array}{cc}\hfill \dot{\widehat{v}}& =k_i^{+}{r}^{+}+k_i^{-}{r}^{-}+k_i^{-}{y}^{-}+k_i^{+}{y}^{+}-2\eta v^{+}{v}^{-}\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill \dot{\widehat{y}}& =b^{+}{v}^{+}+b^{-}{v}^{-}+a^{+}{y}^{+}+a^{-}{y}^{-}-2\eta y^{+}{y}^{-}\hfill \end{array}$$

(92)

Remark 2.

Note that, for a signal $\overline{y}=y^{+}-y^{-}$ and $\widehat{y}=y^{+}+y^{-}$ and applying the same notation to a parameter a such that $\overline{a}=a^{+}-a^{-}$, $\widehat{a}=a^{+}+a^{-}$, we have the equivalences

$$a^{+}{y}^{+}-a^{-}{y}^{-}=\frac{\widehat{a}}{2}\overline{y}+\overline{a}\frac{\widehat{y}}{2}$$

(93)

$$a^{+}{y}^{+}+a^{-}{y}^{-}=\frac{\overline{a}}{2}\overline{y}+\widehat{a}\frac{\widehat{y}}{2}$$

(94)

The correspondences in (93) to (94) allow us to express the rotated dynamics as

$$\begin{array}{cc}\hfill \dot{\overline{v}}& =\frac{\widehat{k_i}}{2}\left(\overline{r}-\overline{y}\right)+\frac{\overline{k_i}}{2}\left(\widehat{r}-\widehat{y}\right)\hfill \end{array}$$

(95)

$$\begin{array}{cc}\hfill \dot{\overline{y}}& =\frac{\widehat{b}}{2}\overline{v}-\frac{\widehat{a}}{2}\overline{y}+\frac{\overline{b}}{2}\widehat{v}-\frac{\overline{a}}{2}\widehat{y}\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill \dot{\widehat{v}}& =\frac{{\widehat{k}}_i}{2}\left(\widehat{r}+\widehat{y}\right)-\frac{\eta }{2}\left({\widehat{v}}^2-{\overline{v}}^2\right)+\frac{{\overline{k}}_i}{2}\left(\overline{r}+\overline{y}\right)\hfill \end{array}$$

(97)

$$\begin{array}{cc}\hfill \dot{\widehat{y}}& =\frac{\widehat{b}}{2}\widehat{v}+\frac{\widehat{a}}{2}\widehat{y}-\frac{\eta }{2}\left({\widehat{y}}^2-{\overline{y}}^2\right)+\frac{\overline{b}}{2}\overline{v}+\frac{\overline{a}}{2}\overline{y}\hfill \end{array}$$

(98)

If the parameters are symmetric, ${\overline{k}}_i=\overline{a}=\overline{b}=0$, we recover the desired I/O system decoupled from the positive and nonlinear dynamics. With mismatching reaction rates, we have instead interconnected inputs from the nonlinear to I/O system that change fundamentally the properties of the CRN, and robustness must be addressed for the complete CRN.

Even if the nonlinearities do not participate in the I/O dynamics directly, $\dot{\overline{v}}$ and $\dot{\overline{y}}$ are no longer decoupled from internal positive dynamics. The Assumption 1 of a represented stable I/O system no longer provides guarantees of stability for the interdependent dynamics. Moreover, we can see that, for an unforced response, $\widehat{r}=0$, the equilibrium conditions:

$$\begin{array}{cc}\hfill 0& =-\frac{\widehat{k_i}}{2}\overline{y}-\frac{\overline{k_i}}{2}\widehat{y}\hfill \end{array}$$

(99)

$$\begin{array}{cc}\hfill 0& =\frac{\widehat{b}}{2}\overline{v}-\frac{\widehat{a}}{2}\overline{y}+\frac{\overline{b}}{2}\widehat{v}-\frac{\overline{a}}{2}\widehat{y}\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill 0& =\frac{{\widehat{k}}_i}{2}\widehat{y}-\frac{\eta }{2}\left({\widehat{v}}^2-{\overline{v}}^2\right)+\frac{{\overline{k}}_i}{2}\overline{y}\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill 0& =\frac{\widehat{b}}{2}\widehat{v}+\frac{\widehat{a}}{2}\widehat{y}-\frac{\eta }{2}\left({\widehat{y}}^2-{\overline{y}}^2\right)+\frac{\overline{b}}{2}\overline{v}+\frac{\overline{a}}{2}\overline{y}\hfill \end{array}$$

(102)

can have non-zero equilibrium since ${\widehat{v}}^{*}>0$ and ${\widehat{y}}^{*}>0$ lead to ${\overline{v}}^{*}\ne 0$ and ${\overline{y}}^{*}\ne 0$, with

$$\begin{array}{cc}\hfill {\overline{y}}^{*}& =-\frac{\overline{k_i}}{\widehat{k_i}}\widehat{y}\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill {\overline{v}}^{*}& =-\frac{\overline{b}}{\widehat{b}}\widehat{v}+\widehat{y}\left(-\frac{\widehat{a}}{\widehat{b}}\frac{\overline{k_i}}{\widehat{k_i}}+\frac{\overline{a}}{\widehat{b}}\right)\hfill \end{array}$$

(104)

The experimental variability of the reaction rates entail that, even in the representation of stable linear feedback systems, the error and uncertainty in the parametrisation of the reaction rates, introduced by the implementation with DSD reactions, can lead to a violation of the symmetry assumptions of the dual-rail methodology and result in unstable circuits [23]. The work by the authors in [22] investigated the impact of such variability and mismatch of the reaction rates in the stability of the mass action kinetics of the CRN, based on robust control techniques based on the structured singular value, adapted to address uncertainty, positivity, and nonlinearities. For small variations in the reaction rates and mismatches in the parametrisation of the dual-rail representation, the method can provide a representative and quantified robustness stability margin for the dual-rail-representation of a feedback system.

Therefore, we know how to quantify the robustness of the CRN to inform how much error and uncertainty introduced by the realisation with DSD reaction networks can be tolerated by the chemical representation. However, we lack experimental assessments of such theoretical results, even with the most basic example of the dual-rail representation of feedback.

The proposed systems with reduced complexity rely on a very few parameters, even considering the complete MAK of the CRN. That means that fewer uncertainty intervals need to be characterised and simpler models for robustness analysis can be generated, while still providing examples where theoretical robustness results could be tested experimentally.

5.5. Leakage Dynamics

Finally, one of the most-problematic issues in the implementation of DSD networks is the existence of leakage, where the output strands are released even in the absence of input. Leakage without exposed toeholds does exist experimentally, albeit at much lower rates, and strands that are designed to be unreactive will unbind and expose domains and toeholds. For example, leakage in the previously unidirectional in the (42) reaction can be represented with

$$\begin{array}{ccc}\hfill Y+\mathit{JoinY}& \underset{k_{leak}}{\overset{c_a\times k_t}{{\displaystyle \rightleftharpoons }}}& JoinY\_1+dY\hfill \end{array}$$

(105)

With a finite unbinding rate $k_{leak}$, irreversible reactions in the constructed CRN are actually implemented with reversible cascades of DSD reactions, and such reversibility can change the effective reaction rate and cause the buffering of the input species [8,13].

Moreover, even if the leakage rate can be assumed to be very low $k_{leak}\ll k_t\times c_i$, the large concentrations of auxiliary strands used in the DSD cascades amplifies the leakage of the output, and its impact on the circuit operation can become relevant [24]. Consider in Example 3 that there are additional and undesired processes that lead to the release of $Y^{\pm }$, where

$$Z^{+}\stackrel{{\rho }^{+}}{\to }{Z}^{+}+Y^{+}$$

(106)

$$Z^{-}\stackrel{{\rho }^{-}}{\to }{Z}^{-}+Y^{-}$$

(107)

In this case, these reactions model the presence of strands $Z^{\pm }$ that interact with auxiliary species at high concentrations, which can be approximated with unimolecular reactions, where ${\rho }^{\pm }\approx C_{max}\times k_{leak}$.

Assumption 3.

Assume that the same leakage processes exist for $Y^{+}$ and $Y^{-}$ and have the same network structure. As the network of unimolecular catalysis and degradation reactions are duplicated, we consider that the leakage processes are also duplicated.

From the construction methodology of the dual-rail representation, the cascades of DSD reactions for catalysis and degradation have a symmetrical structure. That is, the representation always results in pairs of catalysis and degradation reactions with ideally matching reaction rates. See, for example, Figure 5 and Figure 13, where the symmetry of the construction methodology is clear in the graphs of the networks. With the additional leakage reactions from (106) to (107), we then have in the MAK that

$$\begin{array}{cc}\hfill {\dot{v}}^{+}& =k_i^{+}{r}^{+}+k_i^{-}{y}^{-}-\eta v^{+}{v}^{-}\hfill \end{array}$$

(108)

$$\begin{array}{cc}\hfill {\dot{y}}^{+}& =b^{+}{v}^{+}-a^{+}{y}^{+}-\eta y^{+}{y}^{-}+{\rho }^{+}{z}^{+}\hfill \end{array}$$

(109)

$$\begin{array}{cc}\hfill {\dot{v}}^{-}& =k_i^{-}{r}^{-}+k_i^{+}{y}^{+}-\eta v^{+}{v}^{-}\hfill \end{array}$$

(110)

$$\begin{array}{cc}\hfill {\dot{y}}^{-}& =b^{-}{v}^{-}-a^{-}{y}^{-}-\eta y^{+}{y}^{-}+{\rho }^{-}{z}^{-}\hfill \end{array}$$

(111)

and

$$\begin{array}{cc}\hfill \dot{\overline{v}}& =\frac{\widehat{k_i}}{2}\left(\overline{r}-\overline{y}\right)+\frac{\overline{k_i}}{2}\left(\widehat{r}-\widehat{y}\right)\hfill \end{array}$$

(112)

$$\begin{array}{cc}\hfill \dot{\overline{y}}& =\frac{\widehat{b}}{2}\overline{v}-\frac{\widehat{a}}{2}\overline{y}+\frac{\overline{b}}{2}\widehat{v}-\frac{\overline{a}}{2}\widehat{y}+\frac{\widehat{\mathit{\rho }}}{\mathbf{2}}\overline{\mathit{z}}+\frac{\overline{\mathit{\rho }}}{\mathbf{2}}\widehat{\mathit{z}}\hfill \end{array}$$

(113)

$$\begin{array}{cc}\hfill \dot{\widehat{v}}& =\frac{{\widehat{k}}_i}{2}\left(\widehat{r}+\widehat{y}\right)-\frac{\eta }{2}\left({\widehat{v}}^2-{\overline{v}}^2\right)+\frac{{\overline{k}}_i}{2}\left(\overline{r}+\overline{y}\right)\hfill \end{array}$$

(114)

$$\begin{array}{cc}\hfill \dot{\widehat{y}}& =\frac{\widehat{b}}{2}\widehat{v}+\frac{\widehat{a}}{2}\widehat{y}-\frac{\eta }{2}\left({\widehat{y}}^2-{\overline{y}}^2\right)+\frac{\overline{b}}{2}\overline{v}+\frac{\overline{a}}{2}\overline{y}+\frac{\widehat{\mathit{\rho }}}{\mathbf{2}}\widehat{\mathit{z}}+\frac{\overline{\mathit{\rho }}}{\mathbf{2}}\overline{\mathit{z}}\hfill \end{array}$$

(115)

Even in the absence of input $\widehat{r}=0$, the leakage will trigger the production of $\widehat{y}$, and indirectly also $\widehat{v}$, since in the representation of the feedback, the leaked output species $Y^{\pm }$ will trigger the closed-loop response.

An interesting result that needs experimental investigation is to what extent the symmetry of the dual-rail representation helps with the impact of leakage. For illustration, consider Example 3 with ideal $k_i$, b, and a parameters, but with the existence of leakage with asymmetric rates ${\rho }^{+}\ne {\rho }^{-}$. Then

$$\begin{array}{cc}\hfill \dot{\overline{v}}& =k_i\left(\overline{r}-\overline{y}\right)\hfill \end{array}$$

(116)

$$\begin{array}{cc}\hfill \dot{\overline{y}}& =b\overline{v}-a\overline{y}+\frac{\widehat{\mathit{\rho }}}{\mathbf{2}}\overline{\mathit{z}}+\frac{\overline{\mathit{\rho }}}{\mathbf{2}}\widehat{\mathit{z}}\hfill \end{array}$$

(117)

$$\begin{array}{cc}\hfill \dot{\widehat{v}}& =k_i\left(\widehat{r}+\widehat{y}\right)-\frac{\eta }{2}\left({\widehat{v}}^2-{\overline{v}}^2\right)\hfill \end{array}$$

(118)

$$\begin{array}{cc}\hfill \dot{\widehat{y}}& =b\widehat{v}+a\widehat{y}-\frac{\eta }{2}\left({\widehat{y}}^2-{\overline{y}}^2\right)+\frac{\widehat{\mathit{\rho }}}{\mathbf{2}}\widehat{\mathit{z}}+\frac{\overline{\mathit{\rho }}}{\mathbf{2}}\overline{\mathit{z}}\hfill \end{array}$$

(119)

and we see that the dynamics of the I/O system depend on the asymmetry of the leakage reaction rate $\overline{\rho }$ and the difference of leaking strands $Z^{\pm }$. If the initial conditions are the same for $Z^{\pm }$ such that $z^{+}\left(t\right)\approx z^{-}\left(t\right)$, then $\widehat{\rho }\overline{z}\approx 0$, and even if there is a considerable presence of $Z^{\pm }$, if the leakage process is symmetric such that $\overline{\rho }\approx 0$, then $\overline{\rho }\widehat{z}\approx 0$.

We then have that, while the presence of $Z^{\pm }$ can result in serious leakage for each of the components of the signal $y^{\pm }$, with the dual-rail representation of the I/O system, the impact of leakage may be mitigated by reducing the differences between leaked strands and leakage reaction rates. The proposed systems can be used to investigate how the symmetry of the leaking process can be more relevant than the amplitude of the leakage, and if the dual-rail representation can result in systems that are less sensitive to leakage.

6. Conclusions

There currently exist mature theoretical frameworks for the design of linear feedback controllers with CRNs. Systematic procedures and software tools provide a translation to equivalent reactions based on nucleic acid chemistry, which should be amenable to experimental implementation. However, the readiness level of this technology has not followed the theoretical developments, and we are lacking experimental validation of such systems.

The viability of the two proposed constructions is increased by the very few DSD reactions and strands necessary for the chemical representation of linear feedback control systems, with the objective of a near-future experimentation and validation of feedback circuits based on toehold-mediated displacement reactions. For the representation of integral control, we propose a representation simpler than found in the literature, by combining in the same pairs of CRNs the operations of gain, subtraction, and integration. With this, we removed from the implementation two catalytic, two degradation, and one annihilation reaction. The representation of the state feedback candidate was simplified by combining the gain, subtraction, and integration of the first state in the same CRNs. Although a double-integrator is not a stable plant, its representation does not need degradation reactions, and with these choices, four catalytic and one annihilation reaction were removed from the representation.

Both representations capture traditional dynamical features of general classes of linear feedback control systems and can be implemented using six catalytic, two degradations, and one or two annihilation reactions. The reduced number of chemical reactions and the size of the DSD networks puts the proposed candidate constructions within current capabilities of experimental implementation.

Although the representations were simplified to fit viable experimentation with current technical capabilities, the examples still represent challenging feedback control systems, and the circuits are interesting for experimental investigation of the dependence of closed-loop dynamics on toehold design and the impact of the annihilation reactions and internal nonlinear positive dynamics on the stability, persistent equilibrium, and concentration levels of the circuit. Furthermore, the proposed systems are rich enough to investigate spurious effects and dynamics present in DSD networks and how much the differential nature of the dual-rail representation may help mitigate experimental perturbations.