Synthesizing Configurable Biochemical Implementation of Linear Systems from Their Transfer Function Specifications

Abstract

The ability to engineer synthetic systems in the biochemical context is constantly being improved and has a profound societal impact. Linear system design is one of the most pervasive methods applied in control tasks, and its biochemical realization has been proposed by Oishi and Klavins and advanced further in recent years. However, several technical issues remain unsolved. Specifically, the design process is not fully automated from specification at the transfer function level, systems once designed often lack dynamic adaptivity to environmental changes, matching rate constants of reactions is not always possible, and implementation may be approximative and greatly deviate from the specifications. Building upon the work of Oishi and Klavins, this paper overcomes these issues by introducing a design flow that transforms a transfer-function specification of a linear system into a set of chemical reactions, whose input-output response precisely conforms to the specification. This system is implementable using the DNA strand displacement technique. The underlying configurability is embedded into primitive components and template modules, and thus the entire system is adaptive. Simulation of DNA strand displacement implementation confirmed the feasibility and superiority of the proposed synthesis flow.

About the Authors:

Tai-Yin Chiu

Affiliations Department of Physics, National Taiwan University, Taipei, Taiwan, Graduate Institute of Electronics Engineering, National Taiwan University, Taipei, Taiwan

Hui-Ju K. Chiang

Affiliations Graduate Institute of Electronics Engineering, National Taiwan University, Taipei, Taiwan, EPI Lifeware, Inria Paris-Rocquencourt, Rocquencourt, France

Ruei-Yang Huang

Affiliation: Graduate Institute of Electronics Engineering, National Taiwan University, Taipei, Taiwan

Jie-Hong R. Jiang

* E-mail: [email protected]

Affiliations Graduate Institute of Electronics Engineering, National Taiwan University, Taipei, Taiwan, Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan

François Fages

Affiliation: EPI Lifeware, Inria Paris-Rocquencourt, Rocquencourt, France

Introduction

Advancements in synthetic biology have resulted in the development of biochemical systems of increasing complexity that are capable of using living cells as well as cell-free systems. Synthetic biology holds promise for biotechnology, biomedicine, bio-environmental, bioenergy, and other applications. Many computation and control design examples have been demonstrated either in vivo or in vitro. For instance, oscillators [1], toggle switches [2], logic gates [3], band pass filters [4], and analog circuits [5] have been designed and implemented in living cells, while digital circuits [6], neural networks [7], and switchable memories [8] have been demonstrated in cell-free systems. These developments show a clear growing trend in design complexity. Indeed, the more sophisticated systems we can build, the better skills we will have to comprehend biology and to build even more sophisticated systems.

The increasing system complexity renders the necessity of design automation tools. Constructing biochemical systems bottom-up using pre-characterized parts in synthetic biology is analogous to designing electronic systems using pre-designed standard modules. While electronic design automation (EDA) has largely enabled system design under exponential capacity growth over the past five decades thanks to the driving force of Moore’s law, design automation in synthetic biology may play a similar key role in the construction of complex biochemical systems [9]. Computer-aided modeling, simulation, synthesis, and verification are crucial because a biochemical design is often intended to function in a biological context that is often too complex to be fully characterized.

Depending on its target application, a system is often designed for a specific application rather than as a general purpose computing machine. In control applications, linear systems are pervasive due to their simplicity of design and analysis. Any linear control system can be realized with three primitive components: integration, gain, and summation. Realizing linear control with biochemical reactions has been proposed by Oishi and Klavins [10], who showed that three types of chemical reactions (catalysis, degradation, and annihilation) are sufficient to implement the three primitive components. In principle, any polynomial ordinary differential equation can be approximated by chemical reaction networks (CRNs) [11]. Therefore, any linear control system can be built using CRNs.

To realize CRNs, nucleic acids have been exploited as a universal tool for biomolecular computation [12]. By appropriately changing their nucleotide sequences, the interaction between nucleic acids can be precisely controlled and programmed in vitro. Specifically, toehold-mediated or enzyme-free DNA strand displacement (DSD) is a promising approach to perform biological computations. Although the DSD mechanism has been studied since the 1970s [13–17], until recently it was systematically used to build a molecular machine made of DNA and RNA [18]. The DSD method is attractive for several reasons. First, the kinetics of DSD devices can be engineered by controlling toehold binding rates, which may range from 1 to 6 × 10⁶ M⁻¹ s⁻¹ [19–21]. Second, composability [22] can be achieved. In typical implementation schemes, single-stranded DNAs play the role of signals and double-stranded DNAs act as gates, and they work together with the mechanism of toehold-mediated DNA displacement [19]. DNA signals of different domain lengths [23–25] have been proposed. For example, gates can be composed to realize complex systems using a 2-domain signal architecture [25]. Third, the DSD mechanism works autonomously [26] as long as the DNA or RNA fuels are supplied. With these advantages, many biomolecular devices, e.g., [27–29], have been designed using this powerful mechanism. The DSD technology allows computation and interfacing with molecular components in living organisms [30] and shows potential in bio-sensing and control, biomedicine, and other applications.

Despite the advancement of biochemical implementation of linear systems, four challenges remain to be solved. First, the construction of Oishi and Klavins’s system requires that the rate constants of the underlying reactions be carefully matched to achieve the intended integration, gain, and summation functions. This requirement imposes substantial practicality restrictions because, in reality, the reaction rates of available reactions can be limited. Thus, not all gain and summation blocks can be realized. Second, once a system is constructed, its function is fixed and cannot be changed without redesign. However, this fixed functionality can be inadequate for a system reacting to its biochemical environment, which is intrinsically stochastic and often full of uncertainty. Thus, designing systems with dynamic adaptation capabilities is crucial, especially in the biochemical context. Third, the CRN implementations of the gain and summation components are approximative. In essence, the transfer functions of the gain and summation components contain extra poles rather than just the ideal scalar and summation. When a complex system is built from these components, these additional poles might cause system behavior that deviates from its specification and can even lead to unwanted instability. Fourth, an automated flow of the synthesis of CRN from the transfer function specification of a linear system is lacking. Although a linear system can be compiled into a CRN with a direct block-by-block conversion from a block diagram, the block diagram may not be available in the first place and the CRN implementation may only approximate the specification. Specifying a linear system using transfer functions can be more natural than using block diagrams. However, converting a transfer function into a block diagram suitable for CRN realization can be a nontrivial optimization task.

In this paper, we tackle the above challenges. To address the first two challenges, we devise a mechanism to make linear control systems configurable by adding auxiliary species to CRN as control knobs. The concentrations of the auxiliary species can be adjusted not only to compensate for reaction rate mismatch but also to reconfigure a control system to dynamically adapt to environmental changes. Hence, implementing linear control systems in biochemistry can be made more practical. To resolve the last two challenges, we apply parallel decomposition on a given transfer function and express it as a summation of elementary modules, and propose a CRN solution that achieves the exact implementation of the elementary modules. The CRN implementation of a transfer function can be further mapped and realized using the DSD method. The proposed method lends itself to an automated design flow (as depicted in Fig 1) where the linear system to be synthesized is specified by a transfer function, which is decomposed for CRN implementation and further mapped into the DSD reactions. Through simulation using Visual DSD [31, 32], we show the feasibility and superiority of our proposed design automation flow for synthesizing a specified linear system into DSD reactions.

[Figure omitted. See PDF.]

Fig 1. Automated flow of linear system synthesis.

https://doi.org/10.1371/journal.pone.0137442.g001

Methods and Results

Although synthetic biology shows promising potential in biotechnology applications, it poses grand challenges to complex system design. A typical reactive system may involve sensing, actuation, information processing, and other tasks, while a control unit is one of the fundamental ingredients that constitute a reactive system. In this work, we address the design automation problem of linear control systems that are to be realized by biochemical reactions.

Any linear time-invariant (LTI) system can be fully specified by means of a transfer function, which corresponds to the Laplace transform of the impulse response of the system, with all initial conditions set to zero. In essence, once the impulse response of a linear system is known, the output y(t) of the system with any input u(t) can be characterized using the transfer function. However, for a multiple-input multiple-output (MIMO) linear system, the principle of superposition suggest that the total effect on any output due to all the inputs acting simultaneously is obtained by summing up the output effects due to the individual inputs acting alone. Consequently any MIMO system can be specified by a matrix of transfer functions. For simplicity, in this paper we focus on single-input single-output (SISO) systems without loss of generality.

Configurable Primitive Components

A transfer function can be realized by building a network of primitive components consisting of integration, gain, and summation blocks; therefore, implementing an LTI system in biochemical systems amounts to realizing these three basic components with chemical reactions [10]. However, the CRN realizations of these components proposed by Oishi and Klavins [10] are deficient in that rate constants have to be carefully matched and the designed control system cannot be reconfigured. These deficiencies impose a practicality issue that restricts system realization and a flexibility issue that prevents dynamic system adaptation. We can overcome these issues by introducing configurable primitive components as discussed below.

Following Oishi and Klavins [10], we represent a real signal x by the difference (x⁺ − x⁻) between the concentrations of two complementary molecular species x⁺ and x⁻. In this paper, we do not distinguish notationally between a species and its concentration. Moreover, we abbreviate a pair of chemical reactions of complementary specieswith rate constants r⁺ and r⁻, towhere the plus and minus signs in the superscripts ± and ∓ are ordered.

Table 1 shows the chemical reactions and transfer functions of the proposed configurable primitive components, which we will elaborate on in the following sections.

[Figure omitted. See PDF.]

Table 1. Three primitive components, their chemical reactions, and their transfer functions.

https://doi.org/10.1371/journal.pone.0137442.t001

Integration Block.

An integration block takes a signal u(t) as input, and produces a signal , for some constant k ∈ ℝ, as output. The chemical realization of an integration block for k ≥ 0 consists of a pair of catalytic reactionsand an annihilation reactionwhere are the rate constants. The reactions differ from those of Oishi and Klavins [10] in that the auxiliary species x^± are newly added to the catalytic reactions. Both the auxiliary species x^± and the input species u^± serve as catalysts.

With the definition that , the signal y is exactly the integration of signal u as described by the ordinary differential equationfor . Taking the Laplace transform of the above equation, we obtain the transfer function .

Because the concentrations of x⁺ and x⁻ can be controlled externally, in theory it is always possible to design a reaction network to meet any required k. For k < 0, the reactions are the same except that the catalytic reactions should be modified by reversing the complementary species of the output signal y aswhere y^∓ replaces the original y^±.

Gain and Summation Blocks.

A weighted summation block takes a number of input signals u_i(t), i = 1, 2, …, n and produces an output signal for k_i ∈ ℝ. A gain block is a special weighted summation block with only one input u(t) that produces output y(t) = k₁ u(t) for k₁ ∈ ℝ. The gain block with k₁ ≥ 0 can be realized by one pair of catalytic reactions,one pair of degradation reactions,and an annihilation reactionwhere x^± and z^± are auxiliary species. These reactions induce the following kinetic equationsLet and . The mass action of y becomes(1)When the steady state is reached, the changing rate of y (i.e., ) equals zero, which implies . Note that the concentrations of the auxiliary species x^± and z^± are controlled externally. Therefore, it becomes theoretically possible to meet any required k₁ by tuning the concentrations of these auxiliary species. For k₁ < 0, the plus and minus signs of the superscripts of y in the pair of catalytic reactions should be swapped.

For the weighted summation block, the reactions are the same as those for the gain block except that the pair of catalytic reactions becomesBy setting we obtain the equation(2)in the steady state. Similarly, if a scaling factor k_i < 0, we swap the signs in the superscript of y in the reaction corresponding to input u_i.

It is worth noting that in contrast to the integration block, the weighted summation block only approximates the intended weighted summation when the steady state is not yet reached; this is often the case in practice. This approximation can be clearly seen from its transfer function , rather than the ideal . Although this approximation may seem to inevitably cause a chemical reaction implementation to deviate from its specification, we will later show methods that can be used to circumvent this imperfection.

Case Study.

To assess the benefit of the auxiliary species, we perform a proof-of-concept case study on the mass-spring-damper (MSD) system as shown in Fig 2A. The system can be modeled by the equationUsing M = 1 kg, b = 10 N s/m, k = 20 N/m, and F = 1 N, by Laplace transform we derive the transfer functionThe transfer function can be implemented with the block diagram shown in Fig 2B. The proportional-integral (PI) controller to the MSD system is shown in Fig 2C.

[Figure omitted. See PDF.]

Fig 2. Mass-spring-damper system and its block diagrams.

(A) MSD system. (Let F = 1 N, M = 1 kg, b = 10 N s/m, and k = 20 N/m.) (B) Block diagram of the MSD model, where the triangular blocks denote gain functions with their corresponding weights, the rectangular blocks denote integrators, and the circle blocks denote mixers for summation and/or subtraction. (Let A = 2.764, B = 7.236 and C = 0.2236.) (C) Block diagram of the PI-controlled MSD model, where G is the plant shown in (B). (Assume the values of K_P and K_I given in Table 2.)

https://doi.org/10.1371/journal.pone.0137442.g002

The block diagrams are constructed with k₀ = 10 for all the summation and gain blocks, with the exception of the summation blocks shown in red and gain blocks A and B with k₀ = 50. The values of k_i’s are set to αk₀, where the values of α are equal to the weights specified in the corresponding gain blocks. Additionally, we assume and have the same values as k₀, and and are mismatched to k₁ by 10%, where q₁ and q₀ are the rate constants of the catalytic and degradation reactions, respectively, formulated in the prior work [10].

Fig 3A1, 3A2, 3A3, and 3B show the step, impulse, and sinusoidal responses of the MSD system, and the step responses of the PI-controlled MSD system, respectively. Our method achieves better approximation to the ideal cases than the prior method [10]. One of the advantages of our method is that we can match the weight by tuning the concentrations of x^± and z^±. In contrast, no tuning is possible in the prior method [10] to avoid the inexact gain , where q₁ and q₀ are the rate constants of the catalytic and degradation reactions formulated in the prior work [10], respectively, due to the mismatch of the rate constants. (Note that the biochemical implementations have their own optimal K_P and K_I values, shown in Table 2, to approximate the ideal system.)

[Figure omitted. See PDF.]

Fig 3. Responses of MSD systems with and without a PI controller.

The blue, green, and red curves represent the responses in ideal, configurable biochemical implementation, and nonconfigurable biochemical implementation cases, respectively. (A) Step, impulse, and sinusoidal responses of the MSD. (B) Step responses of PI-controlled MSD. (Assume 10% rate mismatch in the MSD system.) (C) Step responses of PI-controlled MSD, where the MSD undergoes a parameter change with b = 40 N s/m and k = 60 N/m, respectively, to induce a gain change of A = 1.561, B = 38.44 and C = 0.0271.

https://doi.org/10.1371/journal.pone.0137442.g003

[Figure omitted. See PDF.]

Table 2. Values of (K_P, K_I) in ideal, configurable and non-configurable implementations of the original and new MSD systems.

https://doi.org/10.1371/journal.pone.0137442.t002

Suppose that the spring and damper of the above MSD system are now replaced with new ones for b = 40 N s/m and k = 60 N/m. Our method can still adapt the PI-controller to the new MSD system without redesigning the PI-controller, whereas the prior method has no such capability. Because we can tune the concentrations of x^± and z^± in biochemical implementation, it is possible for us to adapt (K_P, K_I) to optimal values (40, 60) for the new PI-controlled MSD system, which is in contrast to the original system (15, 20). Fig 3C compares the results with and without such reconfigurability.

DSD Realization of Configurable Primitive Components

We exploit the DNA strand displacement (DSD) technique [24, 33] as an experimental chassis for our synthesis flow, and map the synthesized CRNs to DSD reactions. The simulation is conducted using the Visual DSD tool [34] for validation.

To implement each primitive component, we map the reactions listed in Table 1 to DSD reactions considering the compatibility among the components. We adopt the two-domain DSD method [25] to compile CRNs into nucleic acid-based chemistry to achieve flexible composability. In this method, any free (i.e., unbound) single-stranded DNA species consist of two domains: one toehold domain and one recognition domain. A toehold may initiate the binding between a single-stranded DNA species and a double-stranded DNA species. The two-domain DSD method supports compositionality, which indicates that a block implemented with DSD reactions can be cascaded with other blocks implemented with DSD reactions. Thus, the output species of one block fits the input species required by its downstream blocks.

As a notational convention, we indicate a domain t to be a toehold by a hat “” as . We also use superscript * to indicate the Watson-Crick complement of a domain DNA sequence (e.g., domain t* is sequence CTAG if domain t is sequence GATC). For a single-stranded species, we indicate the 5′ end using a vertical bar “∣” and the 3′ end by an angle bracket “⟩” or “⟨.” We call a single-stranded species of the form | ⟩ (5′ end on the left and 3′ end on the right) and ⟨ | (5′ end on the right and 3′ end on the left) an upper strand and a lower strand, respectively. A double-stranded species with an upper strand |t u⟩ (and a lower strand ⟨t* u*|) is denoted as [[t u]]. A double-stranded species with one or more dangling (or exposed) toeholds is called a gate, and is used to help transduce one or more species into other species. Double-stranded species without any dangling (exposed) toeholds or single-stranded species containing no toeholds are regarded as waste because they lack the ability to react with other species. The gates and species separately used to construct catalysis, degradation, and annihilation CRNs are called fuel.

To ignore the effects of fuel depletion, we maintain all of the separately used species at constant concentrations. Moreover, we set the binding rates of all toeholds to 0.05 nM⁻¹ s⁻¹, which is consistent with the prior work [35]. We assume infinite unbinding rates [32, 36] to ignore the unproductive reactions [32] that result from the situation where a single-stranded species binds to a gate but does not further displace the strand with which it competes for the same recognition domain. Under these settings, we map the catalytic chemical reactionsto the DSD reactions shown in Fig 4, where the domains highlighted in red and green are toeholds. Reactions and can be realized by the same set of DSD reactions except for species renaming (on recognition domains). For brevity, we simply use the same set of species (catalysis_7), (sp_0), and (sp10) to denote the species x⁺, u⁺, and y⁺, respectively, to implement the former reaction and to denote x⁻, u⁻, and y⁻, respectively, for the latter reaction.

[Figure omitted. See PDF.]

Fig 4. DSD realization of catalytic reactions.

The species in green are fuels whose concentrations are fixed to 2000 nM to avoid the effect of fuel depletion. The species sp_0, sp10, and catalysis_7 represent the input u⁺ or u⁻, output y⁺ or y⁻, and catalyst x⁺ or x⁻, respectively. Species catalysis_7 reacts with the gate catalysis to begin the first stage of the reactions followed by a series of reactions involving species sp_0 and catalysis_2. At the end of the first stage, species sp16 is produced and reacts with the gate catalysis_1 to begin the second stage of the reactions. In this stage, species catalysis_7 and sp_0 are yielded to compensate for their consumption in the first stage. Finally, the output species sp10 is generated.

https://doi.org/10.1371/journal.pone.0137442.g004

The set of DSD reactions outlined in Fig 4 consists of two stages. Gate catalysis initiates the first stage of the reactions by reacting with . At the end of the first stage, the species will be produced, and this species will react with gate catalysis_1 to initiate the second stage. During the second stage reactions, the species and will be produced to compensate for their consumption in the first stage.

We note that if the recognition domains of u and x were exchanged within the seven gates and one waste in the four first stage reactions, the entire set of reactions would still represent an implementation of the intended catalytic chemical reactions. However, the new arrangement of followed by would accelerate input (effectively, u⁺ or u⁻ in the catalytic chemical reactions) consumption because the input would directly react with the new catalysis gate. This is in contrast to the original catalysis structure, where is followed by . Therefore, this consumption acceleration of would amplify the bias between the production and consumption rates of . Because the concentration of could not remain constant, the input species could not be treated as a catalyst as expected. As a result, the new set of DSD reactions will not be as effective as the example presented in Fig 4.

To calibrate the values of the rate constants , we conduct an experiment on the catalytic reaction . Let x^± be mapped to , u^± to and y^± to at the concentrations summarized in Table 3. Additionally, we set and maintain the concentrations of all fuel to 2000 nM. However, the concentration of auxiliary species can be obtained by the weight after are determined. Based on the relationship , where are the concentrations of after an elapsed time τ = 1000 seconds, we derive the rate constants . As seen from Table 3, the reaction rates hold almost the same value when the concentrations of and are smaller than 100 nM, but slightly decrease when one of the concentrations approaches 100 nM. Because and will be low concentrations (approximately 10⁻⁶ to 1 nM) in our later experiments, we set nM⁻¹ s⁻¹. Fig 5A and 5B show the concentration of over time under the concentration (in nM) settings of at (0.001, 0.001) and (100, 100), respectively. These results suggest that the DSD reactions produce the output y^± (i.e., |t₂ y⟩) at a constant rate. The linearity of the DSD reactions indicates their adequacy in realizing our proposed configurable primitive components.

[Figure omitted. See PDF.]

Table 3. Experimental data used for measuring the rate constant of DSD catalysis.

https://doi.org/10.1371/journal.pone.0137442.t003

[Figure omitted. See PDF.]

Fig 5. Concentration of over time in DSD realized catalytic reactions.

(A) in nM and (B) in nM.

https://doi.org/10.1371/journal.pone.0137442.g005

Similar to the construction proposed by Yordanov et al. when mapping chemical reactions to DSD reactions [37], we map the degradation chemical reactionsto the DSD reactions shown in Fig 6. For brevity, the DSD reactions shown only correspond to one of the reactions and . The species , labeled deg, is the only gate or fuel used in this network, while the species represents y^±. The products sp3 and sp4 are waste products. In contrast to the DSD reactions proposed by Yordanov et al. [37] where the fuel concentrations are not maintained as constants, we maintain the fuel deg at a constant concentration. Therefore, the reaction can be effectively regarded aswhich is the same form as the above degradation chemical reactions. Thus, the species deg corresponds to z⁺ or z⁻, and the rate constants equal 0.05 nM⁻¹ s⁻¹, which is the binding rate constant of a toehold. By solving the differential equation of the kinetics of y^±, we derive that its concentration at time t is y(t)^± = y(0)^±exp[−(0.05z^±)t], where y(0)^± are the initial concentrations of y(t)^±. Fig 7A and 7B plot y(t)^± for z^± = 1 nM and 100 nM, respectively. In each case, there are three curves corresponding to different initial concentrations of y(0)^± = 10, 50 and 100 nM. These six curves are shown in Fig 7A and 7B and exactly fit the ideal waveforms.

[Figure omitted. See PDF.]

Fig 6. DSD realization of the degradation reaction.

Here, deg serves as fuel with its concentration maintained at a fixed value and corresponds to catalysts z^±. In the reaction, the rate constant equals the toehold binding rate 0.05 nM⁻¹ s⁻¹.

https://doi.org/10.1371/journal.pone.0137442.g006

[Figure omitted. See PDF.]

Fig 7. Simulation of DSD realized degradation reactions.

(A) Concentration of over time under deg = z^± = 1 nM. (B) Concentration of over time under deg = z^± = 100 nM. The three curves in each case correspond to different initial concentrations y(0)^± = 10, 50, and 100 nM.

https://doi.org/10.1371/journal.pone.0137442.g007

Finally, following the construction scheme of Yordanov et al. [37], we map the annihilation chemical reaction(3)to the DSD reactions provided in Fig 8, where species and represent y⁺ and y⁻, respectively. These species correspond to the output species of the primitive components. Two fuels with similar structures are used to make the entire DSD reaction symmetric and balance the consumption rates of and . As shown in Fig 8, the fuel ann first reacts with and then the product sp11 reacts with . In contrast, the fuel ann_1 first reacts with followed by the reaction between the product sp15 and . Therefore, the consumption rates of and are balanced and achieve effective annihilation.

[Figure omitted. See PDF.]

Fig 8. DSD realization of the annihilation reaction.

Here, and correspond to the output species y⁺ and y⁻, respectively. Two fuels (ann and ann_1) with similar structures are used to ensure that the consumption rates of y⁺ and y⁻ are balanced.

https://doi.org/10.1371/journal.pone.0137442.g008

The above DSD reactions for catalysis, degradation, and annihilation form a set of basic elements that enable the construction of primitive components and any proper linear system. We built the proposed configurable primitive components of integration and weighted summation with DSD reactions. The simulation results are shown in Fig 9. Fig 9A shows the output response of the integration component computing with respect to fixed inputs u⁺ = 2 nM and u⁻ = 1 nM. The concentration of y⁻ (green curve) remains at approximately 8 nM due to the effect of the annihilation reaction [Eq (3)], while the concentration of y⁺ (blue curve) grows to 105 nM after 100 seconds. As expected, the output signal y = y⁺−y⁻ grows linearly and reaches 97 nM after 100 seconds. Fig 9B shows the simulation results of the weighted summation component computing y = 2u₁ + u₂. Two sets of inputs and in nM are observed. As expected, the component functions correctly with y = 12 − 7 = 5 (for y⁺ = 2 × 5 + 2 = 12 and y⁻ = 2 × 3 + 1 = 7) nM.

[Figure omitted. See PDF.]

Fig 9. Simulation of DSD realized primitive components.

(A) Concentrations of y⁺, y⁻, and y = y⁺−y⁻ of the integration component computing under constant inputs u⁺ = 2 and u⁻ = 1 nM. (B) Concentrations of y⁺, y⁻, and y = y⁺−y⁻ of the weighted summation component computing y = 2u₁ + u₂ under constant inputs and in nM.

https://doi.org/10.1371/journal.pone.0137442.g009

The above simulation results confirm the successful realization of the integration and weighted summation components using DSD reactions. Because good compatibility is observed among the blocks constructed using the two-domain DSD method, these primitive DSD reactions can be composed to realize sophisticated linear systems. Later we will show that any proper transfer functions can be realized using the DSD realized primitive components.

Implementation of Transfer Functions in Biochemistry

Transfer Function Decomposition.

The transfer function G(s) of an LTI system can always be written as(4)where Y(s) and U(s) are the Laplace transforms of the output y(t) and input u(t) signals, respectively. A transfer function is called strictly proper, as we shall assume, if m < n (i.e., the degree of the numerator polynomial is smaller than the denominator polynomial).

Using the fundamental theorem of algebra, U(s) can be rewritten aswhere α_i, β_i and γ_i are real numbers satisfying for stable LTI systems, l and p are non-negative integers, and n_i and m_i are positive integers. By partial fraction expansion, Eq (4) can be factorized aswhere A_ij, B_ij, and C_ij are real coefficients. Because the terms of and can be realized by cascading j times the elementary modules of and , respectively, we only need to implement the set(5)of elementary modules, where D₁, D₂ and D₃ are real numbers. We refer to the three elementary modules as the “degree-(1, 0),” “degree-(2, 0),” and “degree-(2, 1)” modules, respectively, according to the degrees of their denominator and numerator polynomials.

Below we investigate how to implement these three elementary modules to construct any LTI system with CRNs. First, we provide a naive construction to illustrate inexactness, and then present a refined exact solution.

Naive Implementation of Elementary Modules.

To build a degree-(1, 0) module J(j₁, j₂) with parameters j₁, j₂ with real values, one might use the negative feedback loop consisting of a forward integration block of and a negative feedback of gain j₂ (as shown in Fig 10A). The transfer function of this realization equals , which exactly represents the target form. Recall that the introduction of auxiliary species enables the configurability of primitive components, and thus the values of j₁ and j₂ can be set as desired.

[Figure omitted. See PDF.]

Fig 10. Block diagrams of naive implementation of elementary modules.

(A) Block diagram of the degree-(1, 0) module with parametric weights j₁ and j₂ to match the transfer function coefficients. (B) Block diagram of the degree-(2, 0) module with parametric weights e₁, e₂, and e₃ to match the transfer function coefficients. (C) Block diagram of the degree-(2, 1) module with parametric weights f₁, f₂, and f₃ to match the transfer function coefficients.

https://doi.org/10.1371/journal.pone.0137442.g010

To build a degree-(2, 0) module, one has to increase the degree of the denominator polynomial. Therefore, one might try to use a degree-(1, 0) module for a compositional construction (as shown in Fig 10B). Let the degree-(1, 0) module be J(e₁, e₂) fed to a forward integration block of , and let the negative feedback have gain e₃. Then, the transfer function of the block diagram equals as desired.

To build a degree-(2, 1) module, one might also try to use a degree-(1, 0) module for compositional construction (as shown in Fig 10C). Let the degree-(1, 0) module be J(f₁, f₂) fed forward to the output, and place an integration block of on the negative feedback path. Then, the transfer function of the block diagram equals as desired.

Although the block diagrams of Fig 10 yield the desired transfer functions of the elementary modules, they are problematic due to the imperfection of the weighted summation components realized by the chemical reactions (with transfer function instead of as discussed earlier in the primitive component construction section). The following example demonstrates that a system constructed in this fashion may deviate drastically from the specification and can even be unstable.

Consider the LTI system specified by transfer function(6)By partial fraction expansion, G can be expressed as(7)By defining(8)they correspond to the basic elements in Eq (5) and can be implemented with the block diagrams of Fig 10 by setting (j₁, j₂) = (2, 2.5), (e₁, e₂, e₃) = (1, 3, 3) and (f₁, f₂, f₃) = (1, 2, 2). Assume that k₀ = a for some positive real constant a for all weighted summation primitive components; then, assume that k_i, i = 1, …, n, of a weighted summation primitive component are determined according to the required weights. The transfer functions after the chemical reaction realization become(9)Based on (asymptotic) stability analysis, the implemented system is unstable when a pole of the transfer functions (i.e., a root of the denominator polynomials) falls in the right-half of the s-plane. To avoid instability, the value of a has to be carefully determined. However, choosing a proper value for a can be tedious and sometimes even impossible.

The responses of , , and to the step input with an amplitude of 10⁻⁶ are shown in Fig 11. Specifically, Fig 11A and 11B show the responses of , C and D show the responses of , and E and F show the responses of . As seen from the plots of Fig 11A, 11C, and 11E, the responses are unstable under the chosen small a values, while for plots of B, D, and F, the responses are stable under the chosen large a values. However, even the stable responses may deviate from the ideal transfer function responses to some extent. When a in , and increases to 6, the maximal concentrations of the auxiliary species used in each case are 6 × 3/0.0249 ≈ 720, 6 × 2/0.0249 ≈ 480, and 6 × 2.5/0.0249 ≈ 600 nM, respectively, where the factors 3, 2 and 2.5 are the weights of the gain blocks in each system. As seen in Fig 11, the naive implementations cannot approximate the ideal responses well even under such high concentrations.

[Figure omitted. See PDF.]

Fig 11. Responses of naive implementations to a step input of amplitude 10⁻⁶ nM under different values of parameter a.

(A) Responses of under a = 0.5, 1.0, 1.5. (B) Responses of under a = 4.0, 5.0, 6.0. (C) Responses of under a = 0.5, 1.0, 1.5. (D) Responses of under a = 4.0, 5.0, 6.0. (E) Responses of under a = 0.5, 1.0, 1.5. (F) Responses of under a = 4.0, 5.0, 6.0.

https://doi.org/10.1371/journal.pone.0137442.g011

Exact Implementation of Elementary Modules.

To prevent the non-ideal effect of the weighted summation components, we present a refined method that exactly implements the three elementary modules with chemical reactions.

The degree-(1, 0) module can be realized directly by the CRN transfer function of the primitive gain block. For the other two modules, we take advantage of the non-ideal form of the weighted summation block and devise their corresponding block diagrams, whose CRN transfer functions match the specification. The block diagrams, normal transfer functions, and CRN transfer functions for the three elementary modules are summarized in Fig 12. The CRN transfer function of the degree-(2, 0) module can be calculated from(10)Therefore, the CRN transfer function equalsNote that if the gain and summation blocks are ideal, the normal transfer function should be . Similarly, the CRN transfer function of the degree-(2, 1) module can be derived asAs a result, the responses can be made exact through the control of auxiliary species concentrations.

[Figure omitted. See PDF.]

Fig 12. Block diagrams and their transfer functions for the exact implementation of the elementary modules.

The first, second, and third rows in the table correspond to the implementations of the degree-(1, 0), degree-(2, 0), and degree-(2, 1) modules, respectively. The first, second, and third columns in the table show the block diagrams, normal transfer functions, and CRN transfer functions, respectively. A normal transfer function is obtained directly from its block diagram, whereas a CRN transfer function is obtained when the gain and summation blocks in the block diagram are realized using the gain and summation CRNs in Table 1.

https://doi.org/10.1371/journal.pone.0137442.g012

There are two special cases that need attention. First, if the coefficient β vanishes in the denominators of the transfer functions and of the degree-(2, 0) and degree-(2, 1) modules, respectively, we have to make k₀ = 0 in the CRN transfer functions. Therefore, we need to remove the degradation reactions from the summation blocks. Second, for the negative parameter k₁ or k₂ in Fig 12, as mentioned previously the plus and minus signs of the superscripts of y in the pair of catalytic reactions of the weighted summation CRN should be swapped; for negative k₃, the superscript signs of y in the pair of catalytic reactions of the integration CRN should be swapped. Notice that k₀ can always be non-negative for stable linear systems, whose poles should all be on the left-half of the s-plane. Therefore, with these modifications our method is sufficient to implement any stable linear system.

DSD Realization of Transfer Functions

We map the naive and exact implementations of the transfer function specification of Eq (6) into DSD reactions for comparison. In the simulation, we assume the rate constants nM⁻¹ s⁻¹ and nM⁻¹ s⁻¹ for i = 1, …, n.

The responses of DSD that realized naive implementation of the transfer function specification are shown in Fig 13. Specifically, Fig 13A and 13B show the results realizing the transfer function G₁ in Eq (8), C and D show those realizing G₂, and E and F show those realizing G₃. In the plots of Fig 13A, 13C, and 13E, the realized systems are unstable under the chosen small a values, while in the plots of B, D, and F the realized systems are stable under the chosen large a values. These results are consistent with the CRN responses of Eq (9) in Fig 11.

[Figure omitted. See PDF.]

Fig 13. Responses of DSD realized naive implementations to a step input of amplitude 10⁻⁶ nM under different values of parameter a.

(A) Responses of under a = 0.5, 1.0, 1.5. (B) Responses of under a = 4.0, 5.0, 6.0. (C) Responses of under a = 0.5, 1.0, 1.5. (D) Responses of under a = 4.0, 5.0, 6.0. (E) Responses of under a = 0.5, 1.0, 1.5. (F) Responses of under a = 4.0, 5.0, 6.0.

https://doi.org/10.1371/journal.pone.0137442.g013

The maximum k_i’s in the DSD realizations of , , and are 6.0 × 3 = 18.0, 6.0 × 2 = 12.0, and 6.0 × 2.5 = 15.0, respectively. Because the rate constants , i = 1, …, n are set to 0.0249 nM⁻¹ s⁻¹, the corresponding concentrations of the catalysts x^± in the gain and summation blocks in Fig 10 of , , and are approximately 18.0/0.0249 ≈ 720, 12.0/0.0249 ≈ 480, and 15/0.0249 ≈ 600 nM, respectively. Not only are such high concentrations undesirable, but it is also difficult to find a proper a value to fit the ideal response.

To compare the naive and exact implementations of transfer function specifications under DSD realization, the responses of the targeted transfer functions G₁, G₂ and G₃ are plotted in Fig 14A, 14B, and 14C, respectively. The red curves in Fig 14A, 14B, and 14C correspond to the naive implementations with a = 6; the involved concentrations of catalysts are approximately 720, 480, and 600 nM, respectively. The species with the largest concentrations used in the exact implementations of G₁ and G₂ are both the catalysts x^± in the integration CRNs. Because the weights of the integration blocks in the exact implementations of G₁ and G₂ are both equal to 2.0, the corresponding concentrations of x^± in the catalytic reactions are approximately 80 nM. In contrast, in the exact implementation of G₃ the catalysts z^± in the gain block (where k₀ equals 5.0) have the highest concentration (5.0/0.05 = 100 nM). Compared to the naive implementations involving concentrations as large as 720 nM, the exact implementations require much lower concentrations, and yet fit the ideal responses much better.

[Figure omitted. See PDF.]

Fig 14. Response comparisons on DSD realized systems under the naive and exact implementations of transfer functions G₁, G₂, and G₃.

(A) Responses of G₁ implementations. (B) Responses of G₂ implementations. (C) Responses of G₃ implementations. The curves in red corresponds to naive implementations (under the parameter setting a = 6 in ), while the curves in green correspond to exact implementations.

https://doi.org/10.1371/journal.pone.0137442.g014

For the overall transfer function G = G₁ + G₂ + G₃, the responses of the DSD reactions realized via the naive and exact implementations are compared in Fig 15 against the ideal response. There are two approaches combining the subsystems G₁, G₂, and G₃. One approach is to use a summation block to sum up the outputs Y₁, Y₂, and Y₃ of G₁, G₂, and G₃, respectively; the other approach is to keep Y₁, Y₂, and Y₃ intact without summation by assuming that the output species of Y₁, Y₂, and Y₃ share the same piece of the DNA segment (rather than the entire DNA segment) that corresponds to the intended final product. The result of the former approach is shown in Fig 15A and the latter approach is shown in Fig 15B, which corresponds to the direct superposition of individual G₁, G₂, and G₃ responses. We observe that summing Y₁, Y₂, and Y₃ with an additional summation block may result in slight distortions in comparison with direct superposition. In either case, the exact implementation fits the ideal response much better than the naive implementation.

[Figure omitted. See PDF.]

Fig 15. Response comparisons on DSD realized systems for the naive and exact implementations of transfer function G.

(A) G implemented by a summation block adding G₁, G₂, and G₃, assuming a = 6 for this summation block in both the naive and exact methods. (B) G implemented by superposing G₁, G₂, and G₃, assuming that the DSD output species of G₁, G₂, and G₃ have a common sub-domain that is recognized as the output of G.

https://doi.org/10.1371/journal.pone.0137442.g015

Discussion

Rate Matching and Configurability

The CRN of the integration block proposed by Oishi and Klavins contains the catalytic reactions . The integration is made possible under the assumption that . Accordingly (i.e., output y is a weighted integration of input u). However, the CRNs of the gain block contain one pair of additional degradation reactions . The gain function with weight w_gn is realized by assuming . Thus, not only but also must be satisfied.

To cope with the problem of matching and in the integration CRN, Oishi and Klavins provided a solution by maintaining a large rate constant q_ua in the annihilation reaction so that and do not need to be exactly the same, as long as they are close enough to w_int. The reason is that one of the concentrations of u⁺ and u⁻ would approach zero under fast annihilation, and thus equals either w_int u⁺ or w_int u⁻. Because and are close to w_int, the integration can be reasonably approximated. Similarly, the matching problem in the gain CRN is overcome by asserting not only a large q_ua, but also a large q_ya in . Thus, by requiring that are close enough and the ratios and are close to w_gn, the gain function can be well approximated.

However, the above assumptions of fast annihilations (and ) may not be easily satisfied, especially when these reactions are achieved by composing multiple reactions (e.g., in the DSD realization). Two undesirable situations might occur in composite reactions even if the rate constant q_ua is high. One is that u⁺ and u⁻ both degrade, but the species with the lower concentration converges to a non-zero concentration value due to the restoration of the reactants by other reactions. The other is that although the species with the lower concentration vanishes, the remaining species converges to a concentration value less than ∣u⁺−u⁻∣ because some fuels can react with the remaining species. The approximations may be unsound when these situations occur. Although the effect of inappropriate assumptions might not be serious for a single primitive component, for a complex system composed of many primitive components the approximation would be crude due to error accumulations.

In contrast, because our primitive components are configurable due to the addition of the auxiliary species, we can match not only the reaction rates but also the required weights by tuning the concentrations of auxiliary species. Another advantage provided by the use of auxiliary species is the ability for dynamic system adaptation. As discussed in the case study of the PI-controller, the parameters of a PI-controller may be modified by simply tuning the concentrations of auxiliary species to react to environmental changes. Thus, by using auxiliary species we may fulfill a new system with the same set of reactions by tuning only the concentrations of the auxiliary species. In contrast, without auxiliary species the system would have to be redesigned by finding new catalysts and degradations to match the new weights and rate constants.

Transfer Function Decomposition

A transfer function can be decomposed in many different ways. Two common choices are through parallel decomposition (G = G₁ + G₂) and serial decomposition (G = G₁ ⋅ G₂). In this paper, we adopt the former for the following reasons. First, given a strictly proper transfer function it is not always possible to decompose it as a product of the strictly proper transfer functions. For instance, the transfer function cannot be decomposed serially into a system involving the cascade of two strictly proper subsystems. Because only proper systems where the degree of the numerator polynomial does not exceed that of the denominator polynomial are physically realizable, the only hope for implementing the system physically is to have a product of bi-proper subsystems, where the degrees of the numerator and denominator polynomials are the same. For example, it can be decomposed as . However, if we implement the system with CRNs or the primitives introduced previously, the bi-proper subsystem might be an approximation due to , whose realization might involve a gain and a summation block, both of which are approximations. Second, by cascading subsystems using serial decomposition a system may be susceptible to error accumulation due to realization approximation. In contrast, parallel decomposition avoids the issue of error accumulation and amplification. Third, a strictly proper transfer function can always be decomposed as a summation of strictly proper transfer functions, which can be modeled by the three elementary modules described in Eq (5). Furthermore, we show that the three elementary modules allow exact CRN implementation using our configurable primitive components. For these reasons, parallel decomposition may be preferable to serial decomposition.

Nevertheless, applying parallel decomposition exclusively may not always be the best choice. The reasons are twofold. First, parallel decomposition requires one final summation block to sum up all elementary modules. This summation block, when it is implemented with a CRN, is approximative. Second, a hybrid decomposition strategy may achieve more effective implementation than applying a parallel decomposition. For instance, consider the system of transfer function . If only parallel decomposition is applied, then T is decomposed into . Nevertheless, if serial decomposition is also exploited, then can be rewritten as , which can be implemented by cascading two gain blocks with (k₀, k₁) = (1, 1) and (2, 1), respectively. Therefore, with this rewriting a final summation block only needs to sum up the outputs of the two subsystems rather than three subsystems. Furthermore, the hybrid parallel and serial decomposition strategy may be more cost effective than using parallel decomposition alone. Future work should exploit a good hybrid decomposition strategy and remove the approximative effect due to the final summation block in parallel decomposition.

DSD Reaction Rates

In this paper, we realize all CRNs with DSD reactions. Typically, the rate constants in DSD reactions are approximately 10⁻³ nM⁻¹ s⁻¹ [33]. For example, if one requires the parameter in the degradation reaction to be 100 s⁻¹, then the concentrations of z^± should be of magnitude 10⁵ in nM. However, such a high concentration might be impractical. To alleviate this high concentration requirement, one may try to increase the DSD reaction rates. Zhang et al. constructed and characterized DNA catalytic circuits driven by entropic gains [28]. Based on the entropy effects, a variant called the tethered entropy driven catalytic circuits was introduced [35, 38] to shorten the catalytic cycle and improve the reaction kinetics using localized hybridization reactions, which are achieved by tethering the key species in DSD reactions to increase their local concentrations. The influence of tethering is effectively presented in a speed-up factor λ, whose value can be feasibly reach approximately 10⁵ [35]. Thus, the toehold binding rate constants can be scaled up to λ × 5 × 10⁻⁵ nM⁻¹ s⁻¹, where 5 × 10⁻⁵ nM⁻¹ s⁻¹ is the short toehold binding rate constant in diffusion-based systems [6].

In our DSD simulation, we assume the factor λ = 1000 is available to increase the toehold binding rate constants to approximately 0.05 nM⁻¹ s⁻¹. Effectively, the concentrations of the auxiliary species can be reduced by a factor of 25. Under the toehold binding rate increase, a catalysis has a rate constant approximately 0.025 nM⁻¹ s⁻¹ in our interested range of concentrations under the intended system operation. Thus, the kinetics of the products y^± will now become . If x^± could be tuned to 1000 nM, then we can have an integration block with weight as high as 25.

Although we have reduced the required concentrations of auxiliary species using the speed-up factor λ, there are still some limitations. First, if β is set to 500 in the elementary modules and shown in Fig 12, then we have to prepare the catalysts z^± with large concentrations (approximately 500/0.05 = 10000 nM due to the relation ). In contrast, if , then the concentrations of x^± would be 500/0.025 = 20000 nM. Thus, although we have reduced the concentrations of the auxiliary species to increase the range of coefficients in the transfer functions, further reduction may be needed.

Fuel Supply

DSD is a relatively mature technology for reliable in vitro construction of biochemical systems, and is more practical than in vivo experiments. However, one of the main obstacles for practical DSD experiment is to provide constant supply of fuels [37]. This paper assumes the fuels are large in quantity and ignore the situations where the fuels become deficient. In reality, the fuels are converted to inert waste and one has to continually provide the fuels to the system. Providing stable fuel supply can be challenging in a cellular context. Therefore, it is easier to implement DSD reactions in vitro instead of in vivo.

Author Contributions

Conceived and designed the experiments: TYC RYH JHJ. Performed the experiments: TYC. Analyzed the data: TYC HJC RYH JHJ FF. Contributed reagents/materials/analysis tools: TYC HJC RYH JHJ FF. Wrote the paper: TYC HJC RYH JHJ FF.

Citation: Chiu T-Y, Chiang H-JK, Huang R-Y, Jiang J-HR, Fages F (2015) Synthesizing Configurable Biochemical Implementation of Linear Systems from Their Transfer Function Specifications. PLoS ONE 10(9): e0137442. https://doi.org/10.1371/journal.pone.0137442

References

1. Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000;403:335–338. pmid:10659856

2. Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in Escherichia coli. Nature. 2000;403:339–342. pmid:10659857

3. Tamsir A, Tabor JJ, Voigt CA. Robust multicellular computing using genetically encoded NOR gates and chemical ‘wires’. Nature. 2011;469:212–215. pmid:21150903

4. Sohka T, Heins RA, Phelan RM, Greisler JM, Townsend CA, Ostermeier M. An externally tunable bacterial band-pass filter. Proc Natl Acad Sci USA. 2009;106:10135–10140. pmid:19502423

5. Daniel R, Rubens JR, Sarpeshkar R, Lu TK. Synthetic analog computation in living cells. Nature. 2013;497:619–623. pmid:23676681

6. Qian L, Winfree E. Scaling up digital circuit computation with DNA strand displacement cascades. Science. 2011;332:1196–1201. pmid:21636773

7. Qian L, Winfree E, Bruck J. Neural network computation with DNA strand displacement cascades. Nature. 2011;475:368–372. pmid:21776082

8. Padirac A, Fujii T, Rondelez Y. Bottom-up construction of in vitro switchable memories. Proc Natl Acad Sci USA. 2012;109:E3212–E3220. pmid:23112180

9. Densmore D, Hassoun S. Design automation for synthetic biological systems. IEEE Design & Test of Computers. 2012;29:7–20.

10. Oishi K, Klavins E. Biomolecular implementation of linear I/O systems. IET Syst Biol. 2011;5:252–260. pmid:21823756

11. Klonowski W. Simplifying principles for chemical and enzyme reaction kinetics. Biophys Chem. 1983;18:73–87. pmid:6626688

12. Padirac A, Fujii T, Rondelez Y. Nucleic acids for the rational design of reaction circuits. Curr Opin Biotechnol. 2013;24:575–580. pmid:23265857

13. Holloman WK, Wiegand R, Hoessli C, Radding CM. Uptake of homologous single-stranded fragments by superhelical DNA: A possible mechanism for initiation of genetic recombination. Proc Natl Acad Sci USA. 1975;72:2394–2398. pmid:1094467

14. Green C, Tibbetts C. Reassociation rate limited displacement of DNA strands by branch migration. Nucleic Acids Res. 1981;9:1905–1918. pmid:6264399

15. Robinson BH, Seeman NC. Simulation of double-stranded branch point migration. Biophys J. 1987;51:611–626. pmid:3580487

16. Weinstock PH, Wetmur JG. Branch capture reactions: Effect of recipient structure. Nucleic Acids Res. 1990;18:4207–4213. pmid:2198537

17. Panyutin IG, Hsieh P. The kinetics of spontaneous DNA branch migration. Proc Natl Acad Sci USA. 1994;91:2021–2025. pmid:8134343

18. Yurke B, Turberfield AJ, Mills AP, Simmel FC, Neumann JL. A DNA-fuelled molecular machine made of DNA. Nature. 2000;406:605–608. pmid:10949296

19. Yurke B, Mills AP. Using DNA to power nanostructures. Genetic Programming and Evolvable Machines. 2003;4:111–122.

20. Li Q, Luan G, Guo Q, Liang J. A new class of homogeneous nucleic acid probe based on specific displacement hybridization. Nucl Acids Res. 2002;30:e5. pmid:11788731

21. Zhang DY, Winfree E. Control of DNA strand displacement kinetics using toehold exchange. J Am Chem Soc. 2009;131:17303–17314. pmid:19894722

22. Gyorgy A, Del Vecchio D. Modular composition of gene transcription networks. PLoS Computational Biology. 2014;10:e1003486. pmid:24626132

23. Cardelli L. Strand algebras for DNA computing. In: Deaton R, Suyama A, editors. DNA Computing and Molecular Programming. vol. 5877 of Lecture Notes in Computer Science. Springer Berlin Heidelberg; 2009. p. 12–24. Available from: http://dx.doi.org/10.1007/978-3-642-10604-0_2.

24. Soloveichik D, Seelig G, Winfree E. DNA as a universal substrate for chemical kinetics. Proc Natl Acad Sci USA. 2010;107:5393–5398. pmid:20203007

25. Cardelli L. Two-domain DNA strand displacement. Mathematical Structures in Computer Science. 2013;23:247–271.

26. Green SJ, Lubrich D, Turberfield AJ. DNA hairpins: Fuel for autonomous DNA devices. Biophysical Journal. 2006;91:2966–2975. pmid:16861269

27. Seelig G, Soloveichik D, Zhang DY, Winfree E. Enzyme-free nucleic acid logic circuits. Science. 2006;314:1585–1588. pmid:17158324

28. Zhang DY, Turberfield AJ, Yurke B, Winfree E. Engineering entropy-driven reactions and networks catalyzed by DNA. Science. 2007;318:1121–1125. pmid:18006742

29. Bath J, Turberfield AJ. DNA nanomachines. Nat Nanotechnol. 2007;2:275–284. pmid:18654284

30. Amir Y, Ben-Ishay E, Levner D, Ittah S, Abu-Horowitz A, Bachelet I. Universal computing by DNA origami robots in a living animal. Nat Nanotechnol. 2014;9:353–357. pmid:24705510

31. Phillips A, Cardelli L. A programming language for composable DNA circuits. J R Soc Interface. 2009;6:419–436.

32. Lakin MR, Youssef S, Cardelli L, Phillips A. Abstractions for DNA circuit design. J R Soc Interface. 2012;9:470–486. pmid:21775321

33. Zhang DY, Seelig G. Dynamic DNA nanotechnology using strand-displacement reactions. Nat Chem. 2011;3:103–113. pmid:21258382

34. Lakin MR, Youssef S, Polo F, Emmott S, Phillips A. Visual DSD: A design and analysis tool for DNA strand displacement systems. Bioinformatics. 2011;27:3211–3213. pmid:21984756

35. Chandran H, Gopalkrishnan N, Phillips A, Reif J. Localized Hybridization Circuits. In: Cardelli L, Shih W, editors. DNA Computing and Molecular Programming. vol. 6937 of Lecture Notes in Computer Science. Springer Berlin Heidelberg; 2011. p. 64–83. Available from: http://dx.doi.org/10.1007/978-3-642-23638-9_8.

36. Lakin MR, Parker D, Cardelli L, Kwiatkowska M, Phillips A. Design and analysis of DNA strand displacement devices using probabilistic model checking. J R Soc Interface. 2012;9:1470–1485. pmid:22219398

37. Yordanov B, Kim J, Petersen RL, Shudy A, Kulkarni VV, Phillips A. Computational design of nucleic acid feedback control circuits. ACS Synth Biol. 2014;3:600–616. pmid:25061797

38. Gopalkrishnan N, Garg S, Chandran H, Reif J. Speeding up DNA circuits using localized hybridization; 2011. In Foundations of Nanoscience Meeting (FNANO).