<figure><img src="image_1.png" /></figure>

HAL
open science

# On Estimating Derivatives of Input Signals in Biochemistry

Mathieu Hemery, François Fages

► To cite this version:

Mathieu Hemery, François Fages. On Estimating Derivatives of Input Signals in Biochemistry. CMSB 2023 - 21st International Conference on Computational Methods in Systems Biology, Sep 2023, Luxembourg City, Luxembourg. hal-04154923

HAL Id: hal-04154923

https://inria.hal.science/hal-04154923

Submitted on 7 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On Estimating Derivatives of Input Signals in Biochemistry

Mathieu Hemery and François Fages

Inria Saclay, Lifeware project-team, Palaiseau, France

mathieu.hemery@inria.fr Francois.Fages@inria.fr

**Abstract.** The online estimation of the derivative of an input signal is widespread in control theory and engineering. In the realm of chemical reaction networks (CRN), this raises however a number of specific issues on the different ways to achieve it. A CRN pattern for implementing a derivative block has already been proposed for the PID control of biochemical processes, and proved correct using Tikhonov's limit theorem. In this paper, we give a detailed mathematical analysis of that CRN, thus clarifying the computed quantity and quantifying the error done as a function of the reaction kinetic parameters. In a synthetic biology perspective, we show how this can be used to design error correcting terms to compute online functions involving derivatives with CRNs. In the systems biology perspective, we give the list of models in BioModels containing (in the sense of subgraph epimorphisms) the core derivative CRN, most of which being models of oscillators and control systems in the cell, and discuss in detail two such examples: one model of the circadian clock and one model of a bistable switch.

## 1 Introduction

Sensing the presence of molecular compounds in a cell compartment is a necessary task of living cells to maintain themselves in their environment, and achieve high-level functions as the result of low-level processes of basic biomolecular interactions. The formalism of chemical reaction networks (CRN) [11] is both a useful abstraction to describe such complex systems in the perspective of systems biology [18], and a possible molecular programming language in the perspective of synthetic biology [23, 8].

Sensing the concentration levels of molecular compounds has been well-studied in the domain of signal transduction networks. For instance, the ubiquitous CRN structure of MAPK signaling networks has been shown to provide a way to implement analog-digital converters in our cells, by transforming a continuous input signal, such as the concentration of an external hormone activating membrane receptors, into an almost all-or-nothing output signal according to some threshold value of the input, i.e. using a stiff sigmoid as dose-response input-output function [17].

The analysis of input/output functions fits well with the computational theory of CRNs. In particular, the Turing-completeness result shown in [8] for the
interpretation by Ordinary Differential Equations (ODE) of CRNs, possibly re-
stricted to elementary CRNs using mass-action law kinetics and at most bimolec-
ular reactions, demonstrates the generality of this approach to biomolecular pro-
gramming. Furthermore, it comes with an algorithm to automatically generate
a finite CRN for implementing any computable real function. Such a compiler is
implemented in our CRN modeling software BIOCHAM [4] in several forms, in-
cluding a theoretically more limited but practically more interesting framework
for robust online computation [14].

Sensing the derivative of an input molecular concentration is nevertheless
beyond the scope of this computational paradigm since it assumes that the
input molecular concentrations are stabilized at some fixed values which makes
no sense for computing the derivative. Furthermore, it is well-known that the
derivative of a computable real function is not necessarily computable [20]. We
must thus content ourselves with estimating the derivative of an input with some
error, instead of computing it with arbitrary precision as computability theory
requires.

In control theory and engineering, online estimations of input signal deriva-
tives are used in many places. Proportional Integral Derivative (PID) controllers
adjust a target variable to some desired value by monitoring three components:
the error, that is the difference between the current value and the target, its in-
tegral over a past time slice, and its current derivative. The derivative term can
improve the performance of the controller by avoiding overshoots and solving
some problematic cases of instability.

Following early work on the General Purpose Analog Computer (GPAC) [22],
the integral terms can be implemented with CRNs using simple catalytic syn-
thesis reactions such as $A \rightarrow A + B$ for integrating $A$ over time, indeed $B(T) = \int_{0}^{T} A(t)dt$. Difference terms can be implemented using the annihilation reaction $A_{+} + A_{-} \rightarrow \emptyset$ which is also used in [7, 21, 8] to encode negative values by the difference of two molecular concentrations, i.e. dual-rail encoding. This is at the basis of the CRN implementations of, for instance, antithetic PI controllers presented in [3].

For the CRN implementation of PID controllers, to the best of our knowledge
three different CRN templates have been proposed to estimate derivative terms.
The first one by Chevalier & al. [5] is inspired by bacteria's chemotaxis, but
relies on strong restrictions upon the parameters and the structure of the input
function making it apparently limited in scope. A second one proposed by Alexis
& al. [1] uses tools from signal theory to design a derivative circuit with offset
coding of negative values and to provide analytic expressions for its response. The
third one developed by Whitby & al. [24] is practically similar in its functioning
to the one we study here, differing only on minor implementation details, and
proven correct through Tikhonov's limit theorem. This result ensures that when
the appropriate kinetic rates tend to infinity, the output is precisely the derivative
of the input.

In this paper, we give a detailed mathematical analysis of that third derivative
CRN and quantify the error done as a function of the reaction kinetic param-
eters, by providing a first-order correction term. We illustrate the precision of this analysis on several examples, and show how this estimation of the derivative can be actively used with error-correcting terms to compute elementary math- ematical functions online. Furthermore, we compare our core derivative CRN to the CRN models in the curated part of BioModels.net model repository. For this, we use the theory of subgraph epimorphisms (SEPI) [13, 12] and its imple- mentation in BIOCHAM [4], to identify the models in BioModels which contain the derivative CRN structure. We discuss with some details the SEPIs found on two such models: biomodels 170, one of the smallest eukaryotes circadian clock model [2], and biomodels 318, a model of the bistable switch at the restriction point of the cell cycle [25].

The rest of the article is organized as follow. In Section 2, we provide some preliminaries on CRNs and their interpretation by ODEs. We present the core differentiation CRN in Section 3, in terms of both of some of its different possible biological interpretations, and of its mathematical properties. Section 4 develops the mathematical analysis to bound the error done by that core CRN, and give in Section 5 some examples to test the validity of our estimation and the possibility to introduce error-correcting terms. Section 6 is then devoted to the search of that derivative CRN pattern in BioModels repository and the analysis of those matching in two cases. Finally, we conclude on the perspectives of our approach to both CRN design at an abstract mathematical level, and comparison to natural CRNs to help understanding their functions.

# 2 Preliminaries on CRNs

## 2.1 Reactions and Equations

The CRN formalism allows us to represent the molecular interactions that occur on a finite set of molecular compounds or species, {$X_i$}$_{i}$$_{∈}$$_{1...n}$, through a finite set of formal (bio)chemical reactions, without prejudging their interpretation in the differential, stochastic, Petri Net and Boolean semantics hierarchy [10]. Each reaction is a triplet (R, P, f), also written R $\xrightarrow{f}$ P, where R and P are multisets of respectively reactant and product species in {$X_i$}, and f : $\mathbb{R}_+$ $\to$ $\mathbb{R}_+$ is a kinetic rate function of the reactant species. A CRN is thus entirely described by the two sets of n species and m reactions: {$X_i$}, {$R_s$} $\xrightarrow{f_s}$ {$P_s$}.

The differential semantics of a CRN associates positive real valued molecular concentrations, also noted $X_i$ by abuse of notation, and the following ODEs which define the time evolution of those concentrations:

$$ \frac{dX_i}{dt} = \sum_{s \in S} (P_s(X_i) - R_s(X_i))f_s(X), \quad (1) $$

where $P_s(X_i)$ (resp. $R_s(X_i)$) denotes the multiplicity (stoichiometry) of $X_i$ in the multiset of products (resp. reactants) of reaction s.

In the case of a mass action law kinetics, the rate function is a monomial, $f_s = k_s \prod_{x \in R_s} x$, composed of the product of the concentrations of the reactants
by some positive constant $k_s$. If all reactions have mass action law kinetics, we write the rate constant in place of the rate function $R \xrightarrow{k} P$, and the differential semantics of the CRN is defined by a Polynomial Ordinary Differential Equation (PODE).

From the point of view of the computational theory of CRNs, there is no loss of generality to restrict ourselves to elementary CRNs composed of at most bimolecular reactions with mass action law kinetics. Indeed, [8] shows that any computable real functions (in the sense of computable analysis, i.e. with arbitrary finite precision by a Turing machine), can be computed by such a CRN, using the dual-rail encoding of real values by the difference of molecular concentrations, $x = X_+ - X_-$. While our compiler ensures that the quantity $X_+ - X_-$ behaves properly, it is also important to degrade both of them with an annihilation reaction, $X_+ + X_- \xrightarrow{fast} \emptyset$, to avoid a spurious increase of their concentration. Those annihilation reactions are supposed to be faster than the other reactions of the CRN.

*Example 1.* The first example given in [8] showed the compilation of the cosine function of time, $y = \cos(t)$ in the following CRN:

$$ 
\begin{array}{l@{\quad}l@{\quad}l}
A_p \rightarrow A_p + y_p & A_m \rightarrow A_m + y_m & A_m(0) = 0, A_p(0) = 0 \\
y_m \rightarrow A_p + y_m & y_p \rightarrow A_m + y_p & y_m(0) = 0, y_p(0) = 1 \quad (2) \\
y_m + y_p \xrightarrow{fast} \emptyset & A_m + A_p \xrightarrow{fast} \emptyset &
\end{array}
$$

The last two reactions are necessary to avoid an exponential increase of the species concentration. The associated PODE is:

$$ 
\begin{array}{l@{\quad}l@{\quad}l}
d(A_m)/dt = y_p - fast * A_m * A_p & A_m(0) = 0 & \\
d(A_p)/dt = y_m - fast * A_m * A_p & A_p(0) = 0 & \\
d(y_m)/dt = A_m - fast * y_m * y_p & y_m(0) = 0 & (3) \\
d(y_p)/dt = A_p - fast * y_m * y_p & y_p(0) = 1 &
\end{array}
$$

## 2.2 CRN Computational Frameworks

The notions of CRN computation proposed in [8] and [14] for computing input/output functions, do not provide however a suitable framework for computing derivative functions. Both rely on a computation at the limit, meaning that the output converges to the result of the computation whenever the CRN is either properly initialized [8], or the inputs are stable for a sufficient period of time [14]. To compute a derivative, we cannot ask that the input stay fixed for any period of time as this would imply a null derivative. We want the output to follow « at run time » the derivative of the input.

Our question is thus as follows. Given an input species $X$ following a time course imposed by the environment $X(t)$, is it possible to perform an online computation such that we can approximate the derivative $\frac{dX}{dt}$ on the concentration of 2 output species using a dual-rail encoding?
The idea is to approximate the left derivative by getting back to its very mathematical definition:

$$ \frac{dX}{dt}(t) = \lim_{\epsilon \to 0^+} \frac{X(t) - X(t - \epsilon)}{\epsilon}, \quad (4) $$

but how can we measure $X(t - \epsilon)$?

# 3 Differentiation CRN

## 3.1 Biological intuition using a membrane

One biological intuition we may have to measure a value in a previous time is to use a membrane with a fast diffusive constant. Indeed, if we suppose that the input is the outside species, the inside species equilibrates to follow the concentration of the outside one (the input) but also suffers a lag due to the diffusion. Building upon this simple trick leads to the CRN presented in Fig. 1. As the derivative may be positive or negative, a dual-rail encoding is used for the derivative. This CRN is mainly equivalent to the derivative block proposed in [24] apart from the fact that we suppose (for the sake of clarity) that the input stay positive and no dual-rail encoding is used for it. In the case of a dual-rail encoded input, the two species need to have the same permeability through the membrane, otherwise the delay is not the same for the positive and negative parts.

The delay is thus introduced through a membrane under the assumption that the outside concentration is imposed by the environment. This conveniently explains why the kinetic rates are the same for the two monomials in the derivative of $X_{in}$, but this is not mandatory. Indeed two other settings can be used to construct such a CRN without relying on a membrane. We could use a phosphorylation and a dephosphorylation reactions where $X_{in}$ would be the phosphorylated species. Or we could, as in [24], rely on a catalytic production of $X_{in}$ by $X_{ext}$ and a degradation reaction of $X_{in}$. A drawback of these two other implementations is that they need to be tuned to minimize the difference between the rates of the two monomials in the derivative of $X_{in}$. Otherwise a proportional constant is introduced between $X_{ext}$ and $X_{in}$, and needs to be corrected by adjusting the production rates of $D_+$ and $D_-$.

However, the membrane implementation also has its own drawback as it requires the reaction $X_{ext} \to X_{ext} + D_+$ to occur through the membrane. We may think of a membrane protein $M$ that mediates this reaction ($X_{ext} + M \to X_{ext} + M + D_+$). Then, since its concentration is constant, it can simply be wrapped up in the kinetic constant of the reaction. Which of this three implementations should be chosen may depend on the exact details of the system to be built.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Hypergraph representation of the core differentiation CRN composed of one input, two outputs and one intermediate species. The input $X_{\text{ext}}$ is outside of the membrane and thus present in so large quantity that its concentration is not modified by the dynamics, once it crossed the membrane it is labelled as $X_{\text{in}}$. Each species $X_{\text{ext}}$ (resp. $X_{\text{in}}$) activates the synthesis of its part of the input: $D_+$ (resp. $D_-$). Finally a fast annihilation reaction eliminates both $D_+$ and $D_-$ so that only the highest of the two remains present, depending on the sign of the derivative of $X_{\text{ext}}$.</figcaption></figure>

## 3.2 Core differentiation CRN

Our core differentiation CRN schematized in Fig. 1 is more precisely composed of the following 7 reactions:

$$ 
\begin{array}{l@{\qquad}l}
X_{\text{ext}} \xrightarrow{k_{\text{diff}}} X_{\text{in}} & X_{\text{in}} \xrightarrow{k_{\text{diff}}} X_{\text{ext}} \\
X_{\text{ext}} \xrightarrow{k \cdot k_{\text{diff}}} X_{\text{ext}} + D_{+} & X_{\text{in}} \xrightarrow{k \cdot k_{\text{diff}}} X_{\text{in}} + D_{-} \\
D_{+} \xrightarrow{k} \emptyset & D_{-} \xrightarrow{k} \emptyset \\
D_{+} + D_{-} \xrightarrow{\text{fast}} \emptyset
\end{array}
\qquad (5) $$

The diffusion through the membrane is symmetrical with a constant $k_{\text{diff}}$ and both activations should have the same constant product $k \cdot k_{\text{diff}}$ while the degradation of the outputs should have a rate $k$. We make the assumption that the outside species $X_{\text{ext}}$ is present in large quantity so that its concentration is not affected by the dynamics of the CRN. Under this assumption, the differential
semantics is then the same as the one of the differentiation CRN proposed in [24]:

$$ 
\begin{align} 
\frac{dX_{\text{in}}}{dt} &= k_{\text{diff}}(X_{\text{ext}} - X_{\text{in}}) \nonumber \\
\frac{dD_{+}}{dt} &= kk_{\text{diff}}X_{\text{ext}} - kD_{+} - \text{fast}D_{+}D_{-} \tag{6} \\
\frac{dD_{-}}{dt} &= kk_{\text{diff}}X_{\text{in}} - kD_{-} - \text{fast}D_{+}D_{-} \nonumber
\end{align} 
$$

The derivative is encoded as $D = D_+ - D_-$ and hence obeys the equation (using the two last lines of the previous equation):

$$ 
\begin{align} 
\frac{dD}{dt} &= \frac{dD_+}{dt} - \frac{dD_-}{dt} \nonumber \\
&= kk_{\text{diff}}(X_{\text{ext}} - X_{\text{in}}) - k(D_+ - D_-) \tag{7} \\
\frac{dD}{dt} &= k \left( \frac{X_{\text{ext}} - X_{\text{in}}}{\frac{1}{k_{\text{diff}}}} - D \right) \nonumber
\end{align} 
$$

In the next section, we prove that $X_{\text{in}}$ is equal to $X_{\text{ext}}$ with a delay $\epsilon$, hence giving us our second time point $X(t - \epsilon)$, up to the first order in $\epsilon = \frac{1}{k_{\text{diff}}}$. The fractional part of the last equation is thus precisely an estimate of the derivative of $X_{\text{ext}}$ as defined in Eq. 4, with a finite value for $\epsilon$.

It is also worth remarking that such derivative circuits can in principle be connected to compute higher-order derivatives, with a dual-rail encoded input. It is well known that such estimations of higher-order derivatives can be very sensitive to noise and error, and are thus not reliable for precise computation but may be good enough for biological purposes. We will see a biological example of this kind in Section 6.2 on a simple model of the circadian clock.

# 4 Mathematical analysis of the quality of the estimation

Our first goal is to determine precisely the relation between $X_{\text{in}}$ and $X_{\text{ext}}$ when the later is enforced by the environment. Using the first line of Eq. 6, we obtain by symbolic integration:

$$ X_{\text{in}}(t) = k_{\text{diff}} \int_{0}^{\infty} \exp(-k_{\text{diff}}s) X_{\text{ext}}(t-s) ds, \qquad (8) $$

where we can see that $X_{\text{in}}$ is the convolution of $X_{\text{ext}}$ with a decreasing exponential. This convolution is not without reminding the notion of *evaluation* in the theory of distribution and has important properties of regularisation of the input function. In particular, whatever the input function is, this ensures that the internal representation is continuous and differentiable.

The interesting limit for us is when $k_{\text{diff}} \to \infty$, that is when $\epsilon = \frac{1}{k_{\text{diff}}} \to 0$. In this case, the exponential is neglectable except in a neighbourhood of the current
time and supposing that $X_{\text{ext}}$ is infinitely differentiable<fnref n="1" />, we obtain by Taylor expansion:

$$ 
\begin{aligned}
X_{\text{in}}(t) &= \int_{0}^{\infty} k_{\text{diff}} \exp(-k_{\text{diff}}s) \sum_{n=0}^{\infty} \frac{(-s)^n}{n!} X_{\text{ext}}^{(n)}(t) ds \\
&= \sum_{n=0}^{\infty} \frac{k_{\text{diff}}}{n!} X_{\text{ext}}^{(n)}(t) \int_{0}^{\infty} (-s)^n \exp(-k_{\text{diff}}s) ds
\end{aligned} 
\quad (9) $$

The integral may be evaluated separately using integration by parts and recursion:

$$ 
\begin{aligned}
I_n &= \int_{0}^{\infty} (-s)^n \exp(-k_{\text{diff}}s) ds = -n\epsilon I_{n-1} \\
&= (-1)^n (\epsilon)^{n+1} n!
\end{aligned} 
\quad (10) $$

We thus have:

$$ 
\begin{aligned}
X_{\text{in}}(t) &= \sum_{n=0}^{\infty} \frac{k_{\text{diff}}}{n!} X_{\text{ext}}^{(n)}(t) (-1)^n n! \epsilon^{n+1} \\
&= \sum_{n} (-\epsilon)^n X_{\text{ext}}^{(n)}(t) \\
&= X_{\text{ext}}(t) - \epsilon X'_{\text{ext}}(t) + \epsilon^2 X''_{\text{ext}}(t) + \dots \\
X_{\text{in}}(t) &= X_{\text{ext}}(t - \epsilon) + o(\epsilon^2).
\end{aligned} 
\quad (11) $$

Using Taylor expansion once again in the last equation somehow formalizes our intuition: the concentration of the internal species $X_{\text{in}}$ follows the time course of the external one with a delay equal to the inverse of the diffusive constant $k_{\text{diff}}$. This validates our formulation of the derivative.

Now, it is sufficient to remark that Eq. 7 has exactly the same form as the first line of Eq. 6 that we just study in length. Just replace $X_{\text{ext}}$ by the estimation of the left derivative, $X_{\text{ext}}$ by the output $D$ and the rate constant $k$ instead of $k_{\text{diff}}$. The delay approximation is thus also possible in this step and, introducing the delay $\tau = \frac{1}{k}$, we immediately obtain a precise expression for $D$:

$$ D(t) = \frac{X_{\text{ext}}(t - \tau) - X_{\text{ext}}(t - \epsilon - \tau)}{\epsilon} + o(\epsilon) + o(\tau^2). \quad (12) $$

We can see this as the secant approximation of the derivative of $X_{\text{ext}}$ with a step size $\epsilon$ and a delay $\tau$. Moreover we also know that the residual error on this expression are of first order in $\epsilon$ and second order in $\tau$.

It is well known in the field of numerical computation that the secant method provides a rather poor approximation, but it has the benefit to be the simplest one, and thus gives here a small size derivative circuit. In the hope of improving

<footnote n="1">¹ We also explore in Figures 2D and 3C what a non analyticity of $X_{\text{ext}}$ imply for our model.</footnote>
the precision, one could implement higher-order methods using several "membranes" to access the value of the function on several time points before performing the adapted computation. Such complexation would however also increase the delay between the input and output function.

# 5 Validation on simple examples

## 5.1 Verification of the delay-approximation

In this first subsection, we want to validate the approximation expressed by Eq. 11. For this, we focus on the diffusion part of our CRN: $X_{\text{ext}} \leftrightarrow X_{\text{in}}$. We make numerical simulation for 2 different values of $\epsilon$ and 2 different input functions: a sine wave and an absolute value signals. The second allowing us to see how well the delay approximation works in presence of non analyticity.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Behaviour of the simple diffusion model when fed with different input functions. We have $k_{\text{diff}} = 1$ in panel **A** and $k_{\text{diff}} = 10$ in the three others. The input function of the three first panels is an offset sine wave: $X_{\text{ext}}(t) = 1 + \sin(t)$. Panel **C** is a focus on part of the panel **B**. Panel **D** present as input a shifted absolute value: $X_{\text{ext}}(t) = |t - 2|$ as a study case of non-differentiable function. See the main text for the discussion.</figcaption></figure>

Fig. 2 shows the response of $X_{\text{in}}$ in that different condition. In panel **A**, the kinetic constant is very low so we expect our approximation to fail. Indeed, one can see that in addition to having an important delay, the output is strongly smoothed, this tends to average the variation of the input, bringing back $X_{\text{in}}$ to the average value of the input. In panel **B** the diffusion constant is increased by a factor 10. The delay approximation is now very good and we only expect an error of order $\epsilon^2 = 10^{-2}$ which can be checked with good accuracy on panel **C**. Panel **D** shows a case of a non-differentiable function in which an error of order $\epsilon = 0.1$ is visible shortly after the discontinuity and vanishes in a similar timescale.

## 5.2 Approximation of the derivative

Let us now check the behaviour of the derivative circuit. On Fig. 3, we can see the response of our derivative circuit for a sine wave and an absolute value input functions. In panels **A** and **B** we see that when the first and second order
derivatives of the input are smaller than the kinetic reaction rates, the delay approximation gives a very good picture of the response. On a complementary point of view, the panel C shows that in front of singularity, the system adapts after an exponential transient phase with a characteristic time $\tau = \frac{1}{k}$.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Simulation of the derivative circuit for $k_{\text{diff}} = k = 1$ (panel A) and $k_{\text{diff}} = k = 10$ (panels B and C) when fed with two different input signal: a sine wave (panels A and B) and an offset absolute value (panel C). In A the characteristic time of the diffusion is too large ($\epsilon = 1$) so that the delay approximation fail. In B $\epsilon = 10^{-1}$, the approximation is sounded and the output $D = D_+ - D_-$ correctly follows a cosine (in dark) with a delay of $\tau = 10^{-1}$. In C the singularity makes the derivative harder to compute. As for the diffusion, we see a typical decreasing exponential toward its correct value as shown in the inset that depict the absolute difference between D and its correct value (1 in this case) in logarithmic scale for the time between 2 and 3. (That is precisely the same time as the main figure where the inset is placed.)</figcaption></figure>

## 5.3 Using signal derivatives for online computations

Our main motivation for analyzing the differentiation CRN is to compute a function $f$ of some unknown input signal, $X_{\text{ext}}(t)$, online. that is, given a function $f$, compute a function $f(X_{\text{ext}}(t))$ Yet the differentiation CRN only allows us to approximate the derivative of the input signal. The idea is thus to implement the PODE:

$$ \frac{dY}{dt} = f'(X_{\text{ext}}(t)) \frac{dX_{\text{ext}}}{dt}, \quad Y(0) = f(X(0)) \quad (13) $$

and provide the result online on a set of internal species $Y(t) = Y_+ - Y_-$. This necessitates to compute the function $f'$ and estimate the derivative of the input. Using the formalism developed in [15, 16] we know that there exist an elementary CRN (i.e. quadratic PODE) computing $f'(X_{\text{ext}})$ for any elementary function $f$ and we just have shown that $\frac{dX_{\text{ext}}}{dt}$ can be approximated by the differentiation CRN. Therefore, in principle, any elementary function of input signals can be approximated online by a CRN.

As a toy example, let us consider the square function, $\frac{dY}{dt} = 2X_{\text{ext}}(D_+ - D_-)$, and as input, a sine wave offset to stay positive : $X_{\text{ext}}(t) = 1 + \sin(t)$.
The CRN generated by BIOCHAM according to these principles, to compute the square of the input online is:

$$
\begin{align}
X_{\text{ext}} \xrightarrow{k_{\text{diff}}} X_{\text{in}}, & \qquad X_{\text{in}} \xrightarrow{k_{\text{diff}}} X_{\text{ext}}, \nonumber \\
X_{\text{ext}} \xrightarrow{k_{\text{diff}}k} X_{\text{ext}} + D_{+}, & \qquad D_{+} \xrightarrow{k} \emptyset \nonumber \\
X_{\text{in}} \xrightarrow{k_{\text{diff}}k} X_{\text{in}} + D_{-}, & \qquad D_{-} \xrightarrow{k} \emptyset \tag{14} \\
X_{\text{in}} + D_{+} \xrightarrow{2} X_{\text{in}} + D_{+} + Y_{+}, & \qquad X_{\text{in}} + D_{-} \xrightarrow{2} X_{\text{in}} + D_{-} + Y_{-} \nonumber \\
D_{+} + D_{-} \xrightarrow{\text{fast}} \emptyset, & \qquad Y_{+} + Y_{-} \xrightarrow{\text{fast}} \emptyset \nonumber
\end{align}
$$

The first three lines implement the derivative circuit, the fourth line implements the derivative of Y and the last line provides the dual-rail encoding.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Square computation CRN when fed with an offset sine wave $X_{\text{ext}}(t) = 1+\sin(t)$. The parameters are: $k_{\text{diff}} = k = 10$, fast = $10^6$. Panel A shows a simulation of the naive CRN while panel B shows the corrected one where the derivative is integrated using a delayed input which eliminates the drift presented in the first panel.</figcaption></figure>

The numerical simulation of this CRN is depicted in Fig. 4A. One can see that while it effectively computes the square of the input, it also suffers from a strong drift. To verify if this drift comes from the delay between the input and the output, we can compute analytically the output of our network with our approximation of derivative with a delay (see the full computation in Appendix).

$$
\begin{align}
y(t) &= \int 2x(s)x'(s - \tau)ds \tag{15} \\
&\approx (1 + \sin(t))^2 + \tau t. \nonumber
\end{align}
$$

This is precisely the behaviour that can be seen on the time course of Fig. 4 A. After the integration of 20 time units, the offset is of order 2 which is exactly
what is predicted for a delay $\tau = \frac{1}{k} = 0.1$. Therefore, while it is always possible to get rid of such errors by increasing $k_{diff}$, the identification of the cause of the drift, gives us a potentially simpler path to eliminate it: using a representation of the input that is itself delayed: $X_{in} \leftrightarrow X_{delay}$, and use this delayed signal as the catalyst for the production of $Y_+$ and $Y_-$ in the place of $X_{in}$. This leads to the CRN given in Appendix (Eq. 20) for which numerical integration shows in Fig. 4B that we indeed have get rid of the drift, or said otherwise, the correct implementation for online computation is given by:

$$ \frac{dY}{dt} = f'(X_{\text{ext}}(t - \tau)) \frac{dX_{\text{ext}}}{dt}(t - \tau), \qquad (16) $$

where the delays has to be equal for the two pieces of the derivative.

# 6 Biological examples

## 6.1 BioModels repository

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Number of models in the curated part of BioModels per number of species, given with the number of models having a SEPI reduction to the differentiation CRN and ratio (black curve) between these two quantities.</figcaption></figure>

To explore the possibility that natural biochemical systems already implement a form or another of the core differentiation CRN, one can try to scan the CRN models of the BioModels repository [19]. This can be automated with the general graph matching notion of Subgraph EPImorphism (SEPI) introduced in [13, 12] to compare CRN models and identify model reduction relationships based on their graph structures. SEPI generalizes the classical notion of subgraph isomorphism by introducing an operation of node merging in addition to node deletion. Considering two bipartite graphs of species and reactions, there exists
a SEPI from $G_A$ to $G_B$ if there exists a sequence of mergings<fnref n="2" /> and deletions of nodes in $G_A$ such that the resulting graph is isomorphic to $G_B$.

More precisely, we used the SEPI detection algorithm of BIOCHAM to scan the curated models in Biomodels (after automatic rewriting with well-formed reactions [9]) and check the existence of a SEPI from each model graph to the differentiation CRN graph. Fig. 5 shows that our small differentiation CRN with 4 species is frequently found in large models. It is thus reasonable to restrict to models with no more than 10 species. Table 1 lists the models with no more than 10 species in the 700 first models of BioModels that contain our differentiation CRN. The predominance of models exhibiting oscillatory dynamics, and in particular circadian clock models is striking.

## 6.2 Circadian clock

Model biomodels 170 of the eukaryotes circadian clock proposed by Becker-Weimann & al. [2] is among the smallest models of the circadian clock displaying a SEPI reduction toward our differentiation CRN. Its influence graph is depicted in Fig. 6A, we also display in red the first SEPI found by BIOCHAM, and in green a second one obtained by enforcing the mapping from the PER/Cry species inside the nucleus to the input of the differentiation CRN. Interestingly, this model has the nucleus membrane separating the species mapped to $X_{\text{ext}}$ and the one mapped to $X_{\text{in}}$ in the second SEPI. The oscillatory behavior of this model is shown in panel B.

Now, thinking at the mathematical insight that this relation provides, it is quite natural for a CRN implementing an oscillator to evaluate its own derivative on the fly. Actually, when looking at the natural symmetry of the model, we are inclined to think that this CRN may actually be two interlocked CRNs of the derivative circuit, both computing the derivative of the output of the other, as if a second order derivative circuit was closed on itself. This is something we could easily check by imposing restrictions on the SEPI mapping. Enforcing the nucleus PerCRY protein to be mapped on $X_{\text{ext}}$ gives us the SEPI shown in green in Fig. 6A. To validate the preservation of the function of the derivative CRN given by this SEPI, we can verify that the quantities defined by summing the species that are mapped together are effectively linked by the desired derivative relation. As can be seen in Fig. 6B, the agreement is striking. One can even note that the delay of the chemical derivative is the one predicted by our theory.

The case of Fig. 6C is more complex as this part of the model seems to compute the opposite of the derivative. It is however worth noting that there is absolutely no degree of freedom in our choice of the species used in Fig. 6B and C that are entirely constrained by the SEPI given by BIOCHAM. Taking both SEPI together we see that $B\text{mal1}^{\text{nucleus}}$ and $B\text{mal1}^{\text{cytoplasm}}$ play symmetrical roles, being the input and derivative of the two displayed SEPI. Given that the

<footnote n="2">² A species (resp. reaction) node can only be merged with another species (resp. reaction) node and the resulting node inherits of all the incoming and outcoming edges of the two nodes.</footnote>
<table><caption>**Table 1.** List of models having a SEPI reduction to the differentiation CRN, given with model ID, number of species, number of reactions, and process modeled, among the first 700 models of the curated part of Biomodels with no more than 10 species.</caption><thead><tr><th>Model ID</th><th># Species</th><th># reactions</th><th>Topic</th></tr></thead><tbody><tr><td>0021</td><td>10</td><td>30</td><td>Circadian clock</td></tr><tr><td>0022</td><td>10</td><td>34</td><td>Circadian clock</td></tr><tr><td>0034</td><td>9</td><td>22</td><td>Circadian clock</td></tr><tr><td>0035</td><td>9</td><td>15</td><td>Circadian clock</td></tr><tr><td>0041</td><td>10</td><td>17</td><td>Creatine kinase</td></tr><tr><td>0065</td><td>8</td><td>16</td><td>Operon lactose</td></tr><tr><td>0067</td><td>7</td><td>16</td><td>Circadian clock</td></tr><tr><td>0075</td><td>10</td><td>13</td><td>Phosphoinositide turnover</td></tr><tr><td>0084</td><td>8</td><td>16</td><td>ERK Cascade</td></tr><tr><td>0107</td><td>9</td><td>23</td><td>Cell cycle</td></tr><tr><td>0108</td><td>9</td><td>18</td><td>Superoxide dismutase overexpression</td></tr><tr><td>0170</td><td>7</td><td>17</td><td>Circadian clock</td></tr><tr><td>0171</td><td>10</td><td>27</td><td>Circadian clock</td></tr><tr><td>0179</td><td>7</td><td>17</td><td>Cellular memory</td></tr><tr><td>0185</td><td>8</td><td>20</td><td>Circadian clock</td></tr><tr><td>0206</td><td>9</td><td>22</td><td>Circadian clock</td></tr><tr><td>0206</td><td>8</td><td>15</td><td>Glycolytic oscillations</td></tr><tr><td>0216</td><td>5</td><td>17</td><td>Circadian clock</td></tr><tr><td>0228</td><td>9</td><td>22</td><td>Cell cycle</td></tr><tr><td>0229</td><td>7</td><td>28</td><td>Circadian clock</td></tr><tr><td>0240</td><td>6</td><td>14</td><td>DegU transcriptional regulator</td></tr><tr><td>0257</td><td>8</td><td>19</td><td>Self-maintaining Metabolism</td></tr><tr><td>0262</td><td>9</td><td>14</td><td>AkT Signalling</td></tr><tr><td>0263</td><td>9</td><td>14</td><td>AkT Signalling</td></tr><tr><td>0269</td><td>9</td><td>22</td><td>Hormonal crosstalk in plant</td></tr><tr><td>0318</td><td>7</td><td>17</td><td>Bistable switch</td></tr><tr><td>0355</td><td>9</td><td>17</td><td>Calcium signalling</td></tr><tr><td>0359</td><td>9</td><td>15</td><td>Tissue factor pathway inhibitor</td></tr><tr><td>0360</td><td>9</td><td>15</td><td>Tissue factor pathway inhibitor</td></tr><tr><td>0495</td><td>8</td><td>18</td><td>Phospholipid synthetic pathways</td></tr><tr><td>0530</td><td>10</td><td>17</td><td>Cooperative gene regulation</td></tr><tr><td>0539</td><td>6</td><td>11</td><td>Mixed feedback loop</td></tr><tr><td>0563</td><td>10</td><td>17</td><td>Plant-microbe interaction</td></tr><tr><td>0586</td><td>10</td><td>23</td><td>Genetic oscillatory network</td></tr><tr><td>0587</td><td>10</td><td>23</td><td>Genetic oscillatory network</td></tr><tr><td>0590</td><td>9</td><td>40</td><td>Biosynthesis of pyrimidines</td></tr><tr><td>0615</td><td>4</td><td>34</td><td>Aggregation kinetics in Parkinson's Disease</td></tr><tr><td>0616</td><td>4</td><td>20</td><td>Resolution of inflammation</td></tr><tr><td>0619</td><td>10</td><td>13</td><td>Basic model of Acetaminophen</td></tr><tr><td>0622</td><td>10</td><td>19</td><td>Ubiquitination oscillatory dynamics</td></tr><tr><td>0632</td><td>8</td><td>14</td><td>Cell fate decision</td></tr><tr><td>0665</td><td>7</td><td>13</td><td>Interleukin-2 dynamics</td></tr><tr><td>0696</td><td>9</td><td>23</td><td>Incoherent type 1 feed-forward loop</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Circadian clock model proposed by Becker-Weimann & al. **A** Graph of the model, the grey block indicates the nucleus while the other species are present in the cytoplasm. The first SEPI found by BIOCHAM is displayed in red. A second SEPI is displayed in green where the Per/CRY complex in the nucleus is mapped to the input of our differentiation CRN. **B** Validation of the first (red) SEPI: we display the derivative of the input species here considered as the sum of its part: $X_{\text{ext}} = \text{Bmal1}_{\text{protein}}^{\text{nucleus}} + \text{Bmal1/Clock}_{\text{protein}}^{\text{nucleus}}$ and similarly for the output: $D = D_{+} - D_{-} = \text{Bmal1}_{\text{mRNA}}^{\text{cytoplasm}} + \text{Bmal1}_{\text{protein}}^{\text{cytoplasm}} - \text{Per/CRY}_{\text{protein}}^{\text{cytoplasm}}$. **C** Validation of the second (green) SEPI. As for the previous panel, species that are mapped together are simply summed for this validation. The qualitative matching of the two quantities in these two graphs are a good indication that these SEPI are meaningful.</figcaption></figure>
second SEPI introduces a negative sign, we may see this as:

$$ 
\begin{aligned}
B_{\text{mRNA}}^{\text{cytoplasm}} &= \frac{d}{dt} B_{\text{protein}}^{\text{nucleus}} \\
B_{\text{protein}}^{\text{nucleus}} &= -\frac{d}{dt} B_{\text{mRNA}}^{\text{cytoplasm}}
\end{aligned}
\qquad (17) $$

The solution of this well known equation are the sine and cosine functions, and this perfectly fits the oscillatory behaviour of this CRN. To confirm this hypothesis, we check for the presence of a SEPI from the clock model to the compiled cosine CRN presented in Eq. 2 which is effectively the case. On the other hand, there is no SEPI relation between the compiled cosine and the derivative circuit.

## 6.3 Bistable switch

The model *biomodels* 318 of a bistable switch in the context of the restriction point [25] displays a SEPI toward our derivative circuit. This model, presented in Fig. 7A, study the Rb-E2F pathway as an example of bistable switch where the presence of a (not modeled) growth factor activates the MyC protein, starting the pathway until it reach the E2F factor that constitute the output of the model. Yao & al. show that once E2F reaches a threshold, its activation becomes self sustained hence the notion of switch.

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Analysis of the bistable switch. A Schematic representation of the model and its SEPI, the smaller fonts for the species corresponds to the names used in the model provided online. Two other species are also present in the model (the phosphorylated form of RB and the Rb-E2F complex) but they are not presented in the figure of the article and are deleted during the SEPI, we thus choose to not display them either. B True derivative of MyC (in blue) along with the ones computed by the CRN (in red), see main text for the exact definition of *D*.</figcaption></figure>

The SEPI given by Biocham is worth of interest as it does not merge any species and only three reactions into one leaving all the other either untouched
or deleted, thus indicating that the pattern of the derivative is already well
present. Morevoer, MyC is mapped to the input and E2F to one part of the
output, reinforcing our intuition that the discovered SEPI is closed from the
natural functionning of the CRN.

To conform this, we run the simulation as provided by the models and display
the derivative of the MyC protein against a scaled difference of the D+ and D-
species: D = aRB - bE2F where a and b are positive constant adjusted so that
D goes to 0 at final time and are of the same magnitude as $\frac{dMyC}{dt}$. (This gives
a = 6.3, b = 0.063.) Clearly, D is a delayed and smoothed version of the input
derivative exactly as our derivative device would provide.

7 Conclusion and perspectives

We have presented a mathematical analysis of the core differentiation CRN intro-
duced by Whitby & al. [24]. In particular, we have shown that what is computed
is an approximation of the left derivative given a small time in the past with
a time constant determined by the diffusion constant between the input and
its internal representation: $\epsilon = \frac{1}{k_{\text{diff}}}$ . Moreover, there is a delay $\tau$ due to the
computation time that can also be precisely estimated given the rate of acti-
vation and degradation of the species encoding the derivative: $\tau = \frac{1}{k}$ . We have
shown that such results can be used in some cases to design error-correcting
terms and obtain excellent implementations of functions of input signals using
an approximation of their derivative on the fly.

From a synthetic biology perspective, the derivative CRN may be very rel-
evant in the context of biosensor design, when the test is not be about the
presence of some molecular compounds [6] but on their variation. A derivative
CRN is also needed to construct PID controllers. The derivative control is known
for damping the oscillations around the target of the controller but delays are
also known for producing such oscillations. Being able to determine and quantify
those delays and errors is thus important to optimize the design. This device may
also be used to approximate the derivative of an unknown external input in the
context of online cellular computing. Once again, delay may produce nefarious
artefacts that can easily be avoided when aware of the problem.

Furthermore, using the notion of SEPI to scan the biomodels database, we
were able to highlight a certain number of CRN models that contain the core
differentiation CRN. A high number of these models occur in models presenting
oscillations. We have shown on one such example, a circadian clock model, why
it makes sense for an oscillator to sense its own derivative, and to reproduce
what a mathematician would produce in a more direct way for the most basic
oscillatory function: sine and cosine.

Acknowledgment

This work benefited from ANR-20-CE48-0002 δifference project grant.
References

1. Emmanouil Alexis, Carolin CM Schulte, Luca Cardelli, and Antonis Papachristodoulou. Biomolecular mechanisms for signal differentiation. *IsScience*, 24(12):103462, 2021.

2. Sabine Becker-Weimann, Jana Wolf, Hanspeter Herzel, and Achim Kramer. Modeling feedback loops of the mammalian circadian oscillator. *Biophysical journal*, 87(5):3023–3034, November 2004.

3. Corentin Briat, Ankit Gupta, and Mustafa Khammash. Antithetic integral feedback ensures robust perfect adaptation in noisy biomolecular networks. *Cell Systems*, 2(1):15–26, Jan 2016.

4. Laurence Calzone, François Fages, and Sylvain Soliman. BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge. *Bioinformatics*, 22(14):1805–1807, 2006.

5. Michael Chevalier, Mariana Gómez-Schiavon, Andrew H. Ng, and Hana El-Samad. Design and analysis of a proportional-integral-derivative controller with biological molecules. *Cell Systems*, 9(4):338–353, 2019.

6. Alexis Courbet, Patrick Amar, François Fages, Eric Renard, and Franck Molina. Computer-aided biochemical programming of synthetic microreactors as diagnostic devices. *Molecular Systems Biology*, 14(4), 2018.

7. Péter Érdi and János Tóth. *Mathematical Models of Chemical Reactions: Theory and Applications of Deterministic and Stochastic Models*. Nonlinear science: theory and applications. Manchester University Press, 1989.

8. François Fages, Guillaume Le Guludec, Olivier Bournez, and Amaury Pouly. Strong Turing Completeness of Continuous Chemical Reaction Networks and Compilation of Mixed Analog-Digital Programs. In *CMSB’17: Proceedings of the fiveteen international conference on Computational Methods in Systems Biology*, volume 10545 of *Lecture Notes in Computer Science*, pages 108–127. Springer-Verlag, September 2017.

9. François Fages, Steven Gay, and Sylvain Soliman. Inferring reaction systems from ordinary differential equations. *Theoretical Computer Science*, 599:64–78, September 2015.

10. François Fages and Sylvain Soliman. Abstract interpretation and types for systems biology. *Theoretical Computer Science*, 403(1):52–70, 2008.

11. Martin Feinberg. Mathematical aspects of mass action kinetics. In L. Lapidus and N. R. Amundson, editors, *Chemical Reactor Theory: A Review*, chapter 1, pages 1–78. Prentice-Hall, 1977.

12. Steven Gay, François Fages, Thierry Martinez, Sylvain Soliman, and Christine Solnon. On the subgraph epimorphism problem. *Discrete Applied Mathematics*, 162:214–228, January 2014.

13. Steven Gay, Sylvain Soliman, and François Fages. A graphical method for reducing and relating models in systems biology. *Bioinformatics*, 26(18):i575–i581, 2010. special issue ECCB’10.

14. Mathieu Hemery and François Fages. Algebraic biochemistry: a framework for analog online computation in cells. In *CMSB’22: Proceedings of the twentieth international conference on Computational Methods in Systems Biology*, volume 13447 of *Lecture Notes in BioInformatics*. Springer-Verlag, September 2022.

15. Mathieu Hemery, François Fages, and Sylvain Soliman. On the complexity of quadratization for polynomial differential equations. In *CMSB’20: Proceedings of the eighteenth international conference on Computational Methods in Systems Biology*, Lecture Notes in BioInformatics. Springer-Verlag, September 2020.
16. Mathieu Hemery, François Fages, and Sylvain Soliman. Compiling elementary mathematical functions into finite chemical reaction networks via a polynomialization algorithm for ODEs. In *CMSB'21: Proceedings of the nineteenth international conference on Computational Methods in Systems Biology*, Lecture Notes in Bioinformatics. Springer-Verlag, September 2021.

17. Chi-Ying Huang and James E. Ferrell. Ultrasensitivity in the mitogen-activated protein kinase cascade. *PNAS*, 93(19):10078-10083, September 1996.

18. Hiroaki Kitano. Systems biology: A brief overview. *Science*, 295(5560):1662-1664, March 2002.

19. Nicolas le Novère, Benjamin Bornstein, Alexander Broicher, Mélanie Courtot, Marco Donizelli, Harish Dharuri, Lu Li, Herbert Sauro, Maria Schilstra, Bruce Shapiro, Jacky L. Snoep, and Michael Hucka. BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems. *Nucleic Acid Research*, 1(34):D689-D691, January 2006.

20. J. Myhill. A recursive function defined on a compact interval and having a continuous derivative that is not recursive. *Michigan Mathematical Journal*, 18(2):97-98, May 1971.

21. K. Oishi and E. Klavins. Biomolecular implementation of linear i/o systems. *IET Systems Biology*, 5(4):252-260, 2011.

22. C.E. Shannon. Mathematical theory of the differential analyser. *Journal of Mathematics and Physics*, 20:337-354, 1941.

23. Marko Vasic, David Soloveichik, and Sarfraz Khurshid. CRN++: Molecular programming language. In *Proc. DNA Computing and Molecular Programming*, volume 11145 of LNCS, pages 1-18. Springer-Verlag, 2018.

24. Max Whitby, Luca Cardelli, Marta Kwiatkowska, Luca Laurenti, Mirco Tribastone, and Max Tschaikowski. Pid control of biochemical reaction networks. *IEEE Transactions on Automatic Control*, 67(2):1023-1030, 2021.

25. Guang Yao, Tae Jun Lee, Seiichi Mori, Joseph R Nevins, and Lingchong You. A bistable rb-e2f switch underlies the restriction point. *Nature cell biology*, 10(4):476-482, 2008.

## Appendix: computation of integration with a delay

To prove that the drift of the output is a direct consequence of the delay, we first compute the input and the approximate derivative for our choice of input:

$$ 
\begin{aligned}
x(t) &= 1 + \sin(t) \\
x'(t - \tau) &= \cos(t - \tau) \\
&= \cos(t) + \tau \sin(t) + o(\tau^2)
\end{aligned} 
\qquad (18) 
$$

Then we can compute the output up to the first order:

$$ 
\begin{aligned}
y(t) &= \int 2x(s)x'(s - \tau)ds \\
&= \int 2(1 + \sin(s)) \cos(s)ds + \int 2\tau(\sin(s) + \sin^2(s))ds \\
&= (1 + \sin(t))^2 + 2\tau \int \sin(s) + \sin^2(s)ds \\
y(t) &\approx (1 + \sin(t))^2 + \tau t
\end{aligned} 
\qquad (19) 
$$
Then, to correct the observed drift, we propose to introduce a delay signal and use it in the computation to produce the output species $Y_+$ and $Y_-$, with the following CRN:

$$
\begin{align}
X_{\text{ext}} &\xrightarrow{k_{\text{diff}}} X_{\text{in}}, & X_{\text{in}} &\xrightarrow{k_{\text{diff}}} X_{\text{ext}}, \nonumber \\
X_{\text{in}} &\xrightarrow{k} X_{\text{in}} + X_{\text{delay}}, & X_{\text{delay}} &\xrightarrow{k} \emptyset, \nonumber \\
X_{\text{ext}} &\xrightarrow{k_{\text{diff}}k} X_{\text{ext}} + D_{+} & D_{+} &\xrightarrow{k} \emptyset \tag{20} \\
X_{\text{in}} &\xrightarrow{k_{\text{diff}}k} X_{\text{in}} + D_{-} & D_{-} &\xrightarrow{k} \emptyset \nonumber \\
X_{\text{delay}} + D_{+} &\xrightarrow{2} X_{\text{delay}} + D_{+} + Y_{+} & X_{\text{delay}} + D_{-} &\xrightarrow{2} X_{\text{delay}} + D_{-} + Y_{-} \nonumber \\
D_{+} + D_{-} &\xrightarrow{\text{fast}} \emptyset & Y_{+} + Y_{-} &\xrightarrow{\text{fast}} \emptyset \nonumber
\end{align}
$$