From Electric Circuits to Chemical Networks

TL;DR

This paper presents a systematic method to convert linear electric circuits into equivalent chemical reaction networks, enabling the use of electronic design principles in biochemical system engineering.

Contribution

It introduces a novel approach to translate electric circuit structures into chemical networks, bridging electronics and biochemistry for system design.

Findings

Chemical networks can replicate electric circuit functionalities.

The method allows leveraging electronic device libraries for biochemical network design.

Chemical networks derived are intelligible and comparable to natural biochemical systems.

Abstract

Electric circuits manipulate electric charge and magnetic flux via a small set of discrete components to implement useful functionality over continuous time-varying signals represented by currents and voltages. Much of the same functionality is useful to biological organisms, where it is implemented by a completely different set of discrete components (typically proteins) and signal representations (typically via concentrations). We describe how to take a linear electric circuit and systematically convert it to a chemical reaction network of the same functionality, as a dynamical system. Both the structure and the components of the electric circuit are dissolved in the process, but the resulting chemical network is intelligible. This approach provides access to a large library of well-studied devices, from analog electronics, whose chemical network realization can be compared to natural biochemical networks, or used to engineer synthetic biochemical networks.

Figures

Equations

$$E{\partial }_{t}x=Ax+Bu$$

$$\left(\begin{array}{c}{\partial }_{t}u\\ {\partial }_{t}z\end{array}\right)=D\cdot \left(\begin{array}{c}u\\ z\end{array}\right)+d$$

$${\partial }_{t}i$$

$$iR$$

$${\partial }_{t}i$$

$$v_{out}$$

$$v_{in}$$

$${\partial }_{t}i$$

$$0$$

$$\underset{E}{\underbrace{\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)}}\underset{{\partial }_{t}x}{\underbrace{\left(\begin{array}{c}{\partial }_{t}i\\ {\partial }_t{v}_{out}\end{array}\right)}}$$

$$F_h\left(\begin{array}{c}i\\ v_{out}\end{array}\right)=\left(\begin{array}{c}\frac{v_{in}-i}{1+h}\\ \frac{v_{in}-i-(1+h)v_{out}}{h(1+h)}\end{array}\right)$$

$${\partial }_{t}i$$

$${\partial }_t{v}_{out}$$

$$i$$

$${\partial }_t{i}^{+}$$

$${\partial }_t{v}_{out}^{+}$$

$$v_{in}^{+}$$

$$v_{in}^{-}$$

$$v_{in}^{+}$$

$$v_{in}^{-}$$

$$v_{out}^{-}$$

$$i^{+}+i^{-}$$

$$E{\partial }_{t}x=Ax+b,\text{with }E,A\in {\mathbb{R}}^{n\times n},b\in {\mathbb{R}}^n,$$

$${\partial }_{t}x=E^{-1}Ax+E^{-1}b.$$

$$x[i+1]=x[i]+h{F}_h(x[i]),\text{with }{F}_h(x):=(E-hA{)}^{-1}(Ax+b).$$

$${\partial }_t{x}_h=(E-hA{)}^{-1}\left(A{x}_h+B{u}_h^{⟨0⟩}+hB{u}_h^{⟨1⟩}\right),$$

$$\left(\begin{array}{c}{\partial }_t{u}_h^{⟨0⟩}\\ {\partial }_t{z}_h^{⟨0⟩}\end{array}\right)$$

$$\left(\begin{array}{c}{\partial }_t{u}_h^{⟨1⟩}\\ {\partial }_t{z}_h^{⟨1⟩}\end{array}\right)$$

$$u_h^{⟨0⟩}(0)$$

$$z_h^{⟨0⟩}(0)$$

$$\begin{array}{rl}{\partial }_t{x}^{+}& ={\hat{A}}^{+}{x}^{+}+{\hat{A}}^{-}{x}^{-}+{\hat{b}}^{+}\\ {\partial }_t{x}^{-}& ={\hat{A}}^{+}{x}^{-}+{\hat{A}}^{-}{x}^{+}+{\hat{b}}^{-}\end{array}$$

$$\begin{array}{rl}{\partial }_t{x}^{+}& ={\hat{A}}^{+}{x}^{+}+{\hat{A}}^{-}{x}^{-}+{\hat{b}}^{+}-\gamma Q(x^{+},x^{-})\\ {\partial }_t{x}^{-}& ={\hat{A}}^{+}{x}^{-}+{\hat{A}}^{-}{x}^{+}+{\hat{b}}^{-}-\gamma Q(x^{+},x^{-})\end{array}$$

$$\mathcal{R}$$

$$\cup \{x_j^{+}\text{xlongrightarrow}{\hat{A}}_{i,j}^{+}{x}_j^{+}+x_i^{+}\mid i,j\}\cup \{x_j^{-}\text{xlongrightarrow}{\hat{A}}_{i,j}^{-}{x}_j^{-}+x_i^{+}\mid i,j\}$$

$$\cup \{x_j^{+}\text{xlongrightarrow}{\hat{A}}_{i,j}^{-}{x}_j^{+}+x_i^{-}\mid i,j\}\cup \{x_j^{-}\text{xlongrightarrow}{\hat{A}}_{i,j}^{+}{x}_j^{-}+x_i^{-}\mid i,j\}$$

$$u(t)$$

$${\partial }_{t}u$$

Full text

I Introduction

Living organisms perform a variety of functions that can be described in abstract terms as information processing and regulation. Analogies have been drawn between the biochemical reaction networks that perform such functions and electric circuits of similar nature [1, 2]. The comparison is useful in systems biology, in trying to understand the function of natural systems, and in synthetic biology, in trying to engineer desired functionality in a biochemical context.

An obstacle to exploiting this analogy is that the fundamental components of biochemistry and electronics are very different: the phenomena of resistance, induction, and capacitance based on the interplay between electric and magnetic fields have no immediate parallel in terms of chemical concentrations. Hence, it is not clear how to take systematic advantage of engineering knowledge developed in electronics in understanding biochemical systems. Instead, the search for components of biological network has progressed in different and certainly more appropriate directions [3, 4].

Even within electronics, though, the precise nature of electric components is incidental to the desired function. For example, one may wish to filter the high frequencies of a signal represented by an oscillating voltage. There are countless electric circuits that can perform this function, based on different classes of components in different configurations. Some of those circuits are based on operational amplifiers, which are themselves built from a large network of components to perform an abstract function to which they are incidental. The nature of the fundamental components is inessential, as long as they can be combined to provide a wide range of functionality.

A common way to describe essential function, both in electronics and in biochemistry, is through a system of ordinary differential equations (ODEs). Once the function of a circuit is reduced to this form, it does not matter if the quantities represented are voltages or concentrations: it only matters the way in which they vary over time. Conversely, given an ODE system, one may ask the engineering question of how to realize a circuit (electronic or biochemical) that can perform that function. An early example of this inverse process is the synthesis of mechanical and then electric analog circuits from differential equations [5, 6].

Certain classes of polynomial ODE systems can be systematically turned into chemical reaction networks (CRNs) that obey the same kinetics [7]. Further, CRNs can be compiled systematically into a collection of molecules that can be engineered to obey the kinetics of those reactions [8]. In this paper we wish to go one step further on the front end of this process. Taking advantage of the large libraries of known electric circuits, we wish to take an arbitrary (but linear, for now) electric circuit and show how to turn it into a set of molecules that obey the kinetics of the quantities in that circuit.

The first obstacle we need to confront is that even common electric circuits describe behaviors that go beyond ODEs. Algorithmic approaches for the analysis of linear circuits such as Modified Nodal Analysis produce, in general, systems of differential algebraic equations (DAEs), where the algebraic equations are induced by classical node analysis based on Kirchhoff laws [9]. Hence we first reduce DAEs to ODEs, after which we can apply some further techniques. The second obstacle is to take an ODE that may be about voltages and currents, and turn it into a form that can be interpreted chemically. This means that each variable should only take nonnegative quantities (for concentrations), and that appropriate chemical reactions should be derived about those quantities.

This approach has a dual purpose. From a systems biology point of view, we may want to compare the CRNs derived from electric circuits to the ones occurring in nature. This might help elucidate the function of natural networks. Or, at least, it will provide a spectrum of possible chemical networks of known function, whose structure may not be obvious, therefore broadening our expectations of what is possible chemically. From a systems biology point of view, we take it as given that an abstract CRN can be turned into a collection of molecules. In fact, multiple target architectures are possible, from short oligonucleotides in solution [8] to gene networks [10]. The ability to generate molecular configurations from a vast existing library of (electric, or other) circuits is appealing for systematizing the generations of synthetic organisms. In extending the known techniques from ODEs to DAEs, we extend the scope of potential libraries we can draw from.

Contributions

Our main result is a systematic technique that transforms linear DAE systems into CRNs. Our technique is, to the best of our knowledge, novel. DAE systems are either solved symbolically by relying on index reduction [11], or numerically by relying on numerical methods that compute the trajectories for a given initial condition [12]. In contrast to our approach, index reduction introduces additional derivatives of signals that may not be available to the circuit, while numerical methods do not transform the DAEs into ODEs, which seems necessary for transforming DAE systems into CRNs. We analyze, as an example, an electric high pass filter and provide a chemical reaction network for it, whose function and architecture can be independently interpreted in a biological context.

II Outline of Methods

We start from a linear electric circuit composed of resistors, capacitors and inductors, with variables ranging over voltages and currents, and we systematically derive an equivalent CRN where chemical species (or more precisely their differences) approximate the trajectories of the original variables.

The basis for this process is the so-called Hungarian Lemma [7], which provides a method for converting certain polynomial ODEs into CRNs by converting each monomial on the right hand side of a differential equation into a separate chemical reaction. Polynomial ODEs can represent, exactly, a much broader class of ODEs, including fractional, trigonometric, and exponential terms [13], thus covering a broad range of chemical behavior including Hill kinetics [14].

The Hungarian Lemma, however, has specific requirements. First, the concentrations of the chemical species must be nonnegative, while ordinary ODE variables, and in particular voltages and currents, may be negative. Second, if a monomial has a negative sign, then the differential variable on the left-hand side of the equation must appear as a factor in the monomial. This means, for example, that the ODE $\partial_{t}x=y$ (where the growth rate of $x$ is given by the concentration of $y$, with $\partial_{t}x$ denoting the time derivative of $x$) can be reduced to the reaction $y\rightarrow x+y$. And the ODE $\partial_{t}x=-xy$ can be reduced to the reaction $x+y\rightarrow y$. But the ODE $\partial_{t}x=-y$ (where the decrease in rate of $x$ is given by the concentration of $y$) cannot be reduced: for $x$ to decrease it must appear on the left-hand-side of a chemical reaction, which implies it should appear in a monomial for $\partial_{t}x$ by the law of mass action. An ODE system with no such forbidden negative monomial is said to be Hungarian, and any Hungarian ODE system can be reduced to a CRN (although not uniquely) whose mass action kinetics reproduces the original ODE. We use a technique to reduce a polynomial non-Hungarian ODE to an Hungarian one in twice as many variables, thus allowing us to produce CRNs also for non-Hungarian ODEs. In the same step, we make all trajectories nonnegative so that they can be realized by chemical species.

Some simple electric circuits yield ODE systems that can be converted to CRNs as outlined. Pure resistor circuits yield simple algebraic equations. But more complex circuits yield general DAEs [9], which we must be prepared to handle. Our main technique applies to linear DAE systems of the form

[TABLE]

Here, $x\in\mathbb{R}^{n}$ is the (column) vector of dependent variables, $E\in\mathbb{R}^{n\times n}$ produces a linear combination of their derivatives, $A\in\mathbb{R}^{n\times n}$ produces a linear combination of the variables, and the term $B\in\mathbb{R}^{n\times m}$ is the input matrix, and $u\in\mathbb{R}^{m}$ is the vector of inputs, such as voltage or current sources.

In general, inputs are assumed to be arbitrary, known functions of time. However, for the purposes of converting electric circuits into CRNs, it is necessary to appropriately encode also the inputs as chemical species. In this paper, we will assume that the input vector $u$ can be itself described as the solution of a system of equations. Specifically, we assume that it satisfies an affine ODE system of the form

[TABLE]

for some matrix $D\in\mathbb{R}^{(m+k)\times(m+k)}$ and $\left(u(0),z(0)\right)^{T},d\in\mathbb{R}^{m+k}$. Intuitively, the input $u$ in (1) is part of an ODE solution which may depend on auxiliary ODE variables $z$ that do not appear in the DAE. This is a rather general setting that allows us to encode arbitrary time-varying inputs by approximating them with Fourier series, which can be expressed as solutions of linear ODEs (see Section IV).

Our technique converts the overall system (1)-(2) into an ODE system over the same variables, up to a controllable approximation. That ODE system can then be transformed into a CRN as discussed above.

We now describe in detail the entire process of converting linear electric circuits to CRNs, through a small textbook example involving a single differential equation and a single algebraic equation for the well-known $RL$ (resistor-inductor) circuit in Figure 1 [15]. The analysis of the circuit proceeds as follows. We let $v_{\mathit{in}}$ denote the input voltage (measured with respect to the ground node); the output is the voltage $v_{\mathit{out}}$ across the inductor. By standard node analysis, using Kirchhoff’s current law at each node, we obtain that the three currents are equal, so $i\triangleq i_{T}=i_{R}=i_{L}$. Faraday’s law for the inductor $L$, and Ohm’s law for the resistor $R$, then give us:

[TABLE]

This is a DAE, where (3a) is a differential equation, and (3b) is an algebraic equation. To find its solution, we can replace $v_{\mathit{out}}=v_{\mathit{in}}-iR$ from (3b) in (3a), obtaining (4a), which can be integrated to obtain $i(t)$. From that solution we can then obtain $v_{\mathit{out}}(t)$ via (4b).

[TABLE]

If we briefly assume that $v_{\mathit{in}}$ is a non-negative input source, then (4a) can be converted to a mass action CRN as follows (using chemical species with the same names as our variables):

[TABLE]

where with the symbol $\emptyset$ we have denoted the empty set of products. These two chemical reactions yield (4a) for the evolution of the concentration of species $i$. However, a voltage $v_{\mathit{in}}$ may be negative, which cannot be modeled with chemical concentrations. Additionally, we are interested in the output $v_{\mathit{out}}$, not $i$, and therefore we need to find a way to realize chemically equation (4b) as well.

Because of those difficulties, we cannot make much progress without a more general technique to implement DAEs as CRNs. Our technique involves an approximation, like the one that seems necessary just for algebraic equations, but it can be used in the general case of linear DAEs. We now illustrate it by applying it to the example of Figure 1. We first rearrange our DAE as (5a,5b). Setting $x=(i,v_{\mathit{out}})^{T}$ we have:

[TABLE]

which can be arranged into the form (1):

[TABLE]

where we have fixed the constant vector $b:=Bu$ under the assumption of a constant input source $u:=v_{\mathit{in}}$.

We approximate the DAE, for a parameter $h>0$, by symbolically computing $F_{h}(x)=(E-hA)^{-1}(Ax+b)$, which is related to the numerical backward Euler method. Taking $R=L=1$ for simplicity, we obtain:

[TABLE]

We next use $F_{h}(x)$ as the right-hand side of a new ODE system $\partial_{t}x=F_{h}(x)$, which is such that for $h\rightarrow 0$ the solution of the ODE system (7a)-(7b) converges to the solution of the original DAE system (5a)-(5b) (see Theorem 1).

[TABLE]

Indeed, we can easily see that for $h\rightarrow 0$, (7a) converges exactly to the ODE (4a). As for (7b), this is now a differential equation approximating equation (5b) for $h\rightarrow 0$, where we notice that a small value of $h$ makes (7b) evolve much faster than (7a).

We have reduced a DAE to an ODE, but (7a)-(7b) is not Hungarian because of the $-i$ monomial in (7b). Keeping in mind that we need to deal eventually with non-Hungarian ODEs, we now apply a positivation technique in the style of Oishi and Klavins [16], where each variable is represented as the difference of two non-negative variables:

[TABLE]

Let us now abbreviate $p=1/(1+h)$, $q=1/(h+h^{2})$, and $r=1/h$, and consider the ODE system where we separate the positive and negative monomials of each ODE in (7a)-(7b) into two ODEs:

[TABLE]

The initial conditions for this new system must satisfy $i^{+}_{0}-i^{-}_{0}=i_{0}$ with $i^{+}_{0},i^{-}_{0}\geq 0$, etc. Since differentiation is a linear operator, the solutions of (7a)-(7b) can be recovered as differences from the solutions of (8a)-(8b): $\partial_{t}i^{+}-\partial_{t}i^{-}=\partial_{t}i$ and $\partial_{t}v_{\mathit{out}}^{+}-\partial_{t}v_{\mathit{out}}^{-}=\partial_{t}v_{\mathit{out}}$. Although the goal was to make all variables non-negative, we now also have a Hungarian ODE system because all the monomials in (8a)-(8b) are positive. Hence there is no further difficulty in converting these ODEs to mass action reactions, obtaining the following linear CRN with one reaction for each monomial in (8a)-(8b), and with the parameter $h$ appearing in the reaction rates:

[TABLE]

Here the input $v_{\mathit{in}}^{\pm}$ always acts as a simple catalyst. We note that the CRN implementation does not depend on the actual value of $v_{\mathit{in}}$, which only affects the initial condition of the chemical species that represent $v_{\mathit{in}}^{\pm}$. This decoupling between the CRN implementation of the circuit and that of the input sources carries over to the more general case when the sources are time-varying solutions of the ODEs (2), see Theorem 3 and the subsequent discussion.

Inspecting (9), the chemical species $v_{\mathit{out}}^{\pm}$ and $i^{\pm}$ are involved in autocatalytic cycles. For example both $i^{+}$ and $i^{-}$ grow exponentially over time, while their difference $i$ remains bounded. It is possible to eliminate such exponential growths by adding non-linear dampening reactions to the otherwise linear CRN:

[TABLE]

The first reaction, for example, preserves the difference $i^{+}-i^{-}$, and results in two new identical monomials in the ODEs for $i^{+}$ and $i^{-}$, that then cancel in $\partial_{t}i^{+}-\partial_{t}i^{-}$. Hence that reaction does not change the $i$ solution, but keeps $i^{+}$ and $i^{-}$ bounded.

The network consisting of reactions (9)-(10) is depicted in Figure 2, where for small $h$ we have $1\approx p\ll q\approx r$, and we can take $\gamma=r$. This network has the flavor of an incoherent feedforward motif [4], considering parallel pairs $x^{\pm}\to y^{\pm}$ as activations and cross pairs $x^{\pm}\to y^{\mp}$ as inhibitions, thereby $v_{\mathit{in}}$ activates both $i$ and $v_{\mathit{out}}$, and $i$ ‘incoherently’ inhibits $v_{\mathit{out}}$. Additionally, the motif of mutual catalysis and join degradation around $i^{\pm}$ makes that pair stabilize to a copy of its input $v_{\mathit{in}}^{\pm}$ (in the sense that at steady state $i^{+}-i^{-}=v_{\mathit{in}}^{+}-v_{\mathit{in}}^{-}$) regardless of the value of the rate $p$. This motif is repeated around $v_{\mathit{out}}^{\pm}$. When the input $v_{\mathit{in}}^{\pm}$ remains constant, $i^{\pm}$ becomes a copy of $v_{\mathit{in}}^{\pm}$, and $v_{\mathit{out}}^{\pm}$ becomes a copy of the sum of its two opposite inputs, $v_{\mathit{in}}^{\pm}$ and $i^{\mp}$, and so it converges to a baseline output of $v_{\mathit{out}}^{+}-v_{\mathit{out}}^{-}\approx 0$. When the input $v_{\mathit{in}}$ changes, it affects $v_{\mathit{out}}$ quickly and $i$ slowly, with a delayed inhibition of $v_{\mathit{out}}$ by $i$. It has been shown that feedforward motifs can behave like high-pass filters [17].

The subnetwork in Figure 2 consisting of $v_{\mathit{in}}^{\pm}$, $v_{\mathit{out}}^{\pm}$, and the connecting $q$,$r$,$\gamma$ arcs, is in itself also a low-pass filter. It is exactly what is obtained when replacing the inductor with a capacitor of capacitance $C$ in Figure 1, yielding a well known low-pass filter, and deriving the CRN from it by positivation (with $q=r=1/RC$). The process is simpler in this case, since a single ODE is generated from that circuit, and no approximation via $h\rightarrow 0$ is required.

In summary, we have derived an intelligible chemical reaction network from an electric circuit, and we are guaranteed that it implements the same functionality, as shown in the next section.

III Methods

III-A From DAEs to ODEs

We first show how to convert a linear DAE system into a linear ODE system. To this end we start with a DAE system in the form

[TABLE]

which corresponds to (1) under the assumption of constant inputs $u$. From now on we focus on regular DAE systems, i.e., systems for which there exists no initial condition $x(0)\in\mathbb{R}^{n}$ for which they admit more than one solution. We let $\mathcal{D}$ denote the set of initial conditions for which a regular DAE admits solutions. Elements of $\mathcal{D}$ are called consistent initial conditions (see, e.g., [12, Section 2.1] for details). Obviously, any physically meaningful DAE model has to be regular. In the case of electric circuits, for instance, non-regular DAE systems may arise in the presence of short circuits and other erroneous designs.

If $E$ is invertible, this DAE can be directly recast into a linear ODE system via

[TABLE]

However, in the case of linear electric circuits $E$ is in general not invertible. In this case, a transformation of a DAE system into an ODE system requires index reduction [11, 12], which relies on expensive symbolic computations. Our method consists in circumventing this analysis by considering an explicit scheme arising from numerical methods for the solution of DAE systems.

A numerical method is an algorithm that generates, for a given small time step $h>0$ and initial condition $x(0)$, a sequence $(x[i])_{i}$ such that $x[i]\approx x(ih)$, where $x:[0;\infty)\to\mathbb{R}^{n}$ denotes the true solution of (11). Numerical methods are guaranteed to converge to the true solution $x$ when $h$ approaches zero. That is, that for any error threshold $\varepsilon>0$ and finite time horizon $T>0$, one can find a sufficiently small time step $h>0$ such that $\max_{0\leq i\leq N}\lVert x[i]-x(ih)\rVert\leq\varepsilon$ and $T/h=N\in\mathbb{N}$.

A common numerical method for the solution of DAE systems is the backward Euler method [12, Section 5.2] which is given by

[TABLE]

The function $F_{h}$ is well-defined for sufficiently small $h>0$ because the DAE system is regular [12, Section 2.1].

Noting that the “slope” of the Euler method at point $x[i]$, $(x[i+1]-x[i])/h$, is given by $F_{h}(x[i])$, we expect that the Euler sequence $(x[i])_{i\geq 0}$ will match the solution of the ODE system $\dot{x}_{h}=F_{h}(x_{h})$. That is, we expect that $x_{h}(ih)\approx x[i]$ for all $0\leq i\leq N$. Then, the convergence of the Euler sequence to the DAE solution $x$ would allow us to conclude that $x_{h}(ih)\approx x[i]\approx x(ih)$.

The next theorem is our first main result and proves that this is indeed the case.

Theorem 1**.**

For any $\varepsilon>0$, $x(0)\in\mathcal{D}$ and $T>0$, there exists an $h>0$ such that $\sup_{0\leq t\leq T}\lVert x(t)-x_{h}(t)\rVert\leq\varepsilon$, where $\partial_{t}x_{h}=F_{h}(x_{h})$ and $x_{h}(0)=x(0)$.

We now extend Theorem 1 to systems (1)-(2).

Theorem 2**.**

Given a DAE system via (1) and (2), consider the ODE system

[TABLE]

where $h>0$ is small and the functions $u_{h}^{\langle 0\rangle},u_{h}^{\langle 1\rangle}\in\mathbb{R}^{m}$ satisfy

[TABLE]

with initial conditions

[TABLE]

Then, for any $\varepsilon>0$, $x(0)\in\mathcal{D}$ and $T>0$, there exists an $h>0$ such that $\sup_{0\leq t\leq T}\lVert x(t)-x_{h}(t)\rVert\leq\varepsilon$ if $x_{h}(0)=x(0)$.

Note that $I-hD$ is strictly diagonal dominant and therefore invertible for sufficiently small values of $h$. It can be shown that $(u^{\langle 0\rangle}_{h},v^{\langle 0\rangle}_{h})^{T}$ converges to $(u,v)^{T}$ from (2) as $h\to 0$. Hence, Theorem 2 essentially replaces the constant vector $b$ in $F_{h}(x)=(E-hA)^{-1}(Ax+b)$ by the function $Bu$.

III-B From ODEs to CRNs

We next present a technique that transforms the ODE approximation from Section III-A into a CRN. The approach borrows ideas from [16] that transforms linear ODE systems (i.e., systems without algebraic constraints) into CRNs. We wish to point out, however, that our approach considers state space representation, while [16] acts on the frequency domain.

In the following, let $\partial_{t}x=\hat{A}x+\hat{b}$ denote some ODE system with initial condition $x(0)$.

Proposition 1**.**

Any non-negative quadruple $(\hat{A}^{+},\hat{A}^{-},\hat{b}^{+},\hat{b}^{-})$ satisfying $\hat{A}=\hat{A}^{+}-\hat{A}^{-}$ and $\hat{b}=\hat{b}^{+}-\hat{b}^{-}$ induces the positivation

[TABLE]

If (15) is subject to non-negative $x^{+}(0),x^{-}(0)\in\mathbb{R}^{n}_{\geq 0}$ with $x(0)=x^{+}(0)-x^{-}(0)$, the solution $(x^{+},x^{-})$ remains non-negative and satisfies $x=x^{+}-x^{-}$.

While positivations trivially satisfy the properties of the Hungarian lemma discussed in Section II and can therefore be readily translated into CRNs, they may exhibit divergence even if the original system is bounded, see for instance (8a), which implies that $i^{+}$ and $i^{-}$ diverge. Fortunately, one can apply a correction that leads to bounded positivations.

Proposition 2**.**

Given a positivation $(\hat{A}^{+},\hat{A}^{-},\hat{b}^{+},\hat{b}^{-})$, define the quadratic function $Q(x^{+},x^{-})=(x_{1}^{+}x_{1}^{-},\ldots,x_{n}^{+}x_{n}^{-})^{T}$. Then, for any $\gamma>0$, the corresponding Hungarization is

[TABLE]

If (16) is subject to non-negative $x^{+}(0),x^{-}(0)\in\mathbb{R}^{n}_{\geq 0}$ with $x(0)=x^{+}(0)-x^{-}(0)$, ODE system (16) admits a non-negative solution on $[0;\infty)$ that satisfies $x=x^{+}-x^{-}$. Moreover, if $x$ is bounded on $[0;\infty)$, then so is $(x^{+},x^{-})$.

By applying the law of mass action, it can be easily seen that the $Q$ terms in (16) are captured by the annihilation reactions $x^{+}_{1}+x^{-}_{1}\xlongrightarrow{\gamma}\emptyset,\ldots,x^{+}_{n}+x^{-}_{n}\xlongrightarrow{\gamma}\emptyset$. Combining this with Proposition 2, we arrive at the following statement.

Proposition 3**.**

Define the chemical reactions of Hungarization (16) as

[TABLE]

Then, reactions $\mathcal{R}$ induce, via the law of mass action, the ODE system (16).

For instance, if the positivation is given by (7a-7b), then (9-10) constitute $\mathcal{R}$.

Proposition 3 and Theorem 2 yield our main result.

Theorem 3**.**

Given a DAE system via (1) and (2), let

• •

$\mathcal{H}_{1,h}$* and $\mathcal{H}_{2,h}$ denote the Hungarization of (12) and (13-14), respectively.*

• •

$\mathcal{R}_{1,h}$* and $\mathcal{R}_{2,h}$ refer to the chemical reactions of $\mathcal{H}_{1,h}$ and $\mathcal{H}_{2,h}$, respectively.*

Then, the following holds true.

1. a)

The solution of the union CRN given by $\mathcal{R}_{1,h}\cup\mathcal{R}_{2,h}$ converges to the DAE solution as $h\to 0$. 2. b)

A change of $D$ and $d$ affects the reactions $\mathcal{R}_{2,h}$ but does not alter the reactions $\mathcal{R}_{1,h}$.

As anticipated in Section II, Theorem 3 ensures that a) our encoding is correct up to a controllable error and b) that the CRN implementation of the circuit, $\mathcal{R}_{1,h}$, does not depend on the CRN implementation of the input, $\mathcal{R}_{2,h}$.

Theorems 2 and 3 allow for composition of circuits: a circuit expressed as a DAE (1), with input provided by and ODE (2), yields another ODE (12) that can be supplied as input to a further circuit. The corresponding CRNs can be composed as well.

IV Methods Applied to Example

The $RL$ circuit discussed in Section II is a high-pass filter, attenuating the low frequencies of the input while transmitting the high frequencies to the output. The cutoff frequency is the frequency at which the input signal is attenuated by $\frac{1}{2}$ its power, or equivalently its amplitude is attenuated by $-3\text{dB}\approx\sqrt{\frac{1}{2}}\approx 0.707$. In our circuit, the cutoff frequency is $f_{c}=\frac{R}{2\pi L}\text{Hz}$, where $R$ is the value of the resistance (measured in ohm) and $L$ is the inductance (in henry).

Here we show how we can apply Theorem 3 to the $RL$ circuit with a time-varying input represented by an arbitrary differentiable function. Such a function can be approximated arbitrarily well by a Fourier series $u(t)$ given by

[TABLE]

where $\alpha$, $\beta_{i}$, $\omega_{i}$, and $\gamma_{i}$ are constant parameters. It can be seen that $u(t)$ can be written as the solution of the ODE system:

[TABLE]

with initial conditions $z_{i}(0)=\sin(\gamma_{i})$, $\bar{z}_{i}(0)=\cos(\gamma_{i})$ and $u(0)=\alpha+\sum_{i=1}^{N}\beta_{i}\sin(\gamma_{i})$ with $1\leq i\leq N$, whereby $z_{i}$ and $\bar{z}_{i}$ are auxiliary variables whose solutions give the sinusoidal components of the series. We can recognize the system (17) to be in the required form (2).

As a first example of an input waveform, we consider the ODE system

[TABLE]

With initial conditions $u(0)=0$ and $z(0)=1$, this yields the solution $u(t)=\sin(t)$ for all $t$, whose frequency is the cutoff frequency. Figure 3 shows simulations of the $\sin(t)$ CRN composed with the high-pass filter CRN, taking $h=0.01$ for a sufficiently good approximation, and $\gamma=r$ for a sufficiently fast degradation. As expected, the output $v_{\mathit{out}}$ is attenuated by $-3\text{dB}\approx 0.707$, and its phase is shifted by $45^{\circ}$. The variation of the underlying non-negative variables is shown on the right.

As a second example, in a biological context, high pass filters exhibit perfect adaptation [18]: they adapt to slow but possibly large variations in input stimulus and still react to quick changes. In Figure 4 we supply stepped inputs via an appropriate Fourier series to our filter. At each sudden increase or decrease, the output reacts quickly and then settles back to its original level. The size of the transient response is proportional to the step size, but independent of the level of the input. The adaptation level can be set to any level, not just zero, by adding a constant contribution to $v_{\mathit{out}}^{+}$, so that the output can represent the (positive) concentration of a certain protein.

V Discussion

We have presented a method to convert linear DAEs to CRNs which hinges on a transformation into an approximate linear ODE system with arbitrary accuracy. This is then translated into a set of reactions where the time-course evolutions of the concentrations of the chemical species can be directly related to the original DAE solution.

In principle, any DAE system can be exactly transformed into an ODE system by means of so-called index reduction [12]. However, this relies on symbolic computations. Moreover, the so-obtained ODE system will contain derivatives of the input signal, thus requiring for additional approximations. For instance, by differentiating (3b) and using (3a) we obtain $\partial_{t}v_{\mathit{out}}=\partial_{t}v_{\mathit{in}}-\frac{R}{L}v_{\mathit{out}}$, which together with (3a) is a simple ODE system in the dependent variables with no algebraic equations. This system now depends on the derivative of the input, $\partial_{t}v_{\mathit{in}}$, and would have to be combined with another circuit to supply that derivative from the given input $v_{\mathit{in}}$. In contrast to index reduction our technique does not requires input derivatives. Moreover, while index reduction requires symbolic computations, the matrix inversion at the basis of the construction of the approximate linear ODE system can also be performed using numerical techniques.

The precision of the linear ODE depends on a parameter which, when taken asymptotically small, has the effect of rapidly equilibrating certain components of the ODE system. Therefore, in this respect our approach can be related to quasi steady-state approximation (QSSA, [19, 20]), which applies to semi-explicit DAEs in the special form

[TABLE]

Essentially, QSSA replaces the algebraic constraints $0=A_{2}y+B_{2}u$ with $\varepsilon\partial_{t}y=A_{2}y+B_{2}u$ for some $\varepsilon\approx 0$. However, DAEs of electric circuits are not semi-explicit in general. For instance, the circuit given in Figure 5 yields the DAE system

[TABLE]

which is not semi-explicit because it has two differential variables on the left-hand side.

Extensions of this work are needed to tackle non-linear DAE systems arising from non-linear electronic components such as diodes and transistors. The conversion of polynomial ODEs to Hungarian and positive ones, and thus to CNRs, works essentially unchanged also for non-linear polynomial systems, and further extends to ODE systems including trigonometric and exponential functions, which can model transistors. However, this must be coupled with a general method for converting non-linear DAEs arising from electronic circuits to ODEs.

Proof of Theorem 1 and Theorem 2

Before proving Theorem 1 and 2, we first have to establish some auxiliary results. To allow for a compact notation, we denote in the present section the $i$-th step of the numeric sequence by $x_{i}$ rather than $x[i]$.

Proposition 4**.**

Consider the ODE systems $\partial_{t}{x}=F(x)$ and $\partial_{t}{x}_{h}=F_{h}(x)$ where $F$ and $F_{h}$ are assumed to be Lipschitz continuous on some bounded domain $B\subseteq\mathbb{R}^{n}$ and $L\geq 0$ denotes the Lipschitz constant of $F$. Let us assume further that both ODE systems have solutions on $[0;T]$ which remain in $B$ and that $\sup\{\lVert F(x)-F_{h}(x)\rVert\mid x\in B\}\leq\eta$. Then, if $x(0)=x_{h}(0)$, for all $0\leq t\leq T$ it holds that

[TABLE]

Proof.

We first show a modified version of Gronwall’s inequality. To be more specific, let $\xi_{1}$ and $\xi_{2}$ be positive constants and $v$ a continuous function on $0\leq t<\infty$ such that

[TABLE]

Then, it holds that $v(t)\leq\tfrac{\xi_{2}}{\xi_{1}}(\mathrm{e}^{\xi_{1}t}-1)$. To see this, we first rewrite (20) to

[TABLE]

Since this rewrites to $\tilde{v}(t)\leq\tilde{\alpha}+\int_{0}^{t}\tilde{v}(s)\tilde{w}(s)ds$ for $\tilde{v}(s):=v(s)+\tfrac{\xi_{2}}{\xi_{1}}$, $\tilde{\alpha}:=\tfrac{\xi_{2}}{\xi_{1}}$ and $\tilde{w}(s):=\xi_{1}$, Gronwall’s inequality ensures that $\tilde{v}(t)\leq\tilde{\alpha}\cdot\mathrm{e}^{\int_{0}^{t}\tilde{w}(s)ds}$ and we infer the auxiliary statement. This, in turn, yields

[TABLE]

∎

Proposition 5**.**

Let $E\partial_{t}{x}=Ax+b$ be a regular linear DAE system and let $\mathcal{D}\subseteq\mathbb{R}^{n}$ denote the corresponding set of consistent initial conditions. Then, $\mathcal{D}$ is an affine subspace of $\mathbb{R}^{n}$ and $x+hF_{h}(x)\in\mathcal{D}$ whenever $x\in\mathcal{D}$.

Proof.

To see that $\mathcal{D}$ is an affine subspace of $\mathbb{R}^{n}$, please refer to [12, Section 2.1]. Note further that $x_{i}=x_{i-1}+hF_{h}(x_{i-1})$ defines the backward Euler scheme which is applied to the DAE system $E\partial_{t}{x}=Ax+b$, see [12, Section 5.2]. Consider the BDF-1 scheme [12, Section 5.3] which is given by

[TABLE]

if applied to $E\partial_{t}{x}=Ax+b$. With this, we first observe that

[TABLE]

where the inversion in the last line can always be performed for sufficiently small $h$ because $E\partial_{t}{x}=Ax+b$ is regular. This, in turn, yields

[TABLE]

This shows that the backward Euler scheme and the BDF-1 scheme are identical if applied to $E\partial_{t}{x}=Ax+b$. With this, the statement of the proposition is closely related to [12, Remark 5.25]. To see it, we may assume without loss of generality (see proof of [12, Theorem 5.24]) that $E\partial_{t}{x}=Ax+b$ is such that $A=I$ and $E=N$ for some nilpotent $N$ with $N^{\nu}=0$ and $N^{\nu-1}\neq 0$. It can be easily seen that in such a case the solution is $x\equiv-b$, thus implying in particular that the set of consistent initial conditions is $\mathcal{D}=\{-b\}$. Moreover, the BDF-1 scheme rewrites to

[TABLE]

where the last equivalence is due to the Neumann series and the nilpotency of $N$. This, in turn, implies that

[TABLE]

thus showing that $x_{i}=-b$ whenever $x_{i-1}=-b$. ∎

Proposition 6**.**

Let $E\partial_{t}{x}=Ax+b$ be a regular linear DAE system and let $\mathcal{D}\subseteq\mathbb{R}^{n}$ denote the corresponding set of consistent initial conditions. Then

• •

The solution of $E\partial_{t}{x}=Ax+b$ is contained in $\mathcal{D}$.

• •

There exist $\hat{A}\in\mathbb{R}^{n\times n}$ and $\hat{b}\in\mathbb{R}^{n}$ such that the solution of the ODE system $\partial_{t}{x}=\hat{A}x+\hat{b}$ coincides with that of $E\partial_{t}{x}=Ax+b$ for all $x(0)\in\mathcal{D}$.

• •

Together with $F_{h}(x):=(E-hA)^{-1}(Ax+b)$, where $h>0$, it holds that $F_{h}$ converges uniformly, as $h\to 0$, to $\hat{A}x+\hat{b}$ on any bounded subset of $\mathcal{D}$.

Proof.

The first two points are well-known in the theory of linear DAE systems, see Section [12, Section 2.1] (it is interesting to note that an efficient computation of $\hat{A}\in\mathbb{R}^{n\times n}$ and $\hat{b}\in\mathbb{R}^{n}$ is difficult because it relies on index reduction [11]).

To see third claim, we observe that $x_{i}=x_{i-1}+hF_{h}(x_{i-1})$ defines the backward Euler scheme applied to the DAE system $E\partial_{t}{x}=Ax+b$, see [12, Section 5.2]. We next show that $x_{0}\mapsto\frac{1}{h}(x_{1}-x_{0})$ converges uniformly on any bounded subset of $\mathcal{D}$ to $x_{0}\mapsto\hat{A}x_{0}+\hat{b}$ when $h\to 0$. To this end, we may assume without loss of generality (see discussion after Equation 5.25 in [12]) that the DAE system $E\partial_{t}{x}=Ax+b$ is such that

[TABLE]

where $N$ is such that $N^{\nu}=0$ and $N^{\nu-1}\neq 0$ for some $\nu\geq 1$. This implies that the solution of $E\partial_{t}{x}=Ax+b$ is characterized by a pair of decoupled dynamical systems, namely by the ODE system $\partial_{t}{x}^{\textup{I}}=Jx^{\textup{I}}+b^{\textup{I}}$ and the DAE system $N\partial_{t}{x}^{\textup{II}}=x^{\textup{II}}+b^{\textup{II}}$, where $x=(x^{\textup{I}},x^{\textup{II}})$ and $b=(b^{\textup{I}},b^{\textup{II}})$. Thanks to this, it suffices to consider $x^{\textup{I}}_{1}-x^{\textup{I}}_{0}$ and $x^{\textup{II}}_{1}-x^{\textup{II}}_{0}$ separately.

Since $x_{\textup{II}}\equiv-b_{\textup{II}}$ solves $N\partial_{t}{x}^{\textup{II}}=x^{\textup{II}}+b^{\textup{II}}$, we infer that $\mathcal{D}=\{(x^{\textup{I}},x^{\textup{II}})\mid x^{\textup{II}}=-b^{\textup{II}}\}$. Hence, Proposition 5 shows that $x^{\textup{II}}_{1}-x^{\textup{II}}_{0}=0$ whenever $x_{0}\in\mathcal{D}$.

We next focus on $x^{\textup{I}}_{1}-x^{\textup{I}}_{0}$. Thanks to the fact that $\partial_{t}{x}^{\textup{I}}=Jx^{\textup{I}}+b^{\textup{I}}$, we have to investigate the local truncation error of the backward Euler scheme in the context of a linear ODE system. Despite the fact this is discussed in many books about ODEs, we provide here a proof because most texts do not show that the local truncation error converges uniformly to zero on arbitrarily large compact sets. To this end, we first observe that the Taylor expansion of $x^{\textup{I}}$ around zero yields

[TABLE]

for some $\xi\in(0;h)$. With $\tilde{F}_{h}(x^{\textup{I}}_{0})=(I-hJ)^{-1}(Jx^{\textup{I}}_{0}+b^{\textup{I}})$, the proof of Proposition 5 implies that $\tilde{F}_{h}(x^{\textup{I}}_{0})=\frac{1}{h}(x^{\textup{I}}_{1}-x^{\textup{I}}_{0})$. This, in turn, implies that

[TABLE]

In the case $h\leq 1/(2\lVert J\rVert)$, the Neumann series allows us to deduce that

[TABLE]

with $\lVert\sum_{k=0}^{\infty}(hJ)^{k}\rVert\leq\sum_{k=0}^{\infty}2^{-k}=2$. Moreover, a differentiation of $\partial_{t}{x}^{\textup{I}}=Jx^{\textup{I}}+b^{\textup{I}}$ yields $\ddot{x}^{\textup{I}}=J^{2}x^{\textup{I}}+Jb^{\textup{I}}$. This and the last statement imply the existence of constants $\zeta_{1},\zeta_{2}\geq 0$ that neither depend on $x^{\textup{I}}_{0}$ nor on $h$ and that satisfy

[TABLE]

for all $0\leq h\leq 1$. This shows that $x^{\textup{I}}_{0}\mapsto\frac{1}{h}(x^{\textup{I}}_{1}-x^{\textup{I}}_{0})$ converges uniformly on any bounded set to $x^{\textup{I}}_{0}\mapsto Jx^{\textup{I}}_{0}+b^{\textup{I}}$. ∎

We are in a position to prove Theorem 1.

Proof of Theorem 1.

Let $\hat{A}\in\mathbb{R}^{n\times n}$ and $\hat{b}\in\mathbb{R}^{n}$ be as in Proposition 6 and fix $T>0$ and $x(0)\in\mathcal{D}$. Since the solution of $E\partial_{t}{x}=Ax+b$ solves the linear ODE system $\partial_{t}{x}=\hat{A}x+\hat{b}$, this implies that $x$ exists and is bounded on $[0;T]$. Hence, there exists a closed ball $B_{\rho}(0)$ centered at $0\in\mathbb{R}^{n}$ with radius $\rho>0$ such that $x(t)\in B_{\rho}(0)$ for all $0\leq t\leq T$. Since $B_{\rho}(0)$ is bounded, Proposition 6 ensures that $x$ is contained in $\mathcal{D}$ and that for any $\eta>0$ there exists an $h>0$ such that

[TABLE]

Moreover, Proposition 5 ensures that the solution $x_{h}$ of $\partial_{t}{x}_{h}=F(x_{h})$ is contained in $\mathcal{D}$. By combining the foregoing statements, Proposition 4 yields the claim. ∎

The following auxiliary results are needed for the proof of Theorem 2

Proposition 7**.**

Fix $E,A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times(k+m)}$, $D\in\mathbb{R}^{(k+m)\times(k+m)}$, $d\in\mathbb{R}^{(k+m)}$ and consider the linear DAE system

[TABLE]

Then, $(\hat{E}-h\hat{A})^{-1}(\hat{A}\hat{x}+\hat{b})$ is given by

[TABLE]

Proof.

By relying on the inversion formula for block matrices, we obtain

[TABLE]

Armed with this, we infer that

[TABLE]

and

[TABLE]

A summation of the foregoing statements yields (23). ∎

Corollary 1**.**

Fix an arbitrary consistent initial condition $(x(0),u(0))^{T}\in\mathbb{R}^{n+k+m}$ of the DAE system from Proposition 7. The corresponding ODE approximation is then

[TABLE]

with $u^{\langle 0\rangle}_{h}(0)=u(0)$ and $u^{\langle 1\rangle}_{h}(0)=(I-hD)^{-1}Du(0)$.

Proof.

Follows directly from Proposition 7. ∎

Proof of Theorem 2.

In Theorem 2, replace $u$ with $\hat{u}$, $z$ with $\hat{z}$ and $B$ with $\hat{B}$. Afterwards, apply Corollary 1 to the case where $u:=\binom{\hat{u}}{\hat{z}}\in\mathbb{R}^{k+m}$ and

[TABLE]

∎