Supervised Learning in Adaptive DNA Strand ... - ACS Publications

Apr 25, 2016 - Supervised Learning in Adaptive DNA Strand Displacement Networks. Matthew R. Lakin*,†,‡,§ and Darko Stefanovic. ‡,§. †. Depar...
3 downloads 7 Views 2MB Size
Research Article pubs.acs.org/synthbio

Supervised Learning in Adaptive DNA Strand Displacement Networks Matthew R. Lakin*,†,‡,§ and Darko Stefanovic‡,§ †

Department of Chemical & Biological Engineering, ‡Department of Computer Science, and §Center for Biomedical Engineering, University of New Mexico, Albuquerque, New Mexico 87131, United States S Supporting Information *

ABSTRACT: The development of engineered biochemical circuits that exhibit adaptive behavior is a key goal of synthetic biology and molecular computing. Such circuits could be used for long-term monitoring and control of biochemical systems, for instance, to prevent disease or to enable the development of artificial life. In this article, we present a framework for developing adaptive molecular circuits using buffered DNA strand displacement networks, which extend existing DNA strand displacement circuit architectures to enable straightforward storage and modification of behavioral parameters. As a proof of concept, we use this framework to design and simulate a DNA circuit for supervised learning of a class of linear functions by stochastic gradient descent. This work highlights the potential of buffered DNA strand displacement as a powerful circuit architecture for implementing adaptive molecular systems. KEYWORDS: molecular computing, DNA strand displacement, machine learning, gradient descent, adaptive algorithms

N

incumbent strand bound to the template. The incumbent strand is complementary to the template but shorter, so there is a short single-stranded overhang on one end of the template strand, referred to as a toehold. The toehold is of critical importance because it provides a nucleation point for the interaction between the gate and the invader by holding the invader in proximity to the gate and in the correct orientation for the subsequent reaction to proceed.8 Toeholds are typically 5−8 nucleotides in length, so the binding of the invader to the toehold is a reversible process. (Strand displacement reaction rates are particularly sensitive to toehold length.9,10) Once the invader is bound to the toehold, the sequence commonality between the invader and incumbent strands allows the invader to compete to bind to the template strand. This happens due to the breathing of the base pair at the end of the incumbent−template duplex, which gives the invader strand a chance to replace the briefly opened bond with an invader−template bond. This is known as branch migration because the boundary between the incumbent− template duplex and the invader−template duplex migrates back and forth in a random walk process. When the branch migration reaches the far end of the template strand and the final base pair of the incumbent−template duplex is replaced by an invader− template base pair, the strand displacement reaction is complete. As the resulting invader−template duplex contains no singlestranded toehold, the displaced incumbent strand cannot rebind and, hence, the strand displacement reaction is typically modeled as an irreversible process. Furthermore, there is a thermodynamic

ucleic acid-based molecular computing systems are tools for information processing on the nanoscale. In these systems, computations on molecular input signals are performed by engineered reactions in a dynamic nucleic acid circuit. Such circuits can interface directly with biological and chemical systems at the molecular level, offering promising applications such as the development of autonomous “smart” therapeutics that could sense, diagnose, and treat disease within individual cells without external intervention. This article proposes a molecular circuit design framework capable of adaptive behavior, such as implementing machine learning algorithms, which could have practical applications in long-term autonomous monitoring and control of biological or chemical systems. DNA strand displacement has been proposed as an approach to molecular computing that enables sophisticated computations to be carried out by circuits that are built in a uniform, principled manner by composing computational units.1 DNA strand displacement reactions proceed by nucleic acid hybridization according to the rules of Watson−Crick complementarity, which makes it relatively straightforward to predict their behavior using thermodynamic and kinetic models.2,3 It has been shown that DNA strand displacement reaction networks are universal in the sense that they can faithfully emulate the dynamics of any abstract chemical reaction network.4 Hence, strand displacement networks provide a platform for realizing distributed algorithms, represented as abstract chemical reaction networks, in wet chemistry via a formal compilation process. A number of different schemes to achieve this encoding have been published in recent years.4−7 A basic strand displacement reaction involves a single-stranded invader and a multistranded gate complex that consists of a template strand that is complementary to the invader and an © XXXX American Chemical Society

Special Issue: DNA 21 Received: January 12, 2016

A

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

synthetic gene networks are capable of associative learning.20,21 It has even been shown that Arabidopsis plants carry out simple arithmetic division operations to enable adaptive management of their photosynthetic starch reserves in the face of variable reserves and variable hours of darkness.22 Developing engineered molecular circuits that exhibit adaptive behaviors could enable molecular computing solutions for practical applications that rely on processing inputs over extended periods of time, such as longterm health monitoring and adaptive support of failing cellular regulatory networks, detection of emerging pathogens23 whose signatures mutate over time, and in situ measurement of cellular processes over extended periods of time to answer basic science questions. Adaptive molecular computing networks could also serve as control systems for artificial life systems engineered using biomolecules. Hence, there are both practical and basic science motivations for the development of adaptive molecular computing systems. In this article, we present a DNA circuit architecture that is generally applicable to the implementation of adaptive algorithms using buffered DNA strand displacement gates. Briefly, a buffered DNA strand displacement gate is one where the gates are initially present in an inactive form and must be activated by a precursor input strand.24 This extension enables the inactive gates to be present in excess and called up as required to carry out computation. By designing the circuit so that availability of these activating strands can be controlled by another gate, we can implement buffered gate networks with feedback loops that adjust future circuit behavior in response to past stimuli. As a proof of concept, we provide an example circuit that uses buffered strand displacement gates to implement supervised learning of a class of linear functions using stochastic gradient descent, which is a well-known and well-studied machine learning algorithm.25 We use computational simulations to demonstrate that the DNA circuit correctly implements this learning algorithm and to characterize its learning performance for a range of initial conditions and learning targets as well as in the presence of noise in the inputs. Thus, this article offers a route to an experimental realization of adaptive molecular computing systems using DNA strand displacement. (A preliminary version of this article appeared in the DNA21 conference.26)

bias toward completion of the displacement reaction because the resulting invader−template duplex is a lower-energy structure than the initial incumbent−invader duplex since more base pairs are formed. Toehold exchange is a generalization of the basic strand displacement process in which the gate contains a secondary toehold that is initially sequestered in the incumbent−template duplex and is not directly displaced by the invader strand.9 This means that the incumbent strand must spontaneously unbind when the branch migration process is complete. This secondary toehold provides a pathway for a displaced incumbent strand to rebind to the gate and displace the original invader, thereby providing a means to implement reversible strand displacement reactions. Crucially, a displaced incumbent strand that is the “output” from one gate may diffuse and bind to another gate elsewhere in the system, thereby serving as its “input”. For this to take place, it suffices that the nucleotide sequences of the strands should match correctly. Thus, we can straightforwardly connect strand displacement gates into circuits by a judicious choice of nucleotide sequences. The foregoing principles of DNA strand displacement can be used to create computational gates that accept multiple inputs and produce multiple outputs, and they can be used to connect these gates into networks that perform nontrivial information processing tasks. Previous experimental work has produced strand displacement cascades that mimic the behavior of digital logic circuits11,12 and catalytic cycles.13 A distributed consensus algorithm initially specified as an abstract chemical reaction network was also recently implemented using DNA strand displacement,14 and antibody-conjugated DNA strand displacement cascades have been used to analyze cell surface markers expressed by individual blood cells.15 Those results demonstrate the unique capabilities of DNA strand displacement as a programmable chemical framework. However, all published experimental DNA strand displacement systems have the drawback that they are “one shot” computational devices: an experimental system is prepared and is perturbed by the introduction of input strands, and the reactions within the circuit compute the corresponding output as the system evolves toward equilibrium. Once equilibrium is reached and the circuit output has been read out, typically as a fluorescent signal, that collection of molecules cannot be used to carry out another computation and must be discarded. Thus, these molecular computing systems cannot be used in contexts where the circuit must remain operational to process multiple presentations of the input signal(s) over an extended period of time. In particular, this drawback precludes the implementation of adaptive molecular algorithms that can modify their behavior over time in response to time-dependent input signals. For example, a previous experimental implementation of artificial neural networks using DNA strand displacement16 required the network to be trained in silico, and each prepared instance of the experimental system could be used only once. Implementing adaptive molecular information processing systems is a key challenge for the fields of molecular computing and synthetic biology because such systems are ubiquitous in nature. Adaptive capabilities are critical to enable cells (or synthetic protocells) to respond appropriately to changing environmental conditions and maintain internal homeostasis. Processes such as bacterial chemotaxis are regulated in an adaptive manner by simple desensitization of chemical sensors,17 and associative learning has been demonstrated in paramecia18 and sea slugs.19 Previous theoretical work has shown that

1. RESULTS 1.1. Buffered DNA Strand Displacement Gates. In this section, we present the framework of buffered strand displacement gates that we will use below to implement adaptive algorithms. The basic concept of buffered strand displacement gates was introduced by Cardelli.27 The key idea is that buffered gates that encode a given abstract chemical reaction, say X → Y, are initially present not in an “active” configuration that can immediately accept the input X but rather in an “inactive” configuration that requires an initial step to activate a subset of the inactive gates so that they can accept the input X. We call this pool of “inactive” gates the buf fer and refer to the strand that initiates the activation reaction (also referred to as the “unbuffering” reaction) as an unbuf fering strand. The unbuffering reaction is typically made into an irreversible process by the addition of a sink gate that irreversibly consumes the strand displaced by the unbuffering strand. As a matter of for a buffered reaction gate notation, we will write that emulates the abstract chemical reaction R̅ → P̅ and whose unbuffering reaction is controlled by the unbuffering strand B. B

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

Figure 1. Example of a buffered four-domain DNA strand displacement gate, implementing buffered reaction . (a) The first reaction is the unbuffering reaction, in which an unbuffering strand B binds to a buffered gate and reversibly produces an active input gate. The binding toehold for the X reactant strand is exposed in the active input gate. (b) To make the unbuffering reaction effectively irreversible, an additional sink gate for the strand released by the unbuffering reaction is provided (this gate is not required in standard four-domain chemical reaction network encodings). (c) Reactant strand X can bind to an active input gate via the Xa^ toehold, irreversibly displacing an intermediate strand. (d) The released intermediate strand binds to the output gate and releases all of the products in a single-strand displacement reaction. Here, in addition to product strands X and Y, an additional output, B, is generated, which is a new copy of the unbuffering strand for this gate. Thus, a new gate will be activated from the buffer to replace the gate that was consumed to execute this reaction.

We will use the ⌀ symbol to denote an empty list of reactants or products. The advantage of using buffered strand displacement gates is that, because the inactive gates in the buffer do not interact directly with the input strands, the buffer can contain an excess of inactive gates without affecting the kinetics of the reactions that actually consume the strands representing abstract reactants and produce the strands representing abstract products. The kinetics of the emulated reactions are governed (in part) by the quantity of active gates, that is, by the amount of the unbuffering strand that is initially added. Buffered strand displacement gates are also typically designed so that, in addition to producing product strands, an additional copy of the corresponding unbuffering strand is also released into solution. This means that the population of active gates remains approximately constant over time, so the emulated reaction kinetics remain constant even as the population of gates in the buffer is consumed. This property was used in previous theoretical work to show that buffered strand displacement gates can be used to implement a three-phase DNA oscillator system with more robust long-term kinetics than the corresponding unbuffered circuit.24 Furthermore, the buffer can be replenished periodically without causing a corresponding shift in the emulated reaction kinetics, and if waste products are also removed to

prevent them from accumulating, then the buffered reaction scheme can allow long-running strand displacement computations to be carried out without the additional complication and expense of an open microfluidic system to continuously provide active gates to the system. As a concrete example of a buffered strand displacement that uses gate, Figure 1 presents a buffered gate the reactant strand X to catalyze production of the product strand Y under the control of unbuffering strand B. This example, and indeed the remainder of this article, uses the “four-domain” encoding of abstract chemical reactions developed by Soloveichik4 as a basis for developing a buffered strand displacement scheme. In four-domain encoding, each abstract species X is represented by three domains Xa^, Xb, and Xc^ arranged from 5′ to 3′ in a single strand of the form ⟨h Xa^ Xb Xc^⟩, written here in the syntax of the DSD language.24,28 Briefly, in the DSD language, domains are identified as toeholds using the circumflex symbol and strands are written as a sequence of domains enclosed in angle brackets. The domain h at the 5′ end of the strand is called the “history domain” because it simply records which gate generated that particular strand and plays no role in downstream reactions. Thus, the total “population” of strands representing abstract species X is the sum of the populations of strand of the form ⟨h Xa^ Xb Xc^⟩ for all history domains h. C

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology We may write “?” for a “history domain” whose identity is irrelevant for the purposes of the reaction under discussion. We note, however, that while the identity of the history domain in a given strand is irrelevant to downstream reactions it is the case that history domains for product species must be distinct for each product-producing gate, to avoid crosstalk. Finally, we note that, although here we use the four-domain reaction encoding scheme the basic idea of buffered strand displacement is not tied to any particular reaction encoding scheme and buffered gates could, in principle, be used with any available encoding. Indeed, our previous work on buffered oscillators24 used the three-domain encoding developed by Cardelli.5 1.2. A Buffered Architecture for Adaptive Strand Displacement Networks. In this section, we demonstrate how buffered strand displacement gates can be used to implement adaptive molecular circuits using the simple example of a signal amplifier motif whose gain can be dynamically adjusted. This circuit design motif will be a critical component of the proof-ofprinciple learning circuit design that we will present below. The key property of buffered gate systems is that the kinetics of gate reactions are controlled by the size of the activated gate population, which is directly related to the quantity of unbuffering strands introduced into the system. Releasing fresh unbuffering strands into solution to replace those that were used to activate a gate that has been consumed in the execution of a reaction produces stable, long-term kinetics, as mentioned above. However, here we do not necessarily want reaction kinetics to stay constant over time; rather, dynamically altering reaction kinetics is a powerful and convenient way to implement adaptive systems. In a buffered strand displacement system, we can dynamically increase the rate of a buffered reaction by releasing additional free copies of the corresponding unbuffering strand into the system. (This is a simpler approach to control reaction kinetics than engineering toehold binding energies or developing remote toehold systems.29) Since the unbuffering strands and the strands that encode abstract reaction species have the same domain-level structure, as shown in Figure 1, a given buffered gate may release additional unbuffering strands for other buffered gate populations as part of its output. Thus, one buffered gate can adjust the number of active gates of a second kind, thereby controlling the rate of the second buffered gate’s reaction. 1.3. An Adaptive, Buffered Amplifier. As a concrete example of this approach to implementing adaptive molecular circuits, we now illustrate the use of buffered strand displacement gates to implement an adaptive amplifier whose gain can be dynamically adjusted. Previous work on DNA strand displacement-based amplifiers12,13,16,30,31 has relied on hardcoding the gain of the amplifiers in the initial species concentrations and provides no consideration to reusability or autonomously adjusting the gain of the amplifier. The circuit design motif introduced here will demonstrate a buffered strand displacement system whose result can be controlled by adjusting the provision of unbuffering strands and will also be a key component of our learning circuit design. Our design for a reusable, adaptive, buffered strand displacementbased amplifier consists of two buffered strand displacement gates:

output species y, i.e., x effectively catalyzes production of y. The second gate is a “degradation” gate that removes input strand x from the system by consuming it without releasing any output species (except for a new activating strand B′). A graphical shorthand for this circuit design element is presented in Figure 2a.

Figure 2. Buffered amplifier design and ODE simulation data. (a) The buffered amplifier consists of two buffered reactions: a catalytic reaction that generates the output and a degradation reaction that removes the input from the system. The ratio of the concentrations of active gates from the two buffers corresponds to the gain of the amplifier. (b) Data from an ODE simulation of the buffered amplifier with multiple sequential input additions and dynamic control of amplifier gain. The initial concentration of B′ was 100 nM, and the initial concentration of B was 200 nM, which sets the gain, w, to be 2.0. Addition of x input species (1 nM) at t = 100 s causes the amplifier to generate output species y at a concentration of 2 nM, as expected. At t = 600 s, an additional 300 nM of unbuffering strand B was added, increasing the amplifier’s gain (w) by 3.0. Then, subsequent addition of x (1 nM) at t = 700 s produced a further 5 nM of output species y, showing that the gain had indeed been increased to 5.0, as expected. (Note that the plotted value of w is not a concentration but rather a ratio of concentrations, but the value follows the left axis.)

This gate design uses competition to achieve amplification of the input signal by the desired gain factor, as the ratio between the effective rates of the first and second gates controls the number of output molecules y that each input molecule x can produce before it is irreversibly consumed; therefore, this ratio controls the gain of the amplifier. Thus, in our buffered amplifier design, the ratio of the concentrations of the active catalyst and degradation gates, which is also the ratio of the concentrations of the corresponding unbuffering strands, controls the gain of the amplifier. We shall use the notation [z]0 for the initial concentration of species z and [z]ss for the steady-state concentration of species z (when a steady state exists). Then, it is easy to show that [y]ss =

[B]0 × [x]0 [B′]0

We typically let w =

[B]0 [B ′]0

(3)

and refer to w as the gain or weight of

the amplifier. To illustrate the dynamic gain control capability of this buffered amplifier motif, we carried out ordinary differential equation (ODE) simulations of the two buffered reactions shown above. We used buffered four-domain strand displacement gates, as shown in Figure 1, to implement this circuit motif and simulated it using the beta version of the DSD compiler that includes support for mixing events.32 The results from this simulation are presented in Figure 2b and show that the amplifier produces correct results for several gain settings and that the gain of the amplifier can be dynamically adjusted

The first gate is a “catalyst” gate that accepts input species x as its reactant and produces another copy of x and a copy of D

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

enable the representation of both positive and negative signal values while avoiding the need to insert an annihilator gate to cancel out the species that represent the positive and negative components of each signal in the system. Using annihilator gates is the standard approach to this issue;32,33 however, as we elaborate in the Discussion, that approach can cause problems with sequestration of the signal species by the annihilator gate, which complicates analysis of circuit results. Therefore, here we accept the simultaneous presence of both positive and negative components of each signal and simply use the differential encoding to recover the true represented value of the signals. The predictor subcircuit design is presented in Figure 4; the motif presented therein is replicated once for each input xi.

by directly adding more unbuffering strands before adding input strands to trigger more amplification reactions. 1.4. A Strand Displacement Learning Circuit. In this section, we will present a buffered strand displacement system that can learn linear functions f of the form f (x1 , x 2) = w1 × x1 + w2 × x 2

(4)

where w1 and w2 are real-valued coefficients and x1 and x2 are real-valued inputs. In this article, we restrict ourselves to the two-input case for clarity, although the circuit design that we will present could be extended to process additional inputs. In particular, a bias term could be included by incorporating an additional input signal x0 that is always supplied with the input value x0 = −1 in each training round, as is common in studies of artificial neural networks. We will present a strand displacement system that learns functions of this form using a stochastic gradient descent algorithm.25 Gradient descent is a general solution for many optimization tasks, in which the current weight approximations are adjusted to minimize the squared error over the entire training set in each training round. Stochastic gradient descent is a simplification of gradient descent that considers only a single training instance in each training round, making it more amenable to implementation in a molecular computing system. A molecular computing system that solves this problem must accept input values x1 and x2, compute the predicted output y = w1 ̂ × x1 + w2 ̂ × x 2 based on the current stored weight approximations w1 ̂ and w2 ,̂ compare y with the supplied expected value d = w1 × x1 + w2 × x2, and update the stored weight approximations according to the gradient descent weight update rule wi ̂ ≔ wi ̂ + α × (d − y) × xi

Figure 4. Design schematic for the predictor subcircuit. The concentrations of input species x+i and x−i (for each input signal xi) are copied by buffered fork gates to produce species that serve as inputs to the buffered amplifier motifs that implement that linear function prediction and additional species (ki1±, ki2±) that will be used in the feedback subcircuit. The gains of these amplifiers store the current weight approximations. The amplifiers produce species y+ and y−, such that [y+] − [y−] equals the predicted output value. These species then interact with the d+ and d− species, which encode the expected output value, via four buffered reactions that implement subtraction. When the predictor subcircuit reaches steady state, the concentration difference [d+] − [d−] should equal d − y, i.e., the error in the prediction when compared with the expected output value.

(5)

where α is a (usually small) positive coefficient called the “learning rate”, which controls the size of weight updates. This process can be visualized as a flowchart, as shown in Figure 3.

Here, and henceforth, we will omit the identities of buffer species from figures if their identity is not important when describing the operation of the circuit. The species highlighted in gray (xi+, xi−, d+, and d−) are provided by the user to initiate each training round; their concentrations represent input value xi and expected result d that the user must derive using the target weight values. (If d = 0, then we must add equal concentrations of d+ and d−. Any non-zero concentration is acceptable, as these species must be present to drive the execution of the predictor subcircuit.) Each positive input signal xi+ is “copied” by a buffered fork gate that accepts the input signal and generates four output signals, each of which will have the same overall concentration as the original at steady state. Two of these outputs, xi1+ and xi2+, are for use by the predictor subcircuit, and the remaining two, ki1+ and ki2+, are for use by the feedback subcircuit (as detailed below). Each negative input signal xi− is copied similarly. The key parts of the predictor subcircuit are the buffered strand displacement amplifier motifs. The initial gains of the predictor subcircuit amplifiers encode the initial approximation of each weight value stored in the system (step (1) from Figure 3).

Figure 3. Flowchart illustrating the basic operations carried out by our learning system.

We now present a strand displacement system that implements this learning algorithm using buffered DNA strand displacement gates. Our design can be divided into two subcircuits, as follows. Predictor Subcircuit. The predictor subcircuit comprises steps (2)−(4) from the flowchart in Figure 3 and is based on the adaptive, buffered amplifier motif presented above. The input signals, and all other numeric signals, are represented in a dual rail format using a differential encoding, that is, input signal xi is actually represented by two species, xi+ and xi−, and the value of xi is interpreted as the concentration difference [xi+] − [xi−]. This encoding of numerical signals was used to E

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

and the expected output, the feedback subcircuit must use the value of this discrepancy to update the weight approximations stored in the predictor subcircuit, according to the gradient descent learning rule. This corresponds to step (5) from Figure 3. A design schematic for our feedback subcircuit is presented in Figure 5.

There is one pair of amplifiers per positive input signal and one pair per negative input signal. In each of these pairs, the gain of one amplifier represents the positive component w+ ̂ of i

the corresponding weight approximation wi ,̂ and the gain of the other amplifier represents the negative component wi−.̂ Thus, the positive component of each input is multiplied by both the positive and negative components of the corresponding weight and similarly for the negative component. The action of the amplifiers in the predictor subcircuit is to convert the xi1± and xi2± input species into output species y+ and y−, which represent positive and negative components of the current prediction at this stage in circuit execution. The concentration of the output species that is generated in each case is the corresponding input concentration multiplied by the corresponding stored weight approximation, and the amplifier gates are wired to their outputs such that the sign of the output species is correct with respect to the signs of the input component and the weight component in each case. To complete the execution of the predictor subcircuit, the y+ and y− species, whose concentrations represent the current prediction, interact with the d+ and d− species. The initial concentrations [d+]0 and [d−]0 of the d+ and d− species represent the value of the expected output from the target function when presented with the input values that were provided in the training instance. In the final part of the predictor subcircuit, the y+, y−, d+, and d− species interact via the four buffered reactions shown in the box on the right-hand side of Figure 4. These reactions have the collective effect of subtracting the value of y from the value of d and storing the resulting value in the d signal, i.e., in the steady-state concentration difference [d+]ss − [d−]ss. We implement this operation via four two-input, two-output reactions, in which the d± species catalyze conversion of the y± species to d±, with signs chosen such that the resulting concentrations of d± represent the result of a subtraction. We choose to implement subtraction in this way, rather than using annihilator gates that degrade the positive and negative variants of the species, to avoid sequestration of the remaining species by the annihilator gates, which has been problematic in other work.32 For this reason, it is crucial that the y± species must be the first input strands consumed by these reaction gates, and to ensure correct operation of the subtraction reactions, it is necessary for the result of the subtraction to be stored using the d± species. This is because the additional d± species generated by conversion of the y± species produced by the amplifiers may themselves be needed to catalyze conversion of additional y± species. Thus, when the predictor subcircuit reactions reach steady state, the concentration difference [d+]ss − [d−]ss represents the value of d − y that must be computed according to the algorithm summarized in Figure 3. (Recall that the initial concentration difference [d+]0 − [d−]0 stored the value of d, but it was updated by the subtraction operation of the predictor subcircuit so that [d+]ss − [d−]ss stores the value of d − y at steady state.) If the value of d − y is positive, i.e., if [d+]ss > [d−]ss, then the predicted output value was too small. Conversely, if the value of d − y is negative, i.e., if [d+]ss < [d−]ss, then the predicted output was too large. The overall goal of the learning process is to adjust the stored weight approximations wi ̂ to match the target weight approximations so that [d+]ss = [d−]ss when the predictor subcircuit reactions reach steady state. Feedback Subcircuit. Once the predictor subcircuit has computed the discrepancy between the predicted function output

Figure 5. Design schematic for the feedback subcircuit. (a) Graphical shorthand for a multiamplifier motif that enables one input concentration [x] to be multiplied by another concentration [k] and a scalar scaling factor β. (b, c) Design schematic for the feedback subcircuit. Here, scaling factor β = α × δ, where α is the learning rate constant and δ is the denominator of the weight ratios in the predictor subcircuit. The feedback circuitry from (b) is activated when d+ is left over from the predictor subcircuit, and the feedback circuitry from (c) is activated when d− is left over from the predictor subcircuit.

The first point to note from the gradient descent learning rule is that the feedback subcircuit must take the concentration of d± that denotes the discrepancy from the predictor subcircuit and, for each input, multiply the discrepancy value by the corresponding input value, both of which are concentrations, and by the learning rate constant α. A single buffered amplifier can multiply an input concentration only by a gain factor encoded as a ratio of concentrations. To enable two input concentrations to be multiplied together, we have developed a two-concentration multiplier circuit motif, shown in Figure 5a. In this motif, we assume that no unbuffering strands are initially present for the uppermost amplifier, which will accept input signal x and produce the output species. Input signal k activates an amplifier that produces the unbuffering strand for the catalyst gate from the output-producing amplifier, with gain β. An additional input signal, whose value is constant 1, activates an amplifier that produces the unbuffering strand for the degradation gate from the output-producing amplifier, with gain 1. Thus, these secondary amplifiers preset the gain of the output-producing amplifier to be β × [k], and upon addition of input signal x, the resulting concentration of both output species (y and z) will be β × [x] × [k], as desired. We include two outputs because this variation of the two-concentration F

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

amplifiers from the feedback subcircuit will be primed with the correct gain values before it starts executing. After a further 2000 s, when the predictor subcircuit has completed its execution, we add unbuffering strands for the fork gates from the feedback subcircuit, which triggers execution of the feedback subcircuit. After a further 3500 s, when the feedback subcircuit has finished updating the weight approximations stored in the predictor subcircuit, we set the concentrations of any remaining unbuffered (active) fork gates and output-generating gates in the feedback subcircuit to zero, to reset the state of the feedback subcircuit. We also add the second set of training inputs at this time and iterate until the entire specified sequence of training instances has been presented. We began by investigating the learning performance of the learning circuit for a series of randomly chosen initial and target weight values and randomly chosen training instances. We ran a total of 1000 simulations, with initial and target weight values and input values selected from a uniform random distribution over the interval [−10, 10] and with a fixed learning rate value α = 0.01. Each simulation used a 15-round training schedule, with each training input also selected from a uniform random distribution over the interval [−10, 10]. For our two-input learning circuit, the weight space can be visualized in a twodimensional Cartesian coordinate space, along with the trajectory of how the weights evolve, from their initial values (labeled “Start”) toward the true values (labeled “Target”) over the training period. Figure 6 presents some representative traces that show the evolution of the weight approximations from their starting values toward the target values over the course of their 15-round training schedules. In each plot, the weight trajectory from an ODE simulation of the DNA learning circuit design (red) is overlaid with a trajectory computed using a reference implementation of the learning algorithm from Figure 3, using the gradient descent learning rule with the same initial state and the same training schedule. The close agreements between the weight trajectories for the DNA circuit and the reference implementation shown in Figure 6 are representative of the good agreement between the DNA system and the reference implementation observed in all cases. Furthermore, the weight trajectories correctly home in on the target weight values. These results give us confidence not only that the DNA circuit is capable of learning but also that it is a correct implementation of the stochastic gradient descent learning procedure. We investigated the aggregate learning performance of the system by computing the root-mean-square error (RMSE) of the stored weight approximations with respect to the target weight values at each step of the training schedule for each of the 1000 training simulations. These RMSE values, averaged over the 1000 training simulations, and their standard deviations are plotted in Figure 7. Again, in this “learning curve”, the results from the DNA circuit simulations and the reference implementation are overlaid precisely, indicating that the DNA circuit behaves as specified by the algorithm. The average RMSE values tend toward zero, which shows that the stored weight approximations in the system tend to improve (that is, more closely approximate the target weight values) with the presentation of additional training instances. Furthermore, the RMSE in the weight approximations is almost zero after the 15 training rounds, which shows that the system is capable of learning the target weights well. This also suggests that an experimental implementation of this system could be trained to a reasonable degree of accuracy in a limited time frame.

multiplier is called for in the feedback subcircuit of our twoinput learning circuit design. The feedback subcircuit uses the two-concentration multiplier circuit motif extensively to produce a molecular implementation of the gradient descent weight update rule. Figure 5b,c presents the feedback subcircuit design for the twoinput case. Execution of the feedback subcircuit is initialized by buffered fork gates that copy the d± species into four species da±, db±, dc±, and dd±. We assume that there are initially no unbuffering strands for these fork gates, and once the predictor subcircuit has run to completion, the addition of these unbuffering strands triggers execution of the feedback subcircuit. For each combination of signs for the leftover d± species and the copied input species kij± from the predictor subcircuit, the copied d± species and the copied kij± species serve as inputs to an instance of the two-concentration multiplier circuit motif. As described above, this circuit motif multiplies these values together (and by a scaling factor β = α × δ, where α is the learning rate constant and δ is the denominator of the weight ratios in the predictor subcircuit) and generates additional unbuffering strands for certain amplifier gates from the predictor subcircuit. From the two-input predictor subcircuit design from Figure 4, we see that additional unbuffering strands Bia+ and Bib+ must be generated to increase weight approximation wi ̂ and that additional unbuffering strands Bia− and Bib− must be generated to decrease weight approximation wi .̂ These pairs of species must be generated together so that the pairs of amplifiers from the predictor subcircuit that are initialized with the same weight approximation are updated together. Finally, we note that the Supporting Information contains a model of our learning circuit expressed in the DSD molecular programming language24,28 which formally defines the structure of our molecular learning circuit design. 1.5. Results from Simulations of the Learning Circuit. In this section, we present the results from simulations in which we investigated the learning capabilities and performance of the learning circuit design presented above. We used buffered four-domain strand displacement gates, as shown in Figure 1, to translate our circuit design into a concrete DNA strand displacement system, and we refer the reader to the Methods section for details of how our simulation models were constructed and of how the simulations were executed. Briefly, the learning circuit reactions were compiled using Visual DSD28 and simulated using MATLAB in such a way that “mixing events” could be defined and carried out at certain time points during the simulation. These mixing events simulate the addition of inputs to the system, or the removal of species, by the experimenter at those time points. The initial state of the learning circuit consists of the various buffered gates and their unbuffering strands (with the exceptions of the unbuffering strands for the fork gates and output-generating amplifiers in the feedback subcircuit, as these are added to initiate the feedback computation, as described below). We used an initial buffer size of 1 × 107 nM. The initial weight approximations wi ̂ are encoded as the gain settings of the amplifiers in the predictor subcircuit. After the gates have unbuffered (we wait 500 s for this to take place), the first training inputs are added, which consist of the xi± and d± species, whose concentrations encode the input values and the expected function output, respectively. At the same time, we also add the constant-valued signals that serve as an input to help prime the feedback subcircuit so that the output-generating G

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

Figure 6. Example weight traces from DNA learning circuit simulations (red solid lines), overlaid with corresponding traces from our reference implementation (black broken lines). Each training schedule was 15 rounds in length. Initial and target weights were drawn from a uniform distribution over [−10, 10], as were the input values for each training instance. These traces are representative of the other simulations that we performed, indicating that the DNA circuit functions as intended.

for a range of initial weight values distributed throughout the weight space to build up a picture of the performance of the learning circuit across the weight space. Target weight values were chosen from the interval [−10, 10] at regular intervals of 1 unit, and initial weight values were chosen on the diagonals of the weight space. Results from weight space scanning simulations for several initial weight values are presented in Figure 8. Similar results for additional initial weight values are presented in the Supporting Information, along with plots that illustrate the weight RMSE values across the weight space at different points in the training schedule, which show how the learning proceeds over time. The results from Figure 8 and the Supporting Information show that the final RMSE values are low throughout the weight space. Furthermore, the final RMSE values are independent of the starting weight values. This indicates that the system can learn different weight values without inherent bias. The occurrence of isolated points with higher RMSE may be attributed to buffer exhaustion, which is discussed further below. 1.6. Investigating the Effect of Noise on Learning Circuit Operation. There are numerous potential sources of error in experimental implementations of molecular computing systems. One of the most prominent is inaccuracies in species concentrations due to pipetting error. To study the effect of such inaccuracies on the performance of the learning circuit, we ran additional simulations in which the species concentrations added in the mixing events that corresponded to training instances, that is, the values of x1, x2, and d were each perturbed by a (different) noise term, drawn from a normal distribution with mean 0 and standard deviation σ. For a range of values of σ,

Figure 7. Learning curves from DNA learning circuit simulations (red) and from our reference implementation (black). The average RMSE in the weight approximations was computed after each training round for 1000 training schedules, each 15 rounds in length. Initial and target weights were drawn from a uniform distribution over [−10, 10], as were the input values for each training instance. The error bars are 1 standard deviation above and below the mean. The data from the DNA circuit and reference implementation coincide well, indicating that our DNA circuit design functions as intended.

The learning curve gives a good measure of aggregate learning performance; however, it may obscure any biases in the learning procedure that make it more challenging for the system to learn target weights in particular areas of the weight space or to adjust the stored weight approximations in particular directions through the weight space. To test our learning circuit design for any such biases, we ran a series of additional simulations in which the initial weight values were fixed and the system was trained to learn a number of target weight values evenly distributed across the weight space. This procedure was repeated H

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

Figure 8. Scatter plots showing the distribution of learning performance across a range of target weight values for several different initial weight values. The position of each marker represents a pair of weight values in the two-dimensional weight space, and each marker is color-coded to represent the final RMSE in the weight approximations when attempting to learn those weight values from the weight values marked “Start”. In each case, the average RMSE in the weight approximations was computed after a 15-round training schedule using input values for each training instance drawn from a uniform distribution over [−10, 10]. We observe a low and uniform distribution of final RMSE values, indicating that the system can learn a range of different weight values without inherent bias.

To summarize, these results illustrate that our learning circuit design can still function tolerably well in the face of noise in the training input values. We refer the reader to the Supporting Information for additional results from these simulations in the form of histograms showing the distributions of final weight RMSE values for different values of σ.

we ran 1000 simulations using the same collection of initial and target weights and the same 15-round sequences of training instances as those used in Figure 7. The results from these simulations are summarized in Figure 9, which plots the final weight RMSE values after 15-round training

2. DISCUSSION We have presented buffered DNA strand displacement reactions as a framework for the implementation of adaptive molecular computing systems and used this framework to implement a strand displacement circuit to learn a class of linear functions using the well-studied stochastic gradient descent learning algorithm. We demonstrated, via in silico simulations, that this circuit correctly implements the desired learning algorithm and can learn target weight values in a reasonable time frame. Furthermore, we showed that there is no overall bias for learning different target weight values and that the circuit can function in the face of certain amounts of noise in the training inputs. These findings demonstrate that our molecular learning circuit design has the potential to be realized in a practical, experimental implementation, and we reserve such experimental implementations for future work. The study of learning in (synthetic) molecular systems is a relatively new endeavor; hence, the list of related work is not extensive. Our own prior theoretical work in this area has demonstrated that a biochemical system assuming hypothetical DNAzyme-like reactions can learn a class of linear functions.34 The learning system presented in the current article surpasses this previous work by allowing negative weight values and input values and also alleviates problems with our previous approach that caused poor performance when trying to learn weight values near zero. This is because the stochastic gradient descent learning procedure used in this work, which is well-known and has been extensively studied,25 is symmetric in the positive and negative weight update directions. Furthermore, it is well-known that the learning rate parameter (α) can have a significant effect on the rate of convergence of the learning process, and the fact that our learning circuit is founded on the well-studied gradient descent scheme means that we can predict the effects of modification to our learning scheme, such as the use of procedures

Figure 9. Learning performance in the face of increasing noise in the training inputs. The average RMSE in the weight approximations was computed after completing 1000 training schedules, each 15 rounds in length. The same initial and target weights and training instances were used as those in Figure 7. The error bars are 1 standard deviation above and below the mean. The data indicate that our learning circuit design can learn even in the face of moderate amounts of input noise.

intervals (averaged over all 1000 simulations in each case) against the standard deviation σ of the noise distribution. The values for σ = 0 are the same as the final point (after 15 training rounds) from Figure 7. We observe that the final average values of the RMSE in the stored weight approximations increases with increasing amounts of noise, that is, with increasing values of σ. This is to be expected; however, on the lower end of the noise spectrum (say, with σ ≤ 2), we still observe good learning performance. This suggests that the learning circuit could perform reasonably well in the face of moderate amounts of input noise, which is an important property for future experimental realizations of our learning circuit design. The high standard deviation of the final weight RMSE values, as reported by the error bars in Figure 9, shows that with high amounts of noise the end results of training are highly variable. Again, this is to be expected, as with large amounts of noise the signal being presented by the trainer is obscured and becomes lost in the noise. I

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

strand displacement gate and reveals a secondary toehold that enables a second reactant strand to bind, and so on. The reactant strands typically bind via reversible toehold exchange reactions9 so that the reactants cannot be consumed irreversibly until all required reactants have been consumed by the gate. However, reactants are still temporarily consumed. For example, for a two-input reaction gate, if the first input is present but the second input is not, then a fraction of the first input population will always be temporarily bound to, or sequestered by, the gate. The original version of Soloveichik’s compilation scheme dealt with this issue by adding additional reaction gates to compensate by sequestering all other species similarly, thereby slowing the whole system. Other work32 has dealt with input sequestration by artificially increasing the values of certain input signals. Our desire to avoid such considerations led to the design of the “subtractor” circuit motif shown in Figure 4, which avoids sequestration of the d± species by using it as the second input to each fork gate. If the d± species were the f irst input, then the resulting d± produced by these gates would be sequestered. This solution, while somewhat elegant, has the property that the absolute values of the positive and negative components of the weight approximations and other signals processed by the system grow monotonically over the course of a multiround training session. This can drain other species out of the system by increasing the amount of fuel species that must be used to process subsequent input signals. This problem could be ameliorated by the inclusion of a single annihilator gate in the feedback subcircuit that takes d+ and d− as inputs but produces no output, which could reduce the absolute sizes of these signals. This annihilator gate would need to compete only with the rest of the feedback subcircuit to drain some of the excess d± signals form the system. By employing just a single annihilator gate at a specific part of the circuit, this approach would strike a balance between our approach (with no annihilator gates at all, which leads to monotonically increasing concentrations of the species representing dual rail signal components) and the standard approach to dual rail signal representation in molecular circuits (with annihilator gates added for every dual rail signal, which increases gate counts and introduces more competition for every signal strand). Thus, we would hope to ameliorate the problem causes by monotonically increasing species concentrations without introducing too much additional competition into the system, which could compromise the accuracy of the steady-state results obtained from our circuit. We will explore this alternative and its effects on the performance and accuracy of the learning circuit in future work. Our use of buffered gates to implement adaptive strand displacement circuits allows us to implement systems whose operation can be extended by the periodical replenishment of buffered gate species before they are completely depleted, without adversely affecting the kinetics of the system. This was demonstrated in previous work on the design of oscillator networks using buffered strand displacement reactions.24 The ability to resupply a circuit with additional input reagents is an important consideration for the implementation of molecular learning circuits, which would typically need to be active for an extended period of time. In particular, in an in vivo application of a molecular circuit it may not be possible to provide a continuous stream of new reagents and periodical replenishment, e.g., via a daily dose of some therapeutic containing the reagents, may be the only option. In our simulations, we found that the size of the buffer, that is, the initial quantity of

for reducing the learning rate over time to achieve reliable convergence. Learning behavior has also been studied in high-level artificial chemistries, which have been shown capable of learning to implement Boolean functions,35 perceptron-like classification tasks,36 analog functions,37 and temporal tasks.38 Other groups have also studied associative learning in the context of simulated transcriptional networks20 and evolutionary computational experiments on chemical reaction networks.21 Associative learning is a learning scheme by which simultaneous triggering of inputs strengthens the connection between those inputs, and in this way, the system learns that they are associated. This is a different form of learning from that studied in this article, and it would be interesting to consider whether a strand displacement-based approach similar to ours could be adapted to implement an associative learning circuit. Another form of learning that may be well-suited to compact implementations in molecular systems is “winner-takes-all”, in which multiple agents compete for a shared catalytic resource and those that do not succeed in this competition are eliminated from the system. Designs for winner-takes-all systems have been proposed in in vitro transcriptional networks39,40 and strand displacement systems.40 From a circuit design perspective, our amplifier design draws on previous work on DNA strand displacement-based amplifiers13 and is inspired in particular by the strand displacement classifier circuit introduced by Zhang and Seelig.31 Our contribution is to extend this idea to buffered strand displacement systems that can be reused to classify multiple input signals and to link it to a feedback subcircuit so that learning can take place. Our approach to implementing the amplifier circuit motif is also related to the “ideal gain blocks” developed by Oishi and Klavins,33 with the main difference being that our approach produces a quiescent final state, whereas their approach produces a steady state in which production and degradation of the output species balance. Thus, our approach prevents the supply of input species being constantly drained, as in Oishi and Klavins. However, our approach does require some care to be taken when composing the output from an amplifier circuit motif with downstream circuit elements, as we outline below. Nevertheless, the adaptive strand displacement amplifier motif that we have developed is of practical interest in its own right, and we anticipate that implementing this component in the laboratory will be an important direction for future research. Our work is also related to previous work on DNA-based feedback and control systems32,41 and on (nontrainable) artificial neural networks implemented using DNA strand displacement.16 Finally, although we have assumed components and reactions based on DNA strand displacement in our circuit design, it is worth pointing out that learning behavior could, in principle, be realized in synthetic biomolecular systems based on different chemical frameworks, such as DNAzyme circuits,42−44 DNAzyme−strand displacement hybrids,45−47 ex vivo transcriptional circuits,39,40 and protein-based molecular computing systems.48−50 When designing our learning circuit, we took care to avoid the issue of input sequestration, which can be problematic in DNA strand displacement systems. In all published compilation schemes for abstract chemical reactions into DNA strand displacement reactions, such as that proposed by Soloveichik4 and used here to create the DNA model of our learning circuit, reactions with more than one abstract reactant are implemented by multistep processes in which one reactant strand binds to a J

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology

buffered DNA strand displacement network. This circuit design was programmed in the DSD domain-specific programming language for DNA strand displacement systems.28 The reader is referred to the Supporting Information for the DSD source code listing. This code was compiled using the DSD compiler with the “Inf inite” reaction semantics selected24 to produce the corresponding chemical reaction network model. The initial state of the two-input system consists of 278 species, of which the majority (222) are gate complexes. It is worth noting, however, that many of these species are variants with different combinations of history domains, so the design complexity of the circuit is not as high as that suggested by the raw species counts. Furthermore, the number of species should scale linearly with the number of input signals, so the circuit design presented here could be replicated to learn similar functions of more than two inputs. 3.2. Execution of Learning Circuit Simulations. The compiled DSD model was exported as a MATLAB53 code file that implements an ODE model of the chemical reaction network. We simulated the ODE model of the learning circuit using a custom MATLAB simulation routine that invokes a stiff ODE solver and that also allows mixing events to be executed at certain time points during the simulation. These mixing events simulate the addition of inputs to the system or the removal of species by the experimenter at certain time points.

unbuffered gates waiting to be activated, did not affect the accuracy of the computation performed by the learning circuit. It did, however, govern the number of training rounds that could be conducted before the buffer was completely exhausted, at which point the computation ceased to function correctly. The Supporting Information contains additional results showing that, if the training schedule is extended without a concomitant increase in the initial quantity of buffered gates, we observe the weight RMSE values decreasing as training proceeds correctly, but after a point (typically around 15 training rounds for the buffer size that we used in our simulations), the weight RMSE values begin to increase again as the training is effectively “undone” by continuing to supply the circuit with new training instances that it cannot process correctly. As shown in the Supporting Information, periodically replenishing the supply of inactive buffered gates in the system can counteract this effect and extend the useful life of the learning circuit. The learning circuit design presented in this article uses the stochastic gradient descent weight update rule in conjunction with a linear transfer function, which allows the circuit to learn linear classification functions. However, this design could, in principle, be used to learn other functional forms by simply replacing the predictor subcircuit so that it computes the desired function of the provided training inputs. A major challenge is to implement nonlinear transfer functions such as the Heaviside (step) function, which is used in classical expositions of perceptron learning.51 Previous experimental work on artificial neural network implementations in DNA strand displacement networks16 implemented a traditional neural network using a Heaviside transfer function by using an extended toehold on a “sink” gate that serves as a threshold for the input species and by supplying an amount of fuel for signal amplification equal to the expected output level. This provides a physical limit on the level to which input signals could be amplified. However, this approach is not compatible with a reusable circuit for multiround training because after the first round of training an unknown quantity of that fuel may remain, making it difficult to restore the system to a state such that it will process the second round of training inputs correctly. To solve this problem, it may be necessary to resort to alternative schemes based on continual production and degradation of species, as occurs in transcriptional circuits, because in this context network motifs such as negative autoregulation52 can be used to limit amplified output signals to a particular level. Finally, it is well-known that, even for perceptron-like classifiers using nonlinear transfer functions, a single unit is capable of learning only linearly separable classification functions. To build molecular systems that can learn arbitrary functions, it would be necessary to connect a number of such units into networks. In this context, the network would need to be trained as a whole using backpropagation, which could be achieved by cascading46,47 several of the circuit motifs described here to produce the output and also cascading the feedback signals to implement the training. The development of such networks would go a long way toward realizing the potential of molecular computing for implementing nanoscale systems capable of highly sophisticated information processing on a par with that achieved by naturally occurring biological systems.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acssynbio.6b00009. Model and simulation code as well as results from additional simulations (PDF)



AUTHOR INFORMATION

Corresponding Author

*Phone: +1 505 2773351. Fax: +1 505 2775433. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This material is based upon work supported by the National Science Foundation under grants 1318833, 1422840, 1518861, and 1525553. M.R.L. thanks the New Mexico Cancer Nanoscience and Microsystems Training Center for additional support.



REFERENCES

(1) Zhang, D. Y., and Seelig, G. (2011) Dynamic DNA nanotechnology using strand-displacement reactions. Nat. Chem. 3, 103− 113. (2) Dirks, R. M., Bois, J. S., Schaeffer, J. M., Winfree, E., and Pierce, N. A. (2007) Thermodynamic analysis of interacting nucleic acid strands. SIAM Rev. 49, 65−88. (3) Schaeffer, J. M., Thachuk, C., and Winfree, E. (2015) Stochastic Simulation of the Kinetics of Multiple Interacting Nucleic Acid Strands. Proceedings of the 21st International Conference on DNA Computing and Molecular Programming 9211, 194−211. (4) Soloveichik, D., Seelig, G., and Winfree, E. (2010) DNA as a universal substrate for chemical kinetics. Proc. Natl. Acad. Sci. U. S. A. 107, 5393−5398. (5) Cardelli, L. (2009) Strand Algebras for DNA Computing. Proceedings of the 15th International Conference on DNA Computing and Molecular Programming 5877, 12−24.

3. METHODS 3.1. Construction of Chemical Reaction Models for Learning Circuit Simulations. We used the four-domain strand displacement encoding of abstract chemical reactions from Soloveichik4 to encode our learning circuit design as a K

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology (6) Cardelli, L. (2013) Two-domain DNA strand displacement. Math. Structures Comput. Sci. 23, 247−271. (7) Thachuk, C., Winfree, E., and Soloveichik, D. (2015) Leakless DNA Strand Displacement Systems. Proceedings of the 21st International Conference on DNA Computing and Molecular Programming 9211, 133−153. (8) Yurke, B., Turberfield, A. J., Mills, A. P., Jr., Simmel, F. C., and Neumann, J. L. (2000) A DNA-fuelled molecular machine made of DNA. Nature 406, 605−608. (9) Zhang, D. Y., and Winfree, E. (2009) Control of DNA strand displacement kinetics using toehold exchange. J. Am. Chem. Soc. 131, 17303−17314. (10) Srinivas, N., Ouldridge, T. E., Šulc, P., Schaeffer, J. M., Yurke, B., Louis, A. A., Doye, J. P. K., and Winfree, E. (2013) On the biophysics and kinetics of toehold-mediated DNA strand displacement. Nucleic Acids Res. 41, 10641−10658. (11) Seelig, G., Soloveichik, D., Zhang, D. Y., and Winfree, E. (2006) Enzyme-free nucleic acid logic circuits. Science 314, 1585−1588. (12) Qian, L., and Winfree, E. (2011) Scaling up digital circuit computation with DNA strand displacement cascades. Science 332, 1196−1201. (13) Zhang, D. Y., Turberfield, A. J., Yurke, B., and Winfree, E. (2007) Engineering entropy-driven reactions and networks catalyzed by DNA. Science 318, 1121−1125. (14) Chen, Y.-J., Dalchau, N., Srinivas, N., Phillips, A., Cardelli, L., Soloveichik, D., and Seelig, G. (2013) Programmable chemical controllers made from DNA. Nat. Nanotechnol. 8, 755−762. (15) Rudchenko, M., Taylor, S., Pallavi, P., Dechkovskaia, A., Khan, S., Butler, V. P., Jr., Rudchenko, S., and Stojanovic, M. N. (2013) Autonomous molecular cascades for evaluation of cell surfaces. Nat. Nanotechnol. 8, 580−586. (16) Qian, L., Winfree, E., and Bruck, J. (2011) Neural network computation with DNA strand displacement cascades. Nature 475, 368−372. (17) Alberts, B., Johnson, A., Lewis, J., Raff, M., Roberts, K., and Walter, P. (2008) Molecular Biology of the Cell, 5th ed., Garland Science. (18) Armus, H. L., Montgomery, A. R., and Jellison, J. L. (2006) Discrimination Learning in Paramecia (P. caudatum). Psychol. Rec. 56, 489−498. (19) Brembs, B., Lorenzetti, F. D., Reyes, F. D., Baxter, D. A., and Byrne, J. H. (2002) Operant Reward Learning in Aplysia: Neuronal Correlates and Mechanisms. Science 296, 1706−1709. (20) Fernando, C. T., Liekens, A. M. K., Bingle, L. E. H., Beck, C., Lenser, T., Stekel, D. J., and Rowe, J. E. (2009) Molecular circuits for associative learning in single-celled organisms. J. R. Soc., Interface 6, 463−469. (21) McGregor, S., Vasas, V., Husbands, P., and Fernando, C. (2012) Evolution of associative learning in chemical networks. PLoS Comput. Biol. 8, e1002739. (22) Scialdone, A., Mugford, S. T., Feike, D., Skeffington, A., Borrill, P., Graf, A., Smith, A. M., and Howard, M. (2013) Arabidopsis plants perform arithmetic division to prevent starvation at night. eLife 2, e00669. (23) Morens, D. M., and Fauci, A. S. (2013) Emerging Infectious Diseases: Threats to Human Health and Global Stability. PLoS Pathog. 9, e1003467. (24) Lakin, M. R., Youssef, S., Cardelli, L., and Phillips, A. (2012) Abstractions for DNA circuit design. J. R. Soc., Interface 9, 470−486. (25) Duda, R. O., Hart, P. E., and Stork, D. G. (2001) Pattern Classification, 2nd ed., John Wiley & Sons. (26) Lakin, M. R., and Stefanovic, D. (2015) Supervised learning in an adaptive DNA strand displacement circuit. Proceedings of the 21st International Conference on DNA Computing and Molecular Programming 9211, 154−167. (27) Cardelli, L. (2011) Strand Algebras for DNA Computing. Nat. Comput. 10, 407−428.

(28) Lakin, M. R., Youssef, S., Polo, F., Emmott, S., and Phillips, A. (2011) Visual DSD: a design and analysis tool for DNA strand displacement systems. Bioinformatics 27, 3211−3213. (29) Genot, A. J., Zhang, D. Y., Bath, J., and Turberfield, A. J. (2011) Remote toehold: a mechanism for flexible control of DNA hybridization kinetics. J. Am. Chem. Soc. 133, 2177−2182. (30) Qian, L., and Winfree, E. (2011) A simple DNA gate motif for synthesizing large-scale circuits. J. R. Soc., Interface 8, 1281−1297. (31) Zhang, D. Y., and Seelig, G. (2011) DNA-based fixed gain amplifiers and linear classifier circuits. Proceedings of the 16th International Conference on DNA Computing and Molecular Programming 6518, 176−186. (32) Yordanov, B., Kim, J., Petersen, R. L., Shudy, A., Kulkarni, V. V., and Phillips, A. (2014) Computational Design of Nucleic Acid Feedback Control Circuits. ACS Synth. Biol. 3, 600−616. (33) Oishi, K., and Klavins, E. (2011) Biomolecular implementation of linear I/O systems. IET Syst. Biol. 5, 252−260. (34) Lakin, M. R., Minnich, A., Lane, T., and Stefanovic, D. (2014) Design of a biochemical circuit motif for learning linear functions. J. R. Soc., Interface 11, 20140902. (35) Banda, P., Teuscher, C., and Lakin, M. R. (2013) Online learning in a chemical perceptron. Artif. Life 19, 195−219. (36) Banda, P., Teuscher, C., and Stefanovic, D. (2014) Training an asymmetric signal perceptron through reinforcement in an artificial chemistry. J. R. Soc., Interface 11, 20131100. (37) Banda, P., and Teuscher, C. (2014) Learning Two-input Linear and Nonlinear Analog Functions with a Simple Chemical System. Proceedings of the Unconventional Computation and Natural Computation 8553, 14−26. (38) Banda, P., and Teuscher, C. (2014) An Analog Chemical Circuit with Parallel-Accessible Delay Line for Learning Temporal Tasks. Artificial Life 14: International Conference on the Synthesis and Simulation of Living Systems, 482−489. (39) Kim, J., Hopfield, J. J., and Winfree, E. (2004) Neural network computation by in vitro transcriptional circuits. Adv. Neural Inf. Process. Syst. 17, 681−688. (40) Genot, A. J., Fujii, T., and Rondelez, Y. (2013) Scaling down DNA circuits with competitive neural networks. J. R. Soc., Interface 10, 20130212. (41) Chiu, T.-Y., Chiang, H.-J. K., Huang, R.-Y., Jiang, J.-H. R., and Fages, F. (2015) Synthesizing Configurable Biochemical Implementation of Linear Systems from Their Transfer Function Specifications. PLoS One 10, e0137442. (42) Stojanovic, M. N., Mitchell, T. E., and Stefanovic, D. (2002) Deoxyribozyme-based logic gates. J. Am. Chem. Soc. 124, 3555−3561. (43) Elbaz, J., Lioubashevski, O., Wang, F., Remacle, F., Levine, R. D., and Willner, I. (2010) DNA computing circuits using libraries of DNAzyme subunits. Nat. Nanotechnol. 5, 417−422. (44) Stojanovic, M. N., Stefanovic, D., and Rudchenko, S. (2014) Exercises in Molecular Computing. Acc. Chem. Res. 47, 1845−1852. (45) Brown, C. W., III, Lakin, M. R., Stefanovic, D., and Graves, S. W. (2014) Catalytic molecular logic devices by DNAzyme displacement. ChemBioChem 15, 950−954. (46) Brown, C. W., III, Lakin, M. R., Horwitz, E. K., Fanning, M. L., West, H. E., Stefanovic, D., and Graves, S. W. (2014) Signal propagation in multi-layer DNAzyme cascades using structured chimeric substrates. Angew. Chem., Int. Ed. 53, 7183−7187. (47) Lakin, M. R., Brown, C. W., III, Horwitz, E. K., Fanning, M. L., West, H. E., Stefanovic, D., and Graves, S. W. (2014) Biophysically inspired rational design of structured chimeric substrates for DNAzyme cascade engineering. PLoS One 9, e110986. (48) Wang, J., and Katz, E. (2010) Digital biosensors with built-in logic for biomedical applications−biosensors based on a biocomputing concept. Anal. Bioanal. Chem. 398, 1591−1603. (49) Katz, E., and Privman, V. (2010) Enzyme-based logic systems for information processing. Chem. Soc. Rev. 39, 1835−1857. (50) Privman, V., and Katz, E. (2015) Can bio-inspired information processing steps be realized as synthetic biochemical processes? Phys. Status Solidi A 212, 219−228. L

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX

Research Article

ACS Synthetic Biology (51) Minsky, M., and Papert, S. (1972) Perceptrons: An Introduction to Computational Geometry, 2nd ed., MIT Press, Cambridge, MA. (52) Alon, U. (2006) An Introduction to Systems Biology: Design Principles of Biological Circuits, Chapman & Hall/CRC. (53) The Mathworks, Inc. (2013) MATLAB, release 2013a, The Mathworks, Inc., Natick, MA.

M

DOI: 10.1021/acssynbio.6b00009 ACS Synth. Biol. XXXX, XXX, XXX−XXX