text:
\section{Introduction}
The LHCb collaboration recently reported the observation of a narrow doubly charmed tetraquark candidate $T_{cc}^+$ in the prompt production of proton-proton collisions~\cite{LHCb:2021vvq,LHCb:2021auc}. This is the first observation of the exotic state containing double charm quarks. The $T_{cc}^+$ is observed in the $D^0D^0\pi^+$ mass spectrum, and has mass of approximately 3875 MeV, which is just below the $D^{*+}D^0$ threshold. Using the Breit-Wigner parametrization, the location of the resonance peak relative to the $D^{*+}D^0$ threshold $\delta_m$ and the width $\Gamma$ are determined to be
\begin{eqnarray}\label{Xstates}
		\delta_m &\equiv& m_{T_{cc}^+}-(m_{D^{*+}}+m_{D^0}) \nonumber \\
		&=& -273\pm 61\pm 5^{+11}_{-14}\ \mbox{keV}, \nonumber\\
	    \Gamma &=&	410\pm 165\pm 43^{+18}_{-38}\ \mbox{keV} ,
\end{eqnarray}
respectively.

The long-lived $T_{cc}^+$ particle has the quark content $cc\bar{u}\bar{d}$ and the spin-parity quantum number $J^P=1^+$. The flavor quantum number is absolutely exotic. In fact, the tetraquark state with double heavy quarks has been studied for many years. There are two popular theoretical pictures describing the underlying structures of such tetraquark states, i.e., the compact tetraquark picture and hadronic molecule picture. In the compact tetraquark picture, the doubly heavy state is usually thought to be composed of the compact diquark and anti-diquark~\cite{Ballot:1983iv,Zouzou:1986qh,Carlson:1987hh,Heller:1986bt,Deng:2018kly,Yan:2018gik,Richard:2018yrm,Meng:2020knc,Meng:2021yjr,Qin:2020zlg,Navarra:2007yw,Dias:2011mi,Du:2012wp,Agaev:2019qqn,Feng:2013kea,Braaten:2020nwp,Cheng:2020wxa,Eichten:2017ffp,Guo:2021yws,Weng:2021hje,Kim:2022mpa,Agaev:2021vur,Karliner:2021wju}, while in the hadronic molecule picture, it is composed of a pair of heavy mesons~\cite{Manohar:1992nd,Tornqvist:1993ng,Ericson:1993wy,Janc:2004qn,Ding:2009vj,Molina:2010tx,Li:2012ss,Xu:2017tsr,Tang:2019nwv,Ohkoda:2012hv,Liu:2019stu,Ding:2020dio}. In Refs.~\cite{Tan:2020ldi,Vijande:2009kj,Yang:2009zzp,Yang:2019itm}, the compact tetraquark and hadronic molecule pictures are taken into account simultaneously. According to the LHCb measurements, the rather closeness of the $T_{cc}^+$ mass to the $D^{*+}D^0$ threshold strongly favors the molecular explanation concerning the $T_{cc}^+$ nature~\cite{Chen:2021vhg,Feijoo:2021ppq,Wang:2021yld,Ren:2021dsi,Deng:2021gnb,Albaladejo:2021vln,Dai:2021vgf,Du:2021zzh,Meng:2021jnw,Ling:2021bir,Yan:2021wdl,Fleming:2021wmk,Xin:2021wcr,Agaev:2022ast,Zhao:2021cvg,Ke:2021rxd}. This is similar to the case of the famous $X(3872)$ state, which is widely supposed to be a weakly bound hadronic molecule composed of $D^*\bar{D}$/$D\bar{D}^*$. We refer to Refs.~\cite{ParticleDataGroup:2022pth,Chen:2022asf,Guo:2017jvc,Kalashnikova:2018vkv,Brambilla:2019esw,Guo:2019twa,Chen:2016qju,Esposito:2016noz,Olsen:2017bmm,Lebed:2016hpi,Ali:2017jda,Liu:2019zoy} for a review concerning the recent progress of the doubly heavy tetraquarks and many other exotic hadrons discovered in the last two decades.

The production mechanism of hadrons is closely related to their intrinsic structures. In theoretical study, the production of $T_{cc}^+$ in the $\gamma p$ scattering, $pp$ collisions and heavy ion collisions have been discussed recently~\cite{Huang:2021urd,Braaten:2022elw,Hu:2021gdg,Abreu:2022lfy}. In 2008, Liu and Zhao have ever proposed an experimental scheme for searching for the doubly heavy mesons in the missing mass spectrum of two $\Lambda_c$ final states in nucleon-nucleon collisions~\cite{Liu:2008ck}. However, up to now the $T_{cc}^+$ is only observed in the prompt production of $pp$ collisions. Searching for the $T_{cc}^+$ in more reactions is important for both confirming its existence and understanding its nature. In this work, we propose to search for the $T_{cc}^+$ and its analog $T_{cc\bar{s}}^+$ in the $B_c$ decays.

This paper is organized as follows. After introduction, in Section II we present the formulae of the rescattering amplitudes corresponding to the $B_c^+\to T_{cc}^+\bar{D}^0$ and $B_c^+\to T_{cc}^+\bar{D}^{*0}$ processes. The pertinent numerical results and discussions are also given in Section II. In Section III, the similar processes $B_c^+\to T_{cc\bar{s}}^+\bar{D}^0$ and $B_c^+\to T_{cc\bar{s}}^+\bar{D}^{*0}$ are discussed. A summary is given in Section IV.

\section{Production of $T_{cc}^+$ in $B_c^+$ decays}\label{sec:II}
If the $T_{cc}^+$ is a hadronic molecule, we can expect that it may be produced from the rescattering processes illustrated in Figs.~\ref{feyn:D0} and \ref{feyn:Dstar0}, i.e., the $B_c^+$ meson firstly decays into a charmonium and a charmed meson, the two particles then rescatter into the $\bar{D}^{(*)0}$ and the $T_{cc}^+$ via exchanging a charmed meson. In these rescattering diagrams, two charm quarks and one anti-charm quark are produced in the $B_c^+$ weak decay vertex, which creates a good environment for the formation of the double charm meson $T_{cc}^+$.

The $B_c$ decays also provide a good environment for the formation of the charmonium or charmonium-like states. In Refs.~\cite{Wu:2019vbk, Wang:2007sxa,Wang:2015rcz}, the production of tetraquark candidates $Z_c(3900)$, $Z_c(4020)$ and $X(3872)$ in $B_c$ decays has ever been discussed. The contributions from rescattering processes in weak decays are also widely discussed in the literature~\cite{Colangelo:2002mj,Colangelo:2003sa,Cheng:2004ru,Cheng:2005bg,Lu:2005mx,Wu:2021cyc,Liu:2020orv,Han:2021azw,Ge:2022dsp,Liu:2016onn,Liu:2017vsf}.

\begin{figure}[htb]
		\centering
		\includegraphics[width=0.85\linewidth]{BctoTccD-feyn.pdf}
		\caption{Rescattering diagrams contributing to $B_c^+\to T_{cc}^+\bar{D}^0$. The boxes and dots represent the weak and strong vertices, respectively.}
		\label{feyn:D0}
\end{figure}
\begin{figure}[htb]
		\centering
		\includegraphics[width=0.85\linewidth]{BctoTccDstar-feyn.pdf}
		\caption{Rescattering diagrams contributing to $B_c^+\to T_{cc}^+\bar{D}^{*0}$. The boxes and dots represent the weak and strong vertices, respectively.}
		\label{feyn:Dstar0}
\end{figure}

\subsection{Rescattering amplitudes}

At the quark level, the weak decay process $B_c^+\to M_{c\bar{c}} {D}^{(*)+} $ is induced by the $\bar{b}\to \bar{c} c\bar{d}$ transition, where $M_{c\bar{c}}$ represents a charmonium state. The charm quark in $B_c^+$ is treated as a spectator.
In the naive factorization approach the effective Hamiltonian governing the process reads
\begin{eqnarray}
H_{eff}&=&\frac{G_F}{\sqrt{2}}V_{cb}^* V_{cd}\big[c_1 (\bar{b}c)_{V-A}(\bar{c}d)_{V-A} \nonumber \\
&&+c_2 (\bar{c}c)_{V-A}(\bar{b}d)_{V-A} \big],
\end{eqnarray}
where $(\bar{q}_1 q_2)_{V-A}\equiv \bar{q}_1 \gamma_\mu (1-\gamma_5) q_2 $, $G_F$ is the Fermi constant, $V_{cb}$ and $V_{cd}$ are the CKM matrix elements, and $c_1$ and $c_2$ are the perturbatively calculable Wilson coefficients. Neglecting the contributions from the nonfactorizable, color-suppressed and annihilation terms, the decay amplitude of $B_c^+\to M_{c\bar{c}} {D}^{(*)+} $ can be factorized as
\begin{eqnarray}\label{eq:fac}
&&\mathcal{A}(B_c^+\to M_{c\bar{c}} {D}^{(*)+})  \nonumber \\ &&=\frac{G_F}{\sqrt{2}}V_{cb}^* V_{cd} a_1
\bra{M_{c\bar{c}}} \bar{b}\gamma_\mu (1-\gamma_5) c \ket{B_c^+} \nonumber \\
&&\times \bra{D^{(*)+}} \bar{c}\gamma^\mu (1-\gamma_5) d \ket{0},
\end{eqnarray}
where $a_1=c_1+c_2/N_c$, and $N_c$ is the number of colors. The factorized amplitude can be expressed in terms of the form factors of the transition $B_c^+\to M_{c\bar{c}}$ and the decay constant of $D^{(*)+}$.
The form factors of $B_c$ decaying into the lower charmonia have been well investigated in the literature~\cite{Wang:2008xt,Wang:2009mi,Zhang:2023ypl,Zhu:2017lqu,Wang:2018duy,Harrison:2020gvo,Nobes:2000pm,Ivanov:2005fd,Sun:2008ew,Dhir:2008hh,Kiselev:1999sc,Ivanov:2000aj,Zuo:2006re,Kiselev:2000pp,Du:1988ws,Ebert:2003cn,Hernandez:2006gt,Rui:2016opu,Ke:2013yka,ATLAS:2022aiy,LHCb:2013kwl,CDF:2012ksy,Colangelo:1999zn,Ivanov:2006ni,Dubnicka:2017job,Chang:1992pt}. To estimate the rescattering amplitudes, in this work we employ the relevant numerical results of Refs.~\cite{Wang:2008xt,Wang:2009mi}, where the form factors are calculated by means of the covariant light-front approach. The form factors of $B_c\to$ $J/\psi$ and $\eta_c$ induced by the vector and axial-vector currents are defined by
\begin{eqnarray}
&&\bra{\eta_c(k_1)} V_\mu \ket{B_c(p_0)}= \left( P_\mu -\frac{m_{B_c}^2-m_{\eta_c}^2}{q^2}q_\mu \right) \nonumber \\
&&\times F_1^{B_c \eta_c}(q^2)+ \frac{m_{B_c}^2-m_{\eta_c}^2}{q^2}q_\mu F_0^{B_c \eta_c}(q^2), \\
&&\bra{J/\psi(k_1,\varepsilon_{J}^*)} V_\mu \ket{B_c(p_0)}=-\frac{1}{m_{B_c}+m_{J/\psi}} \nonumber \\
&&\times \epsilon_{\mu\nu\alpha\beta}\varepsilon_{J}^{*\nu}P^\alpha q^\beta V^{B_c J/\psi}(q^2), \\
&&\bra{J/\psi(k_1,\varepsilon_{J}^*)} A_\mu \ket{B_c(p_0)}=i\bigg\{ (m_{B_c}+m_{J/\psi}) \varepsilon_{J}^{*\mu} \nonumber \\ &&\times A_1^{B_c J/\psi}(q^2) -\frac{P\cdot \varepsilon_J^{*}}{m_{B_c}+m_{J/\psi}} P_\mu A_2^{B_c J/\psi}(q^2) \nonumber \\
&&-2m_{J/\psi} \frac{P\cdot \varepsilon_J^{*}}{q^2} q_\mu [ A_3^{B_c J/\psi}(q^2)-A_0^{B_c J/\psi}(q^2)]\bigg\},
\end{eqnarray}
where $P=p_0+k_1$, $q=p_0-k_1$, and $\varepsilon_J$ is the polarization vector of $J/\psi$. The $B_c$ decaying to the $P$-wave charmonia $\chi_{c0}$, $\chi_{c1}$ and $h_c$ form factors are defined by
\begin{eqnarray}
&&\bra{\chi_{c0}(k_1)} A_\mu \ket{B_c(p_0)}= \left( P_\mu -\frac{m_{B_c}^2-m_{\chi_{c0}}^2}{q^2}q_\mu \right) \nonumber \\
&&\times F_1^{B_c \chi_{c0}}(q^2)+ \frac{m_{B_c}^2-m_{\chi_{c0}}^2}{q^2}q_\mu F_0^{B_c \chi_{c0}}(q^2), \\
&&\bra{A(k_1,\varepsilon_{A}^*)} A_\mu \ket{B_c(p_0)}=-\frac{1}{m_{B_c}-m_{A}} \nonumber \\
&&\times \epsilon_{\mu\nu\alpha\beta}\varepsilon_{A}^{*\nu}P^\alpha q^\beta A^{B_c A}(q^2), \\
&&\bra{A(k_1,\varepsilon_{A}^*)} V_\mu \ket{B_c(p_0)}=-i\bigg\{ (m_{B_c}-m_{A}) \varepsilon_{A}^{*\mu} \nonumber \\ &&\times V_1^{B_c A}(q^2) -\frac{P\cdot \varepsilon_J^{*}}{m_{B_c}-m_{A}} P_\mu V_2^{B_c A}(q^2) \nonumber \\
&&-2m_{A} \frac{P\cdot \varepsilon_A^{*}}{q^2} q_\mu [ V_3^{B_c A}(q^2)-V_0^{B_c A}(q^2)]\bigg\},
\end{eqnarray}
where the state $A$ represents the axial vector meson $\chi_{c1}$ or $h_c$.

The second matrix element in the right-hand side of Eq.~(\ref{eq:fac}) is parameterized as
\begin{eqnarray}
&&\bra{D^{+}(q)} \bar{c}\gamma^\mu (1-\gamma_5) d \ket{0}=i f_D q^{\mu}, \\
&&\bra{D^{*+}(q, \varepsilon_{D^*}^*)} \bar{c}\gamma^\mu (1-\gamma_5) d \ket{0}=f_{D^*}m_{D^*}\varepsilon^{*\mu}_{D^*},
\end{eqnarray}
with $f_D$ and $f_{D^*}$ the decay constants of $D$ and $D^*$, respectively.

One of the strong vertices in the rescattering diagrams of Figs.~\ref{feyn:D0} and \ref{feyn:Dstar0} involves a charmonium and a pair of open charm mesons. Taking into account the heavy quark spin symmetry, the effective Lagrangian describing these couplings is given by \cite{Colangelo:2003sa,Guo:2010ak,Casalbuoni:1996pg}
\begin{eqnarray}
&&\mathcal{L}_{\psi}=ig_\psi \braket{
\mathcal{J} \bar{H}_{2a}  \overleftrightarrow{\slashed{\partial}} \bar{H}_{1a}
} +\mbox{H.c.}, \label{lag:psi}\\
&& \mathcal{L}_{\chi}=ig_\chi \braket{
{\chi}^\mu \bar{H}_{2a}  \gamma_\mu \bar{H}_{1a}
} +\mbox{H.c.}, \label{lag:chi}
\end{eqnarray}
with
\begin{eqnarray}
&& \mathcal{J} = \frac{1+{\rlap{v}/}}{2} [\psi^\mu \gamma_\mu-\eta_c\gamma_5]
\frac{1-{\rlap{v}/}}{2} , \\
&& \chi^\mu = \frac{1+{\rlap{v}/}}{2} \bigg[h_{c}^\mu \gamma_5+ \frac{1}{\sqrt{3}}(\gamma^\mu-v^\mu)\chi_{c0}  \nonumber \\
&& +\frac{1}{\sqrt{2}}\epsilon^{\mu\alpha\beta\gamma} v_\alpha \gamma_\beta \chi_{c1\gamma} +\chi_{c2}^{\mu\alpha}\gamma_\alpha
\bigg]\frac{1-{\rlap{v}/}}{2} ,\\
&& H_{1a}  = \frac{1+{\rlap{v}/}}{2}[\mathcal{D}_{a\mu}^*\gamma^\mu-\mathcal{D}_a\gamma_5] , \\
&& H_{2a}  = [\bar{\mathcal{D}}_{a\mu}^*\gamma^\mu+\bar{\mathcal{D}}_a\gamma_5]\frac{1-{\rlap{v}/}}{2} , \\
&&\bar{H}_{1a,2a}=\gamma^0 H_{1a,2a}^{\dag} \gamma^0, \\ &&\mathcal{D}^{(*)}=(D^{(*)0},D^{(*)+},D^{(*)+}_s) ,
\end{eqnarray}
where $A\overleftrightarrow{\partial}_\mu B\equiv A({\partial}_\mu B) -(\partial_\mu A)B$, $H_{2a}$ is the charge conjugate field of $H_{1a}$, and $a$ is the light flavor index. The $S$- and $P$-wave charmonia are collected into the fields $\mathcal{J}$ and $\chi^\mu$, respectively. The pseudoscalar and vector charmed mesons are collected into the fields $H_{1}$ and $H_2$. All of the heavy fields in the above equations contain a factor $\sqrt{m}$ with $m$ the corresponding heavy meson mass. The $\braket{\cdots}$ in Eqs.~(\ref{lag:psi}) and (\ref{lag:chi}) means the trace over Dirac matrices. According to the effective Lagrangian $\mathcal{L}_\psi$ and $\mathcal{L}_\chi$, we can obtain the corresponding vertex function $\mathcal{A}(M_{c\bar{c}}\to \bar{D}^{(*)0} {D}^{(*)0})$ in the rescattering diagram. The relevant vertex functions are collected in Appendix \ref{sec:Appendix}. The coupling constants $g_\psi$ and $g_\chi$ can be estimated by invoking the vector meson dominance arguments~\cite{Colangelo:2002mj,Colangelo:2003sa}. The results are
\begin{eqnarray}
g_\psi &=& \frac{\sqrt{m_\psi}}{2m_D f_{\psi}}, \\
g_\chi &=& \sqrt{\frac{m_{\chi_{c0}}}{3}} \frac{1}{f_{\chi_{c0}}},
\end{eqnarray}
with $f_{\psi}$ and $f_{\chi_{c0}}$ the decay constants of $J/\psi$ and $\chi_{c0}$, respectively.

If we assume that the $T_{cc}^+$ is a pure isoscalar hadronic molecule, its wave function can be written as
\begin{eqnarray}
\ket{T_{cc}^+} = \frac{1}{\sqrt{2}}\left(\ket{D^{*+}D^0} -\ket{D^{*0}D^+}\right) .
\end{eqnarray}
The vertex functions in the rescattering diagrams for the $D^{*+}D^0$ or $D^{*0}D^+$ fusing into the $T_{cc}^+$ read
\begin{eqnarray}
\mathcal{A}(D^{*+}D^0\to T_{cc}^+)&=& +\frac{g_1}{\sqrt{2}} \varepsilon_{T_{cc}^+}^{*} \cdot \varepsilon_{D^{*+}}, \label{eq:g1}\\
\mathcal{A}(D^{*0}D^+\to T_{cc}^+)&=& -\frac{g_2}{\sqrt{2}} \varepsilon_{T_{cc}^+}^{*} \cdot \varepsilon_{D^{*0}}.\label{eq:g2}
\end{eqnarray}
Taking into account the isospin symmetry, the coupling constant $g_1$ should be equal to $g_2$. However, the thresholds of $D^{*+}D^0$ and $D^{*0}D^+$ are located above $m_{T_{cc}^+}$ around 0.3 and 1.7 MeV, respectively, and in the molecular scheme the coupling of the loosely bound state $T_{cc}^+$ with its component is relatively sensitive to the binding energy. Therefore we adopt the similar prescription used in Ref.~\cite{Meng:2021jnw} to account for the difference between $g_1$ and $g_2$. The coupling constants take the form
\begin{eqnarray}
g_{1,2}=\frac{4m_{T_{cc}^+}\sqrt{\pi \kappa_{1,2}}}{\sqrt{\mu_{1,2}}},
\end{eqnarray}
where $\kappa_1=\sqrt{2\mu_1(m_{D^{*+}}+m_{D^0}-m_{T_{cc}^+})}$, $\kappa_2=\sqrt{2\mu_2(m_{D^{*0}}+m_{D^+}-m_{T_{cc}^+})}$, and $\mu_1$ ($\mu_2$) is the reduced mass of $D^{*+}D^0$ ($D^{*0}D^+$).

The decay amplitude of $B_c^+\to T_{cc}^+\bar{D}^{(*)0}$ via one of the rescattering diagrams in Figs.~\ref{feyn:D0} and \ref{feyn:Dstar0} can be expressed in a general form as follows
\begin{eqnarray}\label{eq:loopamp}
&&\mathcal{A}(B_c^+\to T_{cc}^+\bar{D}^{(*)+})=-i\int \frac{d^4 k_3}{(2\pi)^4} \frac{\mathcal{A}(B_c^+\to M_{c\bar{c}} {D}^{(*)+})}{(k_1^2-m_1^2)} \nonumber \\
&&\times \frac{\mathcal{A}(M_{c\bar{c}}\to \bar{D}^{(*)0} {D}^{(*)0})\mathcal{A}(D^{*}D\to T_{cc}^+)}{(k_2^2-m_2^2)(k_3^2-m_3^2)} \mathcal{F}(k_3^2,m_3^2),
\end{eqnarray}
where $k_{1}$, $k_{2}$ and $k_{3}$ ($m_{1}$, $m_{2}$ and $m_{3}$) correspond to the momenta (mass) of intermediate mesons $M_{c\bar{c}}$, $D^{(*)+}$ and ${D}^{(*)0}$, respectively.
The sum over polarization of intermediate state in Eq.~(\ref{eq:loopamp}) is implicit. For the intermediate spin-1 state, the sum over polarization reads
$\sum_{\mbox{pol}}\varepsilon_\mu \varepsilon_\nu^*=-g_{\mu\nu}+k_\mu k_\nu/m^2$. A dipole form factor $\mathcal{F}(k_3^2,m_3^2)$ is introduced in Eq.~(\ref{eq:loopamp}), which takes the form
\begin{eqnarray}\label{eq:formfactor}
\mathcal{F}(k_3^2,m_3^2)=\left(\frac{\Lambda^2-m_3^2}{\Lambda^2-k_3^2}\right)^2.
\end{eqnarray}
This form factor is supposed to parameterize the off-shell effects of the intermediate state and to kill the ultraviolet divergence in the loop integrals. We employ a dipole form factor rather than a monopole form factor because the latter cannot kill the divergence appearing in some rescattering amplitudes. The loop integral is performed by employing the program package \textit{LoopTools}~\cite{Hahn:1998yk}.

\subsection{Numerical analysis}
In this subsection, we give the results from explicit calculations of the rescattering amplitudes. First we list the input parameters used in the numerical calculation. In the weak decay amplitude $\mathcal{A}(B_c^+\to M_{c\bar{c}} {D}^{(*)+}) $, the central values of the CKM matrix elements $|V_{cb}|=0.0408$ and $|V_{cd}|=0.221$ reported by the Particle Data Group (PDG)~\cite{ParticleDataGroup:2022pth} are adopted, and the combination of Wilson coefficients $a_1$ is set to be 1.14~\cite{Ivanov:2006ni}. For the relevant decay constants, we use the following values: $f_D=0.209$ GeV~\cite{ParticleDataGroup:2022pth}, $f_{D^*}=0.245$ GeV~\cite{Becirevic:1998ua}, $f_\psi=0.416$ GeV~\cite{ParticleDataGroup:2022pth} and $f_{\chi_{c0}}=0.510$ GeV~\cite{Colangelo:2002mj}. Concerning the relevant particle mass, the PDG 2022 central values are used~\cite{ParticleDataGroup:2022pth}.

The numerical results of $B_c\to$ $\eta_c$, $J/\psi$, $\chi_{c0,c1}$ and $h_c$ form factors in Table II of Ref.~\cite{Wang:2008xt} and Table I of Ref.~\cite{Wang:2009mi} are employed. In our calculation of the loop integral, as an approximation we do not take into account the $q^2$-dependence of the form factors. The values of these weak decay form factors are fixed at $q^2=m_{D^{(*)+}}^2$, i.e., an on-shell approximation of the weak vertex function is adopted.

There is still a free parameter, i.e., the cutoff energy $\Lambda$, in the dipole form factor Eq.~(\ref{eq:formfactor}). Its explicit value should be determined from the experimental data. As an theoretical prediction, the empirical value of $\Lambda$ is usually set to be larger than the mass of the exchanged particle, and it also depends on the formalism of the form factor. For instance, in Ref.~\cite{Colangelo:2002mj}, one obtains $\Lambda\approx 2.7$ GeV to roughly fit the experimental data of $\mbox{Br}(B^-\to K^-\chi_{c0})$, where the form factor is the monopole type and the exchanged particles are $D$ and $D^*$. This cut-off energy is just around the typical values of the mass of the radially excited states of $D^{(*)}$ mesons. Taking into account the uncertainty induced by the cutoff $\Lambda$, we plot the branching ratio of $B_c^+\to T_{cc}^+ \bar{D}^{(*)0} $ via the rescattering processes as a function of $\Lambda$. The numerical results are displayed in Figs.~\ref{fig:D0Tcc} and \ref{fig:Dstar0Tcc}, where we consider a relatively larger range $2.5$--$5$ GeV for $\Lambda$. The uncertainties from the $T_{cc}^+$ mass are also taken into account.

The thresholds of final states $T_{cc}^+ \bar{D}^{0}$ and $T_{cc}^+ \bar{D}^{*0}$ are around $5739.7$ and $5881.7$ MeV, respectively, which are not far from the $B_c^+$ mass $6274.5$ MeV. Therefore the higher partial-wave amplitudes of the two channels are expected to be highly suppressed by the limited phase space.

From Figs.~\ref{fig:D0Tcc} and \ref{fig:Dstar0Tcc}, we can see that the two branching ratios increase monotonically with $\Lambda$ increasing. Within the cutoff range $2.5$--$5$ GeV, the branching ratio of $B_c^+\to T_{cc}^+ \bar{D}^0 $ is of the order of $10^{-7}$, while that of $B_c^+\to T_{cc}^+ \bar{D}^{*0} $ increases from $\mathcal{O}(10^{-6})$ to $\mathcal{O}(10^{-5})$ with $\Lambda$ increasing. The branching ratio of $B_c^+\to T_{cc}^+ \bar{D}^{*0} $ is much larger than that of $B_c^+\to T_{cc}^+ \bar{D}^0 $. This is because that the $S$-wave decay is allowed in $B_c^+\to T_{cc}^+ \bar{D}^{*0} $ but forbidden in $B_c^+\to T_{cc}^+ \bar{D}^0 $ to conserve the angular momentum. For the moderate cutoff $\Lambda$ around 3 GeV, $\mbox{Br}(B_c^+\to T_{cc}^+ \bar{D}^{*0})$ is of the order of $10^{-5}$, and we can expect it may be detectable in future experiments.

\begin{figure}[tb]
		\centering
		\includegraphics[width=0.8\linewidth]{TccD0bar.pdf}
		\caption{$\Lambda$-dependence of the branching ratio of $B_c^+\to T_{cc}^+ \bar{D}^0 $ via the rescattering processes in Fig.~\ref{feyn:D0}. The band is obtained by taking into account the uncertainties of the $\delta_m$.}
		\label{fig:D0Tcc}
\end{figure}

\begin{figure}[tb]
		\centering
		\includegraphics[width=0.8\linewidth]{TccDstar0bar.pdf}
		\caption{$\Lambda$-dependence of the branching ratio of $B_c^+\to T_{cc}^+ \bar{D}^{*0} $ via the rescattering processes in Fig.~\ref{feyn:Dstar0}. The band is obtained by taking into account the uncertainties of the $\delta_m$.}
		\label{fig:Dstar0Tcc}
\end{figure}

\section{Production of $T_{cc\bar{s}}^+$ in $B_c^+$ decays}
    \begin{figure}[htb]
		\centering
		\includegraphics[width=0.85\linewidth]{BctoTccsD-feyn.pdf}
		\caption{Rescattering diagrams contributing to $B_c^+\to T_{cc\bar{s}}^+\bar{D}^0$.}
		\label{feyn:TccsD}
	\end{figure}

	\begin{figure}[htb]
		\centering
		\includegraphics[width=0.85\linewidth]{BctoTccsDstar-feyn.pdf}
		\caption{Rescattering diagrams contributing to $B_c^+\to T_{cc\bar{s}}^+\bar{D}^{*0}$.}
		\label{feyn:TccsDstar}
	\end{figure}

Inspired by the observation of $T_{cc}^+$, one may guess some other analogs can also exist, such as the doubly heavy states with the strange quark. In this work, we are interested in a $T_{cc}^+$ analog named $T_{cc\bar{s}}^+$. We assume $T_{cc\bar{s}}^+$ is a hadronic molecule composed of $D^{*+}_s D^0$/$D^+_s D^{*0}$, and the wave function is defined as
\begin{eqnarray}
\ket{T_{cc\bar{s}}^+} = \frac{1}{\sqrt{2}}\left(\ket{D^{*+}_s D^0} -\ket{ D^{*0} D^+_s}\right).
\end{eqnarray}
The $T_{cc\bar{s}}^+$ mass relative to the lower $D^{*0}D_s^+$ threshold $\delta_m$ is defined as
\begin{eqnarray}\label{eq:TccWF}
		\delta_m &\equiv& m_{T_{cc\bar{s}}^+}-(m_{D^{*0}}+m_{D^{+}_s}).
\end{eqnarray}
The relation between $T_{cc}^+$ and $T_{cc\bar{s}}^+$ is similar to that between $Z_c(3900)$ and $Z_{cs}(3985)$, which are widely supposed to be the hadronic molecules composed of $D\bar{D}^*/D^*\bar{D}$ and $D_s\bar{D}^*/D_s^*\bar{D}$.
In Ref.~\cite{Ding:2020dio}, the doubly heavy systems composed of a pair of heavy mesons have been systematically studied in a quasi-potential Bethe-Salpeter equation approach. The authors predict that the bound state with $J^P=1^+$ can be found from the $D_s^*D$-$D_s D^*$ coupled channel interactions, i.e., the hadronic molecule $T_{cc\bar{s}}^+$ we defined here may exist.

We are mainly interested in the production of the $T_{cc\bar{s}}^+$ state in this work. In the $B_c^+$ decays, the production mechanism of $T_{cc\bar{s}}^+$ is rather similar with that of $T_{cc}^+$. Replacing the $D^{(*)+}$ in Figs.~\ref{feyn:D0} and \ref{feyn:Dstar0} with $D^{(*)+}_s$, the rescattering diagrams contributing to the $B_c^+\to T_{cc\bar{s}}^+\bar{D}^0$ and $B_c^+\to T_{cc\bar{s}}^+\bar{D}^{*0}$ are illustrated in Figs.~\ref{feyn:TccsD} and \ref{feyn:TccsDstar}, respectively.

In the rescattering diagrams, the effective Hamiltonian governing the weak process $B_c^+\to M_{c\bar{c}} {D}^{(*)+}_s $ reads
\begin{eqnarray}
H_{eff}&=&\frac{G_F}{\sqrt{2}}V_{cb}^* V_{cs}\big[c_1 (\bar{b}c)_{V-A}(\bar{c}s)_{V-A} \nonumber \\
&&+c_2 (\bar{c}c)_{V-A}(\bar{b}s)_{V-A} \big],
\end{eqnarray}
with $|V_{cs}|=0.975$~\cite{ParticleDataGroup:2022pth}. The decay constants in the factorized amplitudes are taken as $f_{D_s}=0.247$ GeV~\cite{ParticleDataGroup:2022pth} and $f_{D^*_s}=0.272$ GeV~\cite{Becirevic:1998ua}. Notice that $B_c^+\to M_{c\bar{c}} {D}^{(*)+}_s $ is a Cabibbo-favored process compared with $B_c^+\to M_{c\bar{c}} {D}^{(*)+} $. Therefore the branching ratio of $B_c^+\to T_{cc\bar{s}}^+\bar{D}^{(*)0}$ via the rescattering processes is also expected to be larger.

Following the rather similar procedures as that in Section \ref{sec:II}, we calculate the rescattering contributions and the numerical results are presented in Figs.~\ref{fig:D0Tccs} and \ref{fig:Dstar0Tccs}. We choose two typical values $\delta_m=100$ keV and $\delta_m=1$ MeV for a weakly bound hadronic molecule to estimate the couplings between $T_{cc\bar{s}}^+$ and its components. For the $B_c^+\to T_{cc\bar{s}}^+\bar{D}^{0}$ process, varying the $\Lambda$ from 2.5 to 5 GeV, the branching ratio increases from $\mathcal{O}(10^{-6})$ to $\mathcal{O}(10^{-5})$ for both $\delta_m=$100 keV and 1 MeV.
For the moderate cutoff $\Lambda$ around 3 GeV, $\mbox{Br}(B_c^+\to T_{cc\bar{s}}^+\bar{D}^{0})$ is of the order of $10^{-6}$.
Compared with $B_c^+\to T_{cc\bar{s}}^+\bar{D}^{0}$, the $B_c^+\to T_{cc\bar{s}}^+\bar{D}^{*0}$ process has a larger branching ratio. The argument is similar as that in Section \ref{sec:II}. For both of the two $\delta_m$ values,   $\mbox{Br}(B_c^+\to T_{cc\bar{s}}^+\bar{D}^{*0})$ increases from $\mathcal{O}(10^{-4})$ to $\mathcal{O}(10^{-3})$ with $\Lambda$ increasing.  For the moderate cutoff $\Lambda$ around 3 GeV, $\mbox{Br}(B_c^+\to T_{cc\bar{s}}^+\bar{D}^{*0})$ is of the order of $10^{-4}$. This is a sizable branching ratio. If the $T_{cc\bar{s}}^+$ state truly exists, it is very likely that we can find it in the $B_c$ decays.

The two lines in Fig.~\ref{fig:D0Tccs} has a crossing point. This is because that the $T_{cc\bar{s}}^+ D^{*0}D_s^+$ coupling is sensitive to the $\delta_m$ when $\delta_m$ is smaller. Employing the similar coupling definitions as those in Eqs.~(\ref{eq:g1}) and (\ref{eq:g2}), The $T_{cc\bar{s}}^+ D^{*+}_s D^0$ coupling $g_1$ and  $T_{cc\bar{s}}^+ D^{*0}D_s^+$ coupling $g_2$ are determined to be $(g_1,g_2)\simeq(7.1\ \mbox{GeV},3.4\ \mbox{GeV})$ when  $\delta_m=100$ keV, and $(g_1,g_2)\simeq(7.8\ \mbox{GeV},6.0\ \mbox{GeV})$ when  $\delta_m=1$ MeV. As a result the interferences among rescattering diagrams for the two cases behave differently. The $\delta_m=100$ keV line is not simply scaled by the $\delta_m=1$ MeV line.

\begin{figure}[htb]
		\centering
		\includegraphics[width=0.8\linewidth]{D0Tccs.pdf}
		\caption{$\Lambda$-dependence of the branching ratio of $B_c^+\to T_{cc\bar{s}}^+ \bar{D}^0 $ via the rescattering processes in Fig.~\ref{feyn:TccsD}. }
		\label{fig:D0Tccs}
\end{figure}
\begin{figure}[htb]
		\centering
		\includegraphics[width=0.8\linewidth]{Dstar0Tccs.pdf}
		\caption{$\Lambda$-dependence of the branching ratio of $B_c^+\to T_{cc\bar{s}}^+ \bar{D}^{*0} $ via the rescattering processes in Fig.~\ref{feyn:TccsDstar}.}
		\label{fig:Dstar0Tccs}
\end{figure}

\section{Summary}
In this work, we study the production of the doubly charmed state $T_{cc}^+$ and its analog $T_{cc\bar{s}}^+$ in $B_c$ decays, which provide a good environment for the formation of the exotic meson containing double charm quarks.  The $T_{cc}^+$ or $T_{cc\bar{s}}^+$ is produced from a charmonium and a charmed meson reascattering via exchanging another charmed meson. The contributions from various rescatterings with different intermediate states are taken into account. The calculation of rescattering amplitudes is performed under the ansatz that $T_{cc}^+$ and $T_{cc\bar{s}}^+$ are weakly bound hadronic molecules. For the moderate cutoff energy, the branching ratios $\mbox{Br}(B_c^+\to T_{cc}^+ \bar{D}^{0})$ and $\mbox{Br}(B_c^+\to T_{cc}^+ \bar{D}^{*0})$ are estimated to be of the order of $10^{-7}$ and $10^{-5}$, respectively.
The $T_{cc\bar{s}}^+$ production in the $B_c^+$ decay is a Cabibbo-favored process. For the moderate cutoff energy $\mbox{Br}(B_c^+\to T_{cc\bar{s}}^+ \bar{D}^{0})$ and $\mbox{Br}(B_c^+\to T_{cc\bar{s}}^+ \bar{D}^{*0})$ are estimated to be of the order of $10^{-6}$ and $10^{-4}$, respectively.

We should also mention that even if $T_{cc}^+$, $T_{cc\bar{s}}^+$ and some other analogs are not hadronic molecules, the production mechanism of the doubly charm mesons proposed in this paper still works. Correspondingly, the couplings between the tetraquark states and open charm mesons need to be modified.

The predicted relatively sizable branching ratios suggest that in future experiments one may search for the $T_{cc}^+$ and its analogs in the $B_c$ decay processes discussed here. Besides, investigating the different production mechanisms of these exotic doubly charmed states is also crucial in revealing their underlying structures.

\begin{acknowledgments}
This work is supported by the National Natural Science Foundation of China (NSFC) under Grants No.~11975165, No.~12235018, No.~12175165 and No.~12075167.
\end{acknowledgments}

\begin{appendix}
\section{The vertex functions $\mathcal{A}(M_{c\bar{c}}\to \bar{D}^{(*)0} {D}^{(*)0})$ in the rescattering amplitudes}\label{sec:Appendix}
Here we give the vertex functions $\mathcal{A}(M_{c\bar{c}} (k_1)\to \bar{D}^{(*)0} (p_1) {D}^{(*)0} (k_3))$ corresponding to a charmonium coupling with a pair of open charm mesons:
\begin{eqnarray}
&&\mathcal{A}(J/\psi\to \bar{D}^{0}  {D}^{0}) \nonumber \\
&&= 2 g_\psi (k_3-p_1)\cdot \varepsilon_J m_D\sqrt{m_\psi}, \\
&&\mathcal{A}(J/\psi\to \bar{D}^{0}  {D}^{*0})
\nonumber \\
&&= -4i g_\psi \epsilon_{\mu\nu\alpha\beta} k_1^\mu p_1^\nu \varepsilon_{D^*}^{*\alpha} \varepsilon_J^\beta  \sqrt{\frac{m_D m_{D^*}}{m_\psi}}, \\
&&\mathcal{A}(J/\psi\to \bar{D}^{*0}  {D}^{0})
\nonumber \\
&&= -4i g_\psi \epsilon_{\mu\nu\alpha\beta} k_1^\mu p_1^\nu \varepsilon_{\bar{D}^*}^{*\alpha} \varepsilon_J^\beta  \sqrt{\frac{m_D m_{D^*}}{m_\psi}},\\
&&\mathcal{A}(J/\psi\to \bar{D}^{*0}  {D}^{*0}) \nonumber \\
&&= 2 g_\psi m_{D^*}\sqrt{m_\psi}[-(k_3-p_1)\cdot \varepsilon_J \varepsilon_{\bar{D}^*}^* \cdot \varepsilon_{D^*}^*  \\ \nonumber
&&+(k_3-p_1)\cdot \varepsilon_{D^*}^* \varepsilon_{\bar{D}^*}^* \cdot \varepsilon_J+(k_3-p_1)\cdot \varepsilon_{\bar{D}^*}^* \varepsilon_{D^*}^* \cdot \varepsilon_J] , \\
\nonumber
&&\mathcal{A}(\eta_c \to \bar{D}^{0}  {D}^{*0}) \nonumber \\
&&= 2 g_\psi (k_3-p_1)\cdot \varepsilon_{D^*}^* \sqrt{m_{\eta_c} m_D m_{D^*}}, \\
\nonumber
&&\mathcal{A}(\eta_c \to \bar{D}^{*0}  {D}^{0}) \nonumber \\
&&= -2 g_\psi (k_3-p_1)\cdot \varepsilon_{\bar{D}^*}^* \sqrt{m_{\eta_c} m_D m_{D^*}}, \\
\nonumber
&&\mathcal{A}(\eta_c \to \bar{D}^{*0}  {D}^{*0}) \nonumber \\
&&= -4i g_\psi \epsilon_{\mu\nu\alpha\beta} k_1^\mu p_1^\nu \varepsilon_{\bar{D}^*}^{*\alpha} \varepsilon_{D^*}^{*\beta} \frac{m_{D^*}}{\sqrt{m_{\eta_c}}} , \\
\nonumber
&&\mathcal{A}(h_c \to \bar{D}^{*0}  {D}^{*0}) \nonumber \\
&&= 2i g_\chi \epsilon_{\mu\nu\alpha\beta} k_1^\mu \varepsilon_{h_c}^\nu \varepsilon_{\bar{D}^*}^{*\alpha} \varepsilon_{D^*}^{*\beta} \frac{m_{D^*}}{\sqrt{m_{h_c}}}, \\
\nonumber
&&\mathcal{A}(h_c \to \bar{D}^{*0} D^0) \nonumber \\
&&= -2 g_\chi \varepsilon_{h_c} \cdot \varepsilon_{\bar{D}^*}^* \sqrt{m_{h_c} m_D m_{D^*}}, \\
\nonumber
&&\mathcal{A}(h_c \to \bar{D}^0 D^{*0}) \nonumber \\
&&= 2 g_\chi \varepsilon_{h_c} \cdot \varepsilon_{D^*}^* \sqrt{m_{h_c} m_D m_{D^*}}, \\
\nonumber
&&\mathcal{A}(\chi_{c0} \to \bar{D}^0 D^0) \nonumber \\
&&= -2 \sqrt{3} g_\chi m_D\sqrt{m_{\chi_{c0}}}, \\
\nonumber
&&\mathcal{A}(\chi_{c0} \to \bar{D}^{*0} D^{*0}) \nonumber \\
&&= \frac{2}{\sqrt{3}}  g_\chi \varepsilon_{\bar{D}^*}^* \cdot \varepsilon_{D^*}^* m_{D^*}\sqrt{m_{\chi_{c0}}}, \\
\nonumber
&&\mathcal{A}(\chi_{c1} \to \bar{D}^{*0} D^0) \nonumber \\
&&= 2\sqrt{2}i g_\chi \varepsilon_{\chi_{c1}} \cdot \varepsilon_{\bar{D}^*}^* \sqrt{m_{\chi_{c1}}m_Dm_{{D}^*}}, \\
\nonumber
&&\mathcal{A}(\chi_{c1} \to \bar{D}^0 D^{*0}) \nonumber \\
&&= 2\sqrt{2}i g_\chi \varepsilon_{\chi_{c1}} \cdot \varepsilon_{D^*}^* \sqrt{m_{\chi_{c1}}m_Dm_{D^*}}.
\end{eqnarray}

\end{appendix}


text:
\section{Introduction}
\label{sec:introduction}
Recent advances in synthetic biology have seen the incorporation of many control engineering design principles into the construction of biomolecular circuits~\cite{Qian2018,Blanchini,Hancock2019,Siami2019}.
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~\cite{Soloveichik2010,Chen2013}.  The capability of such circuits to operate \emph{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~\cite{Hemphill2013}, multi-input logic based on DNA circuitry interacting with native mRNA~\cite{Groves2015}, and reliable strand displacement probes triggered by mRNA being transcribed into cells~\cite{Chatterjee2018}. 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~\cite{Soloveichik2010} 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 \cite{Phillips2009,lakin2011visual,Yordanova}.
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 \cite{Buisman2008}.
The use of the so-called \emph{dual-rail representation} with nucleic-acids~\cite{Seelig2006} 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~\cite{Chen2014a}, 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~\cite{Cacace2012}.
Unimolecular reactions of catalysis and degradation, and bimolecular reactions of annihilation, can then be used to construct CRNs to approximately represent transfer functions~\cite{Chiu2015},  linear feedback systems~\cite{Oishi2011,Yordanova, Paulino2019a}, and nonlinear controllers~\cite{Sawlekar2016}.

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~\cite{Paulino2019}, 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 \emph{in vitro} and \emph{in vivo}, and also highlight some of the associated technical challenges and limitations.

\subsection{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~\cite{Farina2000}, 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}\geq0$.
$\mathbf{M}\geq0$ means all elements $m_{ji}\geq0$, 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$.

\begin{figure}[!t]
	\centerline{\includegraphics[width=1\columnwidth]{figs2/fig_PI_feedback.pdf}}
	\caption{Linear negative feedback interconnection, where each signal $p_{j}$ results from representing a linear operator with chemical reactions.}
	\label{fig: Feedback-interconnection-linear}
\end{figure}
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~\cite{PeterToth;JanosErdi1989}, i.e.,
\begin{eqnarray}
&&a_{1}X_{1}+a_{2}X_{2}\xrightarrow{\gamma } b X_{3} \Rightarrow\dot{x}_{3}=b \gamma  x_{1}^{a_{1}}x_{2}^{a_{2}}\label{eq:CRN}
\end{eqnarray}
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 \eqref{eq:CRN} 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.~\ref{fig: Feedback-interconnection-linear}. 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~\cite{Oishi2011,Chiu2015} using unimolecular reactions of catalysis and degradation, and bimolecular reactions of annihilation, which are respectively given by
\begin{eqnarray}
&&X_i\xrightarrow{\gamma }X_i+X_{j},\quad X_{j}\xrightarrow{\gamma}\emptyset,\quad X_{i}+X_{j}\xrightarrow{\eta}\emptyset\label{crn: elementary CRN}
\end{eqnarray}
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^{+}$.

\section{Representation of linear feedback control systems with chemical reaction networks}\label{sec: dual rail CRN}

Each reaction in~\eqref{crn: elementary CRN} 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~\cite{Soloveichik2010}).
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~\eqref{crn: elementary CRN} 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.
\begin{assumption}\label{assum: rates are equal}
	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$
	\end{assumption}
\begin{defn}\label{defn: io dynamics}
	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$.
\end{defn}
For example, in Fig.~\ref{fig: Feedback-interconnection-linear} we represent the plant with the following set of chemical reactions
\begin{eqnarray}
X_{4}^{\pm} \xrightarrow{k_{1}^\pm}X_{4}^{\pm}+X_{5}^{\pm}\label{eq: CRN plant 1},\
X_{5}^{\pm}\xrightarrow{k_{2}^\pm}\emptyset,\
X_{5}^{+}+X_{5}^{-}\xrightarrow{\eta}\emptyset\label{crn: plant}
\end{eqnarray}
Computing the I/O dynamics from Definition~\ref{defn: io dynamics} under Assumption~\ref{assum: rates are equal} ($k_i^+=k_i^-=k_i$), the nonlinear terms cancel out and the I/O dynamics of~\eqref{eq: CRN plant 1} 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~\cite{Cacace2012}.

The linear feedback system in Fig.~\ref{fig: Feedback-interconnection-linear} is represented by combining~\eqref{crn: plant} with the CRNs of the linear operations of integration $\dot{p}_{3}=k_{0}p_{1}$, gain $p_{2}=k_Pp_{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.~\cite{Oishi2011,Yordanova,Paulino2019a}).

\begin{example}\label{exmp: ODEs CRN}
	 The mass action kinetics of the linear feedback structure in Fig~\ref{fig: Feedback-interconnection-linear} results in the following dynamics
	%
	\begin{subequations}
	\begin{eqnarray}
	\dot{x}_{1}^{\pm}&=&-\gamma^\pm_{3} x_{1}^{\pm}+\gamma^\mp_{2}x_{5}^{\mp}+\gamma^\pm_{1}r^{\pm}-\eta x_{1}^{+}x_{1}^{-} \label{eq: crn dyn nonlin first}\\
	\dot{x}_{2}^{\pm}&=&\gamma^\pm_{4}x_{1}^{\pm}-\gamma^\pm_{5}x_{2}^{\pm}-\eta x_{2}^{+}x_{2}^{-} \label{crn: gain}\\
	\dot{x}_{3}^{\pm}	&=&k^\pm_{0}x_{1}^{\pm}-\eta x_{3}^{+}x_{3}^{-}\\
	\dot{x}_{4}^{\pm}	&=&\gamma^\pm_{6}x_{2}^{\pm}+\gamma^\pm_{7}x_{3}^{\pm}-\gamma^\pm_{8}x_{4}^{\pm}-\eta x_{4}^{+}x_{4}^{-}\label{crn: sum}\\
	\dot{x}_{5}^{\pm} & =&k^\pm_{1}x_{4}^{\pm}-k^\pm_{2}x_{5}^{\pm}-\eta x_{5}^{+}x_{5}^{-}\label{eq: crn dyn nonlin last}
	\end{eqnarray}
	\end{subequations}
Under Assumption~\ref{assum: rates are equal}, the ODEs for $\dot{p}_j=\dot{x}^+_j-\dot{x}^-_j$ and $r=r^+-r^-$ result in the I/O system
\begin{eqnarray}
&&\dot{\mathbf{p}}=\mathbf{A}_{p}\mathbf{p}+\mathbf{B}_{p}r,\quad
y=p_5 \label{eq: I/O system state space}\\
&&
\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
\end{eqnarray}
\end{example}
\begin{remark}
	The I/O dynamics are an approximation to the original linear system in Fig.~\ref{fig: Feedback-interconnection-linear}, since the representations of subtraction, gain, and sum are exact only at steady state. The impact of the respective transient dynamics~(\ref{eq: crn dyn nonlin first}), (\ref{crn: gain}), and (\ref{crn: sum}) can be mitigated by increasing $\gamma_j$.
\end{remark}

\section{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}^{+}$.
\begin{assumption}\label{R11 stable}
	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.

\end{assumption}

\subsection{The dynamics in the natural coordinates are positive and nonlinear}
\begin{defn}\label{defn: nonlinear dynamics}
	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
	\begin{eqnarray}
	&&\dot{\mathbf{x}}=

	\left(\mathbf{A}^{\boxbslash}-\mathbf{D}\left\{ \left|\mathbf{a}\right|\right\}\right)\mathbf{x}
	+\mathbf{Br} -\eta\left(\mathbf{Px}\right)\circ\mathbf{x}\label{eq: non-linear model a + A}
	\end{eqnarray}
	where $\mathbf{A}=\mathbf{A}^{\boxbslash}-\mathbf{D}\left\{ \left|\mathbf{a}\right|\right\}$ and
		%
	%
	\begin{eqnarray}
	&&\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}
	\label{eq: state partition}
	\\
	%
	&&\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]\label{eq: nonlinear flux}
	\end{eqnarray}
\end{defn}
Compared to~\eqref{eq: I/O system state space}, the model in~\eqref{eq: non-linear model a + A} 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}\geq0$ and the degradation rates result in non-positive elements in the diagonal of $\mathbf{D}\left\{\mathbf{a}\right\}$ ($\mathbf{a}\leq0$).
$\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}\leq0$, and $\mathbf{A}^{\boxbslash}\mathbf{x}+\mathbf{Br}\geq0$.
With $\mathbf{B}\geq0$, $\mathbf{r}\geq0$, and $\mathbf{g}\left\{ \mathbf{x}\right\} =-\eta\left(\mathbf{Px}\right)$, the following Lemma~\ref{lm: positivity} shows that the nonlinear dynamics in their natural coordinates in~\eqref{eq: non-linear model a + A} are non-negative.
\begin{lemma}\label{lm: positivity}
	For a vector function $\mathbf{g}\left\{ \mathbf{x}\right\}$, if $\mathbf{M}\in\mathcal{M}$, $\mathbf{v}\geq0$,
	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.
\end{lemma}
\begin{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}\geq0$ and the trajectory remains in $\mathbb{R}^{+}_0$. \qed
\end{pf}
Rewriting~\eqref{eq: non-linear model a + A} according to the partition in (\ref{eq: state partition}) so that
\begin{eqnarray}
\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}\label{eq: decomposed nlin dyn}
\end{eqnarray}
we have matrices $\mathbf{A}$ and $\mathbf{B}$ structured into
\begin{subequations}\label{eq: matrix decomposed nonlinear model}
\begin{eqnarray}
&&\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)^{\boxbslash}+\mathbf{D}\left\{ \mathbf{a}_{1}^{\pm}\right\} ,
\;\mathbf{a}_{1}^{\pm}\leq0,
\; \mathbf{A}_{2}^{\pm}=\left(\mathbf{A}_{2}^{\pm}\right)^{\boxbslash}
\end{eqnarray}
\end{subequations}
$\mathbf{A}_{j}^{\pm}\in \mathcal{M}$,  and from Definition~\ref{defn: nonlinear dynamics} 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}^{-}$.

\begin{figure}[t]
	\begin{centering}
		\includegraphics[width=0.8\columnwidth]{figs2/fig_pq_cascade}
		\par\end{centering}
	\caption{Interconnection between the I/O dynamics and the underlying positive dynamics in the rotated coordinates. The dashed connection is absent with the nominal symmetric parameterisation from Definition~\ref{defn: symmetrical matrices}.
	}
	\label{fig: cascade}
\end{figure}

\subsection{The positive nonlinear dynamics are unobservable in the I/O dynamics of the linear representation}

\begin{defn}\label{defn: rotated coordinates}
	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
	\begin{eqnarray}
	&&\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}\label{eq: Psi map}
	\end{eqnarray}
	%
	%
\end{defn}
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~\ref{sec: dual rail CRN} and the remaining nonlinear positive dynamics $\dot{ \mathbf{q}}$, and infer their interconnections.
The rotated dynamics are then given by
	\begin{eqnarray}
	\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}\label{eq: rotated non linear dynamics}
	\end{eqnarray}
	%
\begin{remark}\label{rem: diagonal R22}
From the structures in Definition~\ref{defn: nonlinear dynamics} and~\eqref{eq: matrix decomposed nonlinear model} (recall that $\mathbf{a}_{1}^{\pm}\leq0$), we have that
\begin{eqnarray}
	&&\begin{aligned}
	&\mathbf{R_{22}}=
	\frac{\left(\mathbf{A_{1}^{+}+A_{1}^{-}+A_{2}^{+}+A_{2}^{-}}\right)^{\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}^{+}+A_{1}^{-}-A_{2}^{+}-A_{2}^{-}}\right)^{\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}^{+}-A_{1}^{-}-A_{2}^{+}+A_{2}^{-}}\right)^{\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}^{+}-A_{1}^{-}+A_{2}^{+}-A_{2}^{-}}\right)^{\boxbslash}}{2}
	-\frac{\mathbf{D}\left\{ \left|\mathbf{a_{1}^{+}}\right|-\left|\mathbf{a_{1}^{-}}\right|\right\}}{2}
	\end{aligned}\nonumber
\end{eqnarray}
	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}$.
\end{remark}
\begin{defn}\label{defn: symmetrical matrices}
	%

	Considering the condition of perfectly identical reaction rates from Assumption~\ref{assum: rates are equal}, 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}$.

\end{defn}
\begin{proposition}\label{lem: unobservable nominal dynamics}
	For the nominal symmetrical parameterisation in Definition~\ref{defn: symmetrical matrices}, the nonlinear dynamics are unobservable in the I/O system, due to the serial structure of the nominal rotated dynamics given by
	%
	\begin{subequations}\label{eq: rotated nominal}
	\begin{eqnarray}
	\dot{\mathbf{p}}&=&  \bar{\mathbf{R}}_{11}\mathbf{p}+\mathbf{W}_{p}\bar{\mathbf{B}}\mathbf{r}\label{eq: rotated nominal linear model}\\
	\dot{\mathbf{q}}&=&  \bar{\mathbf{R}}_{22}\mathbf{p}+\mathbf{W}_{q}\bar{\mathbf{B}}\mathbf{r}+\frac{\eta}{2}\mathbf{p}\circ\mathbf{p}-\frac{\eta}{2}\mathbf{q}\circ\mathbf{q}
	\label{eq: rotated nominal nonlinear model}
	\end{eqnarray}
	\end{subequations}
\end{proposition}
\begin{pf}
	Applying Definition~\ref{defn: symmetrical matrices} to the matrices in Remark~\ref{rem: diagonal R22}, 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~(\ref{eq: rotated nominal linear model}-\ref{eq: rotated nominal nonlinear model}) (illustrated in Fig.~\ref{fig: cascade}) means that $\mathbf{p}$ evolves independently of $\mathbf{q}$, making $\mathbf{q}$ unobservable in any output of the I/O dynamics
	\qed
	%
\end{pf}

\section{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.
\begin{defn}\label{defn: cascaded system}
	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.

\end{defn}
Cascaded strand displacement reactions are well suited to systematically build large computational and logic gate circuitry~\cite{Qian2011c}.
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~\ref{assum: rates are equal} and~\ref{R11 stable}, and from Remark~\ref{rem: diagonal R22}, 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.~\ref{fig: Feedback-interconnection-linear} 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.~\ref{fig_metzler_part}a.
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$ \commentarxiv{(Lemma~\ref{lem: equilibrium L nominal} 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.
\begin{remark}\label{lem: irreducible dyn}
	The interdependent evolution of all the states is reflected in the \emph{irreducibility} \cite{Farina2000} 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.
\end{remark}
\begin{proposition}\label{thm: general solution of eq.}
	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}$.
\end{proposition}
\begin{pf
	From the equilibrium condition for each coordinate $j$ we take the non-negative roots
	\begin{eqnarray}
	&&kq_{j}^{2}+|m_{jj}|q_{j}-\sum_{i\neq j}m_{ji}q_{i}=0\\
	&&\Rightarrow q_{j}=\frac{1}{2k}\left(-|m_{jj}|+\sqrt{m_{jj}^{2}+4k\sum_{i\neq j}m_{ji}q_{i}}\right)\geq 0\label{eq: roots feedback eq}
	\end{eqnarray}
	%
	\textbf{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.

	\textbf{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~\eqref{eq: roots feedback eq}. Note that even if $m_{jj} =0$, then $q_{j}^{*}>0$.

	Combining \textbf{i)} and \textbf{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$.\qed
\end{pf}
\begin{figure}[t]
	\begin{centering}
		\includegraphics[width=1\columnwidth]{figs2/fig_linear_network} \par
		\hspace{0.25\columnwidth}\textbf{a)}  \hspace{0.45\columnwidth} \textbf{b)}
	\end{centering}
	\caption{\textbf{a)} Network of catalysis (full line) and degradation (dashed) reactions for Example~\ref{exmp: ODEs CRN}. Without the negative feedback the system is a cascade of reactions (black arrows) from inputs $r^\pm$ to the outputs $x_5^\pm$ and $\mathbf{A}_{2}^\pm=0$. \textbf{b)} Introducing the negative feedback from $\mathbf{A}_{2}^{\pm}\geq0$ (red arrows), introduces positive feedback between positive systems.}
	\label{fig_metzler_part}
\end{figure}
\begin{example}\label{exmp: 2 state example}
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
	\begin{eqnarray}
	&&\dot{x}  =-d_{1}x-\boldsymbol{c_{2}y}+u, \quad \dot{y}  =-d_{2}y+c_{1}x
	%
	%
	\\&&\Rightarrow \begin{cases}
	\begin{aligned}
	U^{\pm}\xrightarrow{1}U^{\pm}+X^{\pm},\:X^{\pm}\xrightarrow{d_{1}}\emptyset,\:X^{+}+X^{-}\xrightarrow{k}\emptyset\\
	%
	X^{\pm}\xrightarrow{c_{1}}X^{\pm}+Y^{\pm},\:Y^{\pm}\xrightarrow{d_{2}}\emptyset,\:Y^{+}+Y^{-}\xrightarrow{k}\emptyset\\
	%
	Y^{\pm}\xrightarrow{c_{2}}Y^{\pm}+X^{\mp},\:U^{+}+U^{-}\xrightarrow{k}\emptyset
	\end{aligned}\label{crn: example 2}	\end{cases}\\
&&\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]\label{eq: example 2 R matrices}
\end{eqnarray}
\end{example}
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
\begin{eqnarray}
\begin{aligned}
k^{3}q_{1}^{4}+2k^{2}d_{1}q_{1}^{3}&+\left(d_{1}^{2}+c_{2}d_{2}\right)kq_{1}^2& \\
&+c_{2}\left(d_{2}d_{1}-c_{2}c_{1}\right)q_1=0&
\end{aligned}
\end{eqnarray}
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.
\begin{remark}\label{prop: negative vs positive feedback}
	%
	Note that the use of $\bar{\mathbf{A}}_2$ to represent negative feedback in the I/O dynamics~(\ref{eq: I/O system state space}) 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~\eqref{eq: rotated nominal nonlinear model} with $\bar{\mathbf{R}}_{22}=\bar{\mathbf{A}}_1 + \bar{\mathbf{A}}_2$

\end{remark}
In~\eqref{eq: example 2 R matrices},
$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. \commentarxiv{Further details in the Appendix}.%
\begin{remark}\label{rem: experimental q>0}
	%
	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}}^*\geq0$ 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.
\end{remark}

\ifarxiv \else  \fi

\section{Stability }

We begin by proving the following lemma which is applicable to the unforced dynamics of~\eqref{eq: non-linear model a + A} and~\eqref{eq: rotated nominal nonlinear model}.
\begin{lemma}\label{lm: stable if M hurwitz}
	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$.
\end{lemma}
\begin{pf}
From the stability of Metzler matrices~\cite{Farina2000}, $\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 $\qed
\end{pf}
With $\mathbf{g}\left\{ \mathbf{x}\right\} =-\mathbf{Px}$, Lemma~\ref{lm: stable if M hurwitz} ensures that if the network of catalysis and degradation reactions is stable, $\mathbf{A}\in\mathcal{H}$, the bimolecular reactions cannot destabilise~\eqref{eq: non-linear model a + A}.
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.~\ref{fig_metzler_part}b.
Since we cannot stabilise non-negative systems with positive gains $\mathbf{A}_2^\pm\geq 0$~\cite{Roszak2009}, 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~\cite{Cacace2012}, the presence of the bimolecular reactions are sometimes sufficient for stabilisation, even if  $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$.

 %

\subsection{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~\ref{defn: symmetrical matrices}, it is possible to provide guarantees for stability and boundedness.
\begin{proposition}\label{thm: GAS cascaded}
	The cascaded systems from Definition~\ref{defn: cascaded system} representing stable I/O dynamics, have GAS unforced nonlinear dynamics, for $\mathbf{x}>0$.
	%
\end{proposition}
\begin{pf}

	From Remark~\ref{rem: diagonal R22}, 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~\ref{lm: stable if M hurwitz} ensures $\mathbf{\dot{q}}=\mathbf{\bar{R}_{22}}\mathbf{q}-\frac{\eta}{2}\mathbf{q}\circ\mathbf{q}$
	is GAS at $\mathbf{q}= 0$.\qed
\end{pf}
\begin{remark}
	We can apply Proposition~\ref{thm: GAS cascaded} 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.
\end{remark}
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}$.
\begin{lemma}
	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.

	\label{lem: unstable nom q=0}
\end{lemma}
\begin{pf}
From applying the Frobenius-Perron theorem to Metzler matrices~\cite{Farina2000,Horn2012matrix},  $\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.

	\qed
\end{pf}
The IPR of a stable system using only linear positive systems is therefore not guaranteed to be stable~\cite{Cacace2012}. However, for the nonlinear positive dynamics~\eqref{eq: rotated nominal}, we can still ensure boundedness with the following result.
\begin{lemma}\label{lem: bounded nominal q}
	For $\mathbf{M}\in\mathcal{M}$ , $\mathbf{q}\left(0\right)>0$,
	and a bounded input $\mathbf{v}\geq0$, 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)$
\end{lemma}
\begin{pf}
		Lemma~\ref{lm: positivity} guarantees that the trajectories are nonnegative for
	$\mathbf{q}\left(0\right)>0$. If $\mathbf{M}\in\mathcal{H}$ is Hurwitz, Lemma~\ref{lm: stable if M hurwitz} 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
	%
	\begin{eqnarray*}
		&\dot{V}_1\left\{ \mathbf{q}\right\}  &
		  = \mathbf{1}^{T}\mathbf{M}\mathbf{q}+\mathbf{1}^{T}\mathbf{v}+\mathbf{1}^{T}\mathbf{D\left\{ \mathbf{q}\right\} }\mathbf{g}\left\{ \mathbf{q}\right\} \\
		 && = \mathbf{1}^{T}\mathbf{M}\mathbf{q}+\|\mathbf{v}\|_{1}+\mathbf{q}^{T}\mathbf{g}\left\{ \mathbf{q}\right\} \\

&&  		\leq \|\mathbf{M}\mathbf{q}\|_{1}+\|\mathbf{v}\|_{1}-k\mathbf{q}^{T}\mathbf{q}\\
		 && \leq \sqrt{N}\|\mathbf{M}\|_{2}\|\mathbf{q}\|_{2}+\|\mathbf{v}\|_{1}-k\|\mathbf{q}\|^2_{2}
	\end{eqnarray*}
	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$.
	\qed
\end{pf}
Applying Lemma~\ref{lem: bounded nominal q} with $\mathbf{g}\left\{ \mathbf{q}\right\} =-\frac{\eta}{2}\mathbf{q}$ to the unforced dynamics in~\eqref{eq: rotated nominal nonlinear model} we have $\|\mathbf{q}\|_{2}<\eta^{-1}2\sqrt{N}\|\bar{\mathbf{R}}_{22}\|_{2}$
\commentarxiv{(see illustration in Fig. \ref{fig: q0 Vdot} in Appendix)}.
In general, Lemma~\ref{lem: bounded nominal q} is not applicable to the nonlinear dynamics (\ref{eq: non-linear model a + A}), due to the matrix $\mathbf{P}$. Moreover, it relies on the assumption of a stable $\bar{\mathbf{R}}_{11}$.
\begin{proposition}\label{prop: stability for nominal}
	Consider the nominal dynamics in~(\ref{eq: rotated nominal linear model}-\ref{eq: rotated nominal nonlinear model}), with the symmetrical parameterisation from Assumption~\ref{assum: rates are equal}. Under Assumption~\ref{R11 stable}, the I/O dynamics~\eqref{eq: rotated nominal linear model} are stable, and the concentrations in the complete CRN are bounded and can be scaled down with a faster annihilation reaction rate $\eta$.
\end{proposition}
\begin{pf}
	Assumption~\ref{R11 stable} ensures the trajectories of $\mathbf{p}$ are bounded.
	%

	We can treat $\mathbf{p}$ as an additional input to the system~\eqref{eq: rotated nominal nonlinear model} and apply Lemma~\ref{lem: bounded nominal q} 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$.
	%
\qed
\end{pf}
	The same feedback responsible for a stable I/O linear dynamics can result in $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$ (see Remark~\ref{prop: negative vs positive feedback}).
	%
	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~\cite{Cacace2012}, where we would need $\bar{\mathbf{A}}\in\mathcal{H}$ for the IPR to be stable. Lemma~\ref{lem: bounded nominal q} lifts this constraint, albeit at the cost of a positive equilibrium.
	%
\begin{remark}\label{thm: q>0 for irreducible}
	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.
\end{remark}
\begin{remark}
	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~\cite{DeLeenheer2001,Roszak2009}.
	%
	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.
\end{remark}

\subsection{Local stability with asymmetrical parameterisation from experimental variability}

The construction of the I/O dynamics in (\ref{eq: I/O system state space}) assumes the symmetrical parameterisation in Definition~\ref{defn: symmetrical matrices}.
For a parametric analysis of the I/O system $\dot{\mathbf{p}}=\mathbf{A}_{p}\mathbf{p}+\mathbf{B}_{p}r$ Assumption~\ref{assum: rates are equal} still holds.
Hence, as long as the I/O linear dynamics are stable, Proposition~\ref{prop: stability for nominal} guarantees that the nonlinear dynamics are bounded.

Once we (realistically) allow that all the parameters in Example~\ref{exmp: ODEs CRN} can vary independently, we get an asymmetric parameterisation that deviates from Assumption~\ref{assum: rates are equal}.
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~\eqref{eq: I/O system state space} and~\eqref{eq: rotated nominal linear model})
\begin{eqnarray}
\dot{\mathbf{p}} &
=&\mathbf{R}_{11}\mathbf{p}+\mathbf{W}_{p}\mathbf{B}\mathbf{r}+\mathbf{R}_{12}\mathbf{q}
\end{eqnarray}
\begin{remark}\label{rem: analysis of interconnection}
	%
	With experimental variability, we lose the serial structure from (\ref{eq: rotated nominal nonlinear model}), and the I/O linear system and the underlying positive dynamics become interconnected (dashed connection in Fig. \ref{fig: cascade}).
	%
	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 (\ref{eq: non-linear model a + A}).
	%
\end{remark}
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 (\ref{eq: non-linear model a + A}) results in the followin
\begin{eqnarray}
\dot{\mathbf{s}} & =\left(\mathbf{A}+\eta\mathbf{J}\left\{ \mathbf{x}^{*}\right\} \right)\mathbf{s}+\mathbf{B}\mathbf{r}_{e}=\mathbf{A}_{s}\mathbf{s}+\mathbf{B}\mathbf{r}_{e}
\end{eqnarray}
If $\alpha\left\{ \mathbf{A}_s\right\} <0$ then the
system is locally exponentially stable around the equilibrium~\cite{Khalil2015}.
The equilibrium $\mathbf{x}^{*}=0$ is stable if and only if $\mathbf{A}\in\mathcal{H}$, in agreement with Lemma~\ref{lm: stable if M hurwitz}.
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.

\commentarxiv{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.}

\section{Analysis of an example nucleic acid feedback control system }
\begin{table}[t]
	\caption{ Nominal parameters for the example, and an asymmetrical parameterisation case which results in unstable dynamics.
	\label{tab: parameters}

	\centering{}
	\begin{tabular}{|p{55pt}|p{50pt}|p{100pt}|}
		\hline
		Parameter & Nominal & Asymmetrical case \\
		\hline
		$k_{1}^{\pm}$ & $0.001/$s & $0.00132/$s\\ \hline
		$k_{2}^{\pm}$ & $0.001/$s & $k_2^+= 0.001320 $/s, \par $k_2^-=0.000680$/s\\ \hline
		$\gamma_{i}^{\pm}$, \par $i=1,2,3$ & $0.004/$s & $\gamma_1^\pm=\gamma_2^{\pm}=0.00528/$s, \par $\gamma_3^\pm=0.00272/$s  \\ \hline
		$\gamma_{i}^{\pm}$, \par $i=6,7,8$ & $0.008/$s & $\gamma_{6}^{\pm}=\gamma_{7}^{+}=0.01056/$s, \par $\gamma_{7}^-=\gamma_{8}^{\pm}=0.00544/$s \\ \hline
		$\gamma_{i}^{\pm}$, \par $i=4,5$ & $4\times10^{-6}/$s & $\gamma_4^+=2.72\times 10^{-6}/$s,  \par  $\gamma_4^-=\gamma_5^{\pm}=5.28\times 10^{-6}/$s\\ \hline
		$k_{P}^{\pm}=\gamma_{4}^\pm /\gamma_{5}^\pm$ & $1$ & $k_{P}^+=0.5152$, $k_{P}^-=1$\\ \hline
		$k_{0}^{\pm}$ & $0.00045/$s & $0.000594/$s \\
		\hline
		$\eta$ & $5\times10^{5}/\text{M}/$s & $5\times10^{5}/\text{M}/$s  \\ 		\hline

	\end{tabular}
	\label{tab1}
\end{table}

To illustrate the application of the above results, we now analyse the feedback system given in Example~\ref{exmp: ODEs CRN}. We first consider the nominal parameterisation in Table~\ref{tab: parameters}, and analyse the dynamics in the natural coordinates $x^{\pm}_j$ in~(\ref{eq: non-linear model a + A}) and the I/O linear dynamics from~\eqref{eq: I/O system state space}.
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.~\ref{Ref: fig: Nominal x} where $\mathbf{p}$ and $\mathbf{q}$ are recovered with~\eqref{eq: Psi map}. The output $y=p_{5}$ tracks successfully the reference $r$, while $ \mathbf{q}\geq0$ reveals the underlying dynamics.
Since $\bar{\mathbf{R}}_{22}\notin\mathcal{H}$ the origin is unstable (Lemma~\ref{lem: unstable nom q=0}), 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$.

\begin{table}[h]
	\caption{Poles with maximum real part, for the I/O and linearised dynamics, for the nominal and asymmetrical parametrisations. }
	\label{tab: poles}

	\centering{}
	\begin{tabular}{p{45pt}| p{120pt} | p{36pt}}
		\hline
		Matrix $\mathbf{M}$& Poles corresponding to $\alpha\left\{\mathbf{M}\right\}$ & Stability \\
		\hline
		$ \bar{\mathbf{R}}_{11} $  &$-3.96 \times 10^{-6}$ & $\bar{\mathbf{R}}_{11}\in\mathcal{H}$\\
		%
		$ \bar{\mathbf{A}}_{s}$ 	&  $-3.96\times 10^{-6}$& $ \bar{\mathbf{A}}_{s}\in\mathcal{H}$\\ \hline
		%
		$ \mathbf{R}_{11} $  & 	 $-5.23\times 10^{-6}$& $\mathbf{R}_{11}\in\mathcal{H}$\\
		%
		$\mathbf{A}_{s} $  &
		$\mathbf{+3.16\times 10^{-5}}\pm i\mathbf{1.26\times 10^{-3}}$& $\mathbf{A}_{s}\notin\mathcal{H}$\\
		\hline

	\end{tabular}
\end{table}

Table~\ref{tab: poles} 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~\ref{tab: parameters}.
Perturbing the unforced nonlinear dynamics for this case around its equilibrium ($\mathbf{r}=0$),
results in the unstable response of Fig.~\ref{fig:unstable x xi zeta}.
\begin{figure}[]
	\centerline{\includegraphics[width=1\columnwidth]{figs2/fig_xpq_time_nominal}}
	\caption{Response of the CRN for the nominal parameterisation
		($x_{j}^{\pm}(0) > 0$) to a sequence of reference steps $r$.

	}\label{Ref: fig: Nominal x}

\end{figure}
\begin{figure}[]
	\centerline{\includegraphics[width=0.9\columnwidth]{figs2/fig_x_asym_time}}

	\caption{Trajectories in natural coordinates $x_{i}^{\pm}$ when parameterised with the asymmetrical rates from Table~\ref{tab: parameters} ($\mathbf{r}=0$ ).
	}
	\label{fig:unstable x xi zeta}
\end{figure}
\begin{figure}[]
	\centerline{\includegraphics[width=1\columnwidth]{figs2/fig_dec_p_q_time}}

	\caption{Simulation of the rotated dynamics of $\dot{ \mathbf{p}}$
		and $\dot{ \mathbf{q}}$ with decoupled matrix $\mathbf{R}$ where $\mathbf{R}_{21}=\mathbf{R}_{12}=0$.
	} \label{fig: independent xi and zeta}
\end{figure}
The poles in Table~\ref{tab: poles} 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.~\ref{fig: independent xi and zeta}, 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.

\section{Stability of the controller implementation with DSD reactions }

It remains to verify whether the stability properties of Example~\ref{exmp: ODEs CRN} 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 stran
, 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~\cite{Phillips2009,Zhang2018}.
Tuning the affinities of the toeholds, based on the base-pair affinities and the nucleotides sequences~\cite{Zhang2018}, allows the mapping of the desired reaction rates for the CRN into the DSD implementation.

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

\begin{figure}[t]
	\centering
	\includegraphics[width=1\columnwidth]{figs2/fig_DSD_xp_time_nominal.pdf}
	\caption{Simulation of the DSD reactions in VisualDSD for the symmetrical nominal system, with $\mathbf{x}\left(0\right)=0 $ nM and a sequence of steps on $\mathbf{r}$.
	}\label{fig: DSD nominal}

\end{figure}
\begin{figure}[t]
	\centerline{\includegraphics[width=0.95\columnwidth]{figs2/fig_DSD_xpq_time_nominal_aux.pdf}}
	\caption{Time history of concentrations of auxiliary species involved in the annihilation reactions.
		\label{fig: DSD nominal aux} }

\end{figure}
\begin{figure}[t]
	\centering
	\includegraphics[width=0.95\columnwidth]{figs2/fig_DSD_x_time_ucase.pdf}
	\caption{Simulation in VisualDSD, for the asymmetrical destabilising parameterisation, with $\mathbf{r}=0$

		\label{fig: DSD ucase}}

\end{figure}
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~\cite{lakin2011visual}.
The translation of the CRN system follows the construction proposed in~\cite{Soloveichik2010}, with a two-domain programmming structure~\cite{cardelli2013two}. The model was parameterised with the unstable parameterisation of~Table~\ref{tab: parameters}, 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.~\ref{fig: DSD nominal}, where $p_5$ tracks the step inputs of $r$. After $6\times 10^{6}$ s, the system converges to the positive equilibrium, and Fig.~\ref{fig: DSD nominal aux} 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.~\ref{fig: DSD ucase} shows that the parameterisation which destabilises the CRN also destabilises the DSD implementation, emphasising the practical relevance of the stability results.

\section{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.

\section*{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.

redpajama_set_name: RedPajamaArXiv