Chemical reaction network designs for asynchronous logic circuits

Chemical reaction networks (CRNs) are a versatile language for describing the dynamical behaviour of chemical kinetics, capable of modelling a variety of digital and analogue processes. While CRN designs for synchronous sequential logic circuits have been proposed and their implementation in DNA demonstrated, a physical realisation of these devices is difficult because of their reliance on a clock. Asynchronous sequential logic, on the other hand, does not require a clock, and instead relies on handshaking protocols to ensure the temporal ordering of different phases of the computation. This paper provides novel CRN designs for the construction of asynchronous logic, arithmetic and control flow elements based on a bi-molecular reaction motif with catalytic reactions and uniform reaction rates. We model and validate the designs for the deterministic and stochastic semantics using Microsoft’s GEC tool and the probabilistic model checker PRISM, demonstrating their ability to emulate the function of asynchronous components under low molecular count.


Introduction
Chemical Reaction Networks (CRNs) are traditionally used to capture the behaviour of inorganic and organic chemical reactions in a well-mixed solution. Recently, a paradigm shift in the scientific community has seen the use of CRNs extend to that of a high-level programming language for molecular computing devices (Cook et al. 2009), where the fundamental computational process differs from conventional digital electronics in that it involves transformation of input chemicals into output via reaction rules, as opposed to processing discrete signals (voltage bands) interpreted as Boolean values. Several digital and analogue circuits (Magnasco 1997;Soloveichik et al. 2008) have been designed in CRNs and their computational power studied (Soloveichik et al. 2010;. It has also been demonstrated in principle that any CRN can be physically realised in DNA (Soloveichik et al. 2010;Cardelli 2010;. CRNs are therefore particularly attractive as a programming language for use in nanotechnology and biomedical applications, where it is difficult to integrate traditional electronics. Chemical systems can store and process information in several ways. We focus on finite systems of molecules interacting in a well-mixed solution under mass-action kinetics and emulate Boolean circuits by encoding information through molecular concentrations reaching a particular threshold. The computation proceeds by transforming input species concentrations into outputs according to the reactions of a finite CRN. It is known that the computational power of CRNs is affected by the choice of the semantics, deterministic or stochastic. In particular, assuming a small probability of error, (finite) stochastic CRNs have been shown to be Turing universal (Soloveichik et al. 2008). The deterministic semantics interprets the reactions as a system of differential equations, which describe the evolution of the system as a vector of realvalued species concentrations over time . The stochastic semantics, on the other hand, views the state of the system as a vector of (non-negative) integer molecular counts and state transitions as a reaction which has a non-zero probability of occurring (Cook et al. 2009). The stochastic evolution of the system over time is obtained as a solution of the Chemical Master Equation (CME) (Kampen 1992). It is well known that the deterministic semantics is not accurate for small populations. While the stochastic semantics is exact, it is infeasible for large molecular counts. One scalable alternative is the Linear Noise Approximation, which is a real-valued approximation of the CME (Cardelli et al. 2015). The correctness of the behaviour of a circuit described by a finite CRN can be analysed by inspecting its stochastic and deterministic evolution over time. In addition, techniques such as model checking can be employed to analyse the temporal ordering of events.
While CRN designs for synchronous sequential logic circuits have been proposed, to mention (Magnasco 1997;Soloveichik et al. 2010Soloveichik et al. , 2008, a physical realisation of these devices is challenging because of their reliance on a clock to synchronise events in order to ensure the correct temporal order of the phases of the computation. Clocks are difficult to make, since they arise from unique conditions of chemical concentrations and kinetic constants, and must control a large number of events. In electronics, an alternative circuit design technology is asynchronous sequential logic (Spars and Furber 2002;Myers 2004), which instead of a clock relies on handshaking protocols to synchronise events. Asynchronous circuits are widely used for low-power microprocessor designs, e.g., by ARM, though require a larger circuit area. The key component is the Muller C-element, which is used to synchronise multiple independent processes in a manner insensitive to the delays on wires and individual components. To ensure Turing completeness of asynchronous circuits, we also require an isochronous fork in addition to the Muller C-element. An isochronous fork is a component which produces a fan-out of signals that reach the target at virtually the same time. This assumption is difficult to achieve in conventional electronics, because of the need to make the wires the same length, but is straightforward in chemical kinetics because of the well-mixed assumption.
This paper provides novel CRN designs for the construction of an asynchronous computing device based on a bi-molecular reaction motif inspired by the Approximate Majority network (Angluin et al. 2008;Cardelli and Csikász-Nagy 2012). The motif employs catalytic reactions to achieve bistable switching of molecular concentrations, which emulates high and low voltage signals in digital electronics. All components are produced with simple reactions and uniform reaction rates, where we assume a well-mixed solution under mass action kinetics, and are independent of a universal clock. Moreover, any design provided in this paper could in principle be realised as a DNA strand displacement device (Cardelli 2010). We work with the dual-rail design methodology and employ a variant of the diagrammatic language of Cardelli (2014) to represent the designs at the high level. Starting from the Muller C-element, we design the main components of a complete asynchronous computing device in terms of CRNs in a principled way, including logic gates, control flow and basic arithmetic, as well as more complex structures such as queues. We validate the designs by exploring their time evolution for all possible combinations of inputs using Microsoft's Visual GEC tool, 1 with the latter also approximated using an experimental implementation of the Linear Noise Approximation (LNA) of Cardelli et al. (2015) provided by Visual GEC that offers better scalability. We use the LNA to highlight a flaw with a key design component. Further, we demonstrate the correct behaviour of the circuits against temporal logic specifications with the probabilistic model checker PRISM 2 (Kwiatkowska et al. 2011). Our designs constitute the first feasible implementation of asynchronous computational components as CRNs, and are relevant for a multitude of applications in synthetic biology and biosensing.
This paper is an extended version of the conference paper Cardelli et al. (2016).

Related work
The computational power of Chemical Reaction Networks, viewed as a programming language for engineering biochemical systems, has been studied by a number of authors, to mention Cook et al. (2009) and . There are a number of ways in which chemical systems can encode and process information. This includes simulating Boolean circuits, where information is encoded in binary form using high and low concentrations similarly to this paper, e.g. Magnasco (1997), Soloveichik et al. (2010) and Soloveichik et al. (2008), as well as geometric arrangements, for example self-assembly (Rothemund et al. 2004) and molecular walkers (Dannenberg et al. 2015) not considered here. Researchers have also investigated the power of CRNs to model distributed algorithms (Angluin et al. 2008).
Regarding synchronous logic circuits, much of the work to date considered abstract CRN schemes. One exception is Silva and McClenaghan (2004), where a system of actual chemical reactions is given, together with a precise molecular implementation for gates complete with a thermodynamic analysis of how the system would evolve, though only for simple gate designs. In Cook et al. (2009) we see the construction and composition of simple logic gates based upon catalytic reactions, but they do not mention control flow or systematic component design in a dual rail setting. In Senum and Riedel (2011) the authors propose CRNs for an inverter, an incrementer, a decrementer and a copier; their designs are based on two rate constants, ''fast'' and ''slow'', and thus are not rate-independent, in contrast to the designs presented here.
CRNs can also be viewed as computing functions over reals or Booleans. A single CRN computes a function over a finite domain, which is analogous to Boolean circuits in the sense that any given circuit computes only on inputs of a particular size (Soloveichik et al. 2008). An implementation of dual-rail logic gates that are rate-independent is given in Chen et al. (2014). In contrast, our designs are composable and capable of performing non-trivial computation.
Since the behaviour of CRNs is asynchronous, a fact evident through their equivalence with Petri net models (Cook et al. 2009), the main difficulty with programming them is the need to control the order of reactions. In Cook et al. (2009) it is suggested that this ''uncontrollability'' can be handled by changing rate constants, an idea followed up in Napp and Adams (2013), where CRN designs for basic arithmetic are given based on two rate constants, ''fast'' and ''slow''. Our designs, on the other hand, exploit the asynchrony of the underlying CRN model and work with uniform rates.
Designs for the Muller C-element, though not the remaining components of an asynchronous device, have been constructed from genetic logic gates (Nguyen et al. 2010) and a genetic toggle switch (Nguyen et al. 2007), but we are not aware of any other nanoscale designs for asynchronous circuits. Soloveichik et al. (2010) shows that any CRN, including those presented in this paper, can theoretically be implemented as a DNA Strand Displacement device. These devices have been demonstrated in the lab , and thus provide an indication of experimental feasibility of our designs.

Chemical reaction networks
A Chemical Reaction Network (CRN) C ¼ ðK; RÞ is a pair of finite sets, where K is a set of chemical species and and R is a set of reactions. jKj denotes the size of the set of species. Reactions in R describe how species interact. Formally, a reaction s 2 R is a triple s ¼ ðr s ; p s ; k s Þ, where r s 2 N jKj is the vector of molecular counts of the reactants, p s 2 N jKj is the vector of molecular counts of the products and k s 2 R [ 0 are the coefficient associated to the rate of the reaction. We assume ordering of species within vectors is alphabetical. Given a reaction s 1 ¼ ð½1; 1; 0 T ; ½0; 0; 2 T ; k 1 Þ, where Á T is the transpose of a vector, we often refer to it as s 1 : A, B and C are generic species.
In this paper we are only concerned with uni-molecular reactions, i.e. those which have only one reactant, and bimolecular, i.e. those with two reactants. The ''reversible reaction'' notation A þ B 2C is a shorthand for the two where k 1 and k 2 are not necessarily equal. We assume that the system is well stirred, that is, the probability of the next reaction occurring between two molecules is independent of the location of those molecules, at fixed volume V and temperature; under these assumptions a configuration or state x 2 N jKj of the CRN is given by the number of molecules of each species. Given a configuration x we define z ¼ x N , where N ¼ V Á N A is the volumetric factor, V is the volume and N A Avogadro's number. We write xðk i Þ for the number of molecules of k i in the configuration x and zðk i Þ ¼ xðk i Þ N to denote the concentration of k i in the same configuration.
We will sometimes distinguish between CRNs with different initial configurations, and to this end define a chemical reaction system (CRS) as a tuple S ¼ ðK; R; x 0 Þ where ðK; RÞ is a CRN and x 0 2 N jKj represents its initial configuration, and we sometimes use the terms CRN and CRS interchangeably.
Diagrammatic CRN notation To better visualize a CRN C ¼ ðK; RÞ as a circuit, we employ a directed multiedge graph ðK; EÞ based upon a fragment of the diagrammatic notation for influence graphs Cardelli (2014). K, the nodes, represent the species of the CRN C and edges E are derived from reactions R as follows. A reaction is represented as a directed multi-edge with sets of species as source and target. Each edge is either a pointed arrow (") or a rounded arrow ( ), with the source represented by the flat edge and the target represented by the arrow head. A reaction that produces a species as a product is connected to it by a directed edge.
All reactions within our diagrams are catalytic and are bi-molecular reactions of the form X þ Y ! X þ Z, meaning that Y is transformed to Z and X is a catalyst, that is, X influences the transformation of Y to Z. The edges represent that a source species is catalytic to a target reaction.
Example 1 (CRN Diagrams Example) We illustrate the flexibility of the diagrammatic notation with three CRN examples. Figure 1a shows the CRN with the single reaction C þ A ! C þ B. The ball ( ) indicates that C is a catalyst to the reaction A ! B represented by the arrow ("). A species can act as both a reactant or catalyst in the same reaction, and similarly for a product. In Fig. 1b we depict the CRN with two reactions fA þ A ! A þ B; Bþ B ! B þ Ag, in which both A and B catalyse themselves. Figure 1c depicts the CRN with reaction set fA þ B ! B þ B; C þ B ! B þ Bg, with species B catalysing multiple reactions, which is represented by a multiheaded edge with multiple s. This CRN transforms species A and C into species B.
Dual-rail representation We represent a Boolean circuit with inputs I and outputs O, denoted B(I, O), as follows. Firstly, a Boolean variable b ¼ f0; 1g could be encoded in a single species X, where 0 would be encoded as E½jXj ¼ 0 and 1 as E½jXj ! M, where E[|X|] denotes an expectation of the number of molecules of X and M is a molecular population threshold.
The CRN computes by transforming an input concentration into an output concentration, which reaches the appropriate level upon convergence. However, since absence of molecules cannot be measured, we employ dual-rail methodology and represent every Boolean variable with two species, denoted X hi ; X lo . Just like we cannot represent both 0 and 1 on an electrical wire, we restrict our CRNs such that either E½jX hi j ! M or E½jX lo j ! M, but not both, can be present when a CRN has stabilised and no further reactions occur. We consider a high concentration output as correct if E½jX hi j [ 0:8maxðX hi Þ À 1SDðX hi Þ, where 0.8 is a threshold normalised between the values [0, 1], max() is a function that returns the maximum molecular concentration of a species within the CRN and 1SD computes 1 standard deviation from the mean concentration of E½X hi . 1SD returns 0 under the deterministic semantics. Similarly, we consider a low concentration as correct if E½jX lo j\0:2 Ã maxðX lo Þ þ 1SDðX lo Þ.
For simplicity, we apply the dual rail methodology only to the variables in the input and output sets I and O. The circuits may contain additional variables, which will be considered internal and assumed not to catalyse with any species outside of the CRN circuit. We will encode these with single species and use the naming convention of referring to these internal species as k; k 1 ; Á Á Á ; k i . When composing two circuits B 1 ðI 1 ; O 1 Þ and B 2 ðI 2 ; O 2 Þ in series, we define their composition as a circuit BðI 1 ; O 2 Þ, in which all variables in O 1 [ I 2 have been made internal.
Example 2 (Dual-Rail CRN 'Motif') We introduce a simple two reaction CRN which forms a 'motif' common to all our CRN circuit diagrams. The CRN is given by the Fig. 2a, where the input set I contains X lo and X hi and the output set contains Y lo ; Y hi . In this CRN the input X lo or X hi influences the reaction Y lo Y hi to produce as output the same Boolean value. We also include a diagrammatic CRN showing the composition of two such motifs in series. Here the inputs are X hi ; X lo and outputs Y hi ; Y lo .
Example 3 (Dual-rail CRN with Internal (k) Species) The CRN shown in Fig. 3 represents a circuit B(I, O) with internal species k which does not belong to the set I [ O. Therefore dual rail methodology is not used for k as it is not catalytic to any species outside of this CRN circuit. I comprises X lo ; This CRN converts Y hi and Y lo to k, assuming a non-zero initial concentration of molecules, unless there is either X lo or X hi present, in which case a deadlock occurs. We use the term conversion to mean the occurrence of a reaction where there are nonzero molecular counts of reactants present.
Deterministic semantics Let C ¼ ðK; RÞ be a CRN. The net change associated to a reaction s 2 R is defined by t s ¼ p s À r s . The deterministic semantics models the concentration of the species in K over time as a set of autonomous polynomial first order differential equations (ODEs): which species A and B are both reactants, products and catalysts. c CRN with reactions fA þ B ! B þ B; C þ B ! B þ Bg, which demonstrates that B can be catalytic to multiple reactions where function U : R ! 0 ! R K describes the behaviour of the system assuming a continuous state-space semantics, and therefore UðtÞ 2 R jKj is the vector of the species concentrations at time t and F is simply the derivative of U with respect to time. Assuming t 0 ¼ 0, the initial condition is Uð0Þ ¼ ½x 0 , where x 0 is the initial configuration (vector of concentrations of molecules) of the CRN. It is well known that the deterministic semantics may be imprecise for low molecular counts, but is accurate in the limit for high populations Kampen (1992). However, the deterministic semantics produces the same proportion of molecules, regardless of total concentration.

Stochastic semantics
The stochastic semantics is represented through a continuous-time Markov chain (CTMC), whose transient evolution can be given via the Chemical Master Equation (CME) (Kampen 1992). Let C ¼ ðK; RÞ be a CRN. The propensity rate a s of a reaction s is a function of the current configuration of the system x such that a s ðxÞdt is the probability that a reaction event occurs in the next infinitesimal interval dt. We assume mass action kinetics, therefore a s ðxÞ ¼ k s where r i;s is the i-th component of the vector r s and jr s j ¼ P jKj i¼1 r i;s Anderson and Kurtz (2011). We define a time-homogeneous CTMC ðX C ðtÞ; t 2 R ! 0 Þ with state space Q N jKj as follows. Given x 0 2 Q, where x 0 is the initial configuration of the system, then describes the stochastic evolution of the molecular populations of each species in C at time t. For x 2 Q, we define P ðtÞ ðxÞ ¼ PðXðtÞ ¼ xjXð0Þ ¼ x 0 Þ, where x 0 is the initial configuration. The CME describes the time evolution of X as: The solution of the CME is computed through numerical simulation or discretisation techniques such as uniformization , and is generally feasible only for small populations. The CTMC is often represented as a Q Â Q rates matrix, which can be viewed as a state transition graph and subjected to model checking against temporal logic properties (Kwiatkowska et al. 2011).
Linear noise approximation The LNA approximates the CTMC as a continuous-state Gaussian process, given in the form of a set of ODEs that describe the time evolution of expectation and variance of the species. The error of approximation is dependent upon the volumetric factor N, the structure and the rates of the CRN. Given a CRN C ¼ ðK; RÞ with initial configuration x 0 2 N jKj and in a system of volume size N, we define the stochastic process where U is the deterministic process given in Eq. 1, and Z is a zero-mean Gaussian process (since we assume the initial condition is a fixed value), and with covariance C[Z(t)] described by the solution of the following ODEs with initial condition C½Zð0Þ ¼ 0: Fig. 2 a The CRN 'motif', in which the input X lo or X hi influences the reaction Y lo Y hi to produce as output the same Boolean value. b CRN obtained by composing two 'motif's depicted in a in series. The inputs are now X hi ; X lo , outputs Y hi ; Y lo , and the remaining species have been made internal through renaming with k.

Fig. 3
Example CRN with internal k species, which converts outputs Y hi and Y lo to k, assuming a non-zero initial concentration of molecules, unless there is input X lo or X hi present Chemical reaction network designs for asynchronous logic circuits 113 dC½ZðtÞ dt ¼F 0 ðUðtÞ; C½ZðtÞÞ where J F ðUðtÞÞ is the Jacobian of FðUðtÞÞ, J T F ðUðtÞÞ its transpose, and Expected value and covariance matrix of Y(t) are completely characterized by UðtÞ and C[Z(t)] since E½YðtÞ ¼ NUðtÞ and C½YðtÞ ¼ ffiffiffiffi The LNA requires solving a number of ODEs quadratic in the number of species (Cardelli et al. 2015) and is therefore a scalable alternative to the solution of the CME. In contrast to the deterministic semantics, which considers average concentrations, the LNA does not compromise stochasticity.
Tool support A number of software tools exist for examining the behaviour of CRNs. We employ Microsoft's Visual GEC, which provides a programming language, LBS, for designing and simulating a given CRN under the deterministic or stochastic semantics, including also the LNA approximation of the stochastic semantics. The tool is capable of producing plots of expected or average species concentrations over time. This functionality is used extensively within this paper to validate our circuit designs.
In addition, Visual GEC exports models to the probabilistic model checker PRISM (Kwiatkowska et al. 2011), which then enables verification of the induced continuous-time Markov chain against temporal logic properties. We use PRISM to verify the correctness of the temporal ordering of events occurring as the CRN circuit executes.

Example 4 (CRN Validation Under Different Semantics)
We show the operation of our dual-rail CRN 'motif' given in Example 2 under the deterministic and stochastic semantics. The input configuration is jX hi j ¼ 10 molecules and output jY lo j ¼ 10 molecules. Figure 4 demonstrates that, after 0.4 s, the CRN stabilises reaching the concentrations of jY hi j ¼ 10 and jY lo j ¼ 0 as desired. With regards to simulations provided, concentrations (given in nanomolars) are directly correlated to concentrations as we assume the volumetric factor N is fixed.

Asynchronous hardware
Asynchronous computation is a model of computation that relies on transitions via local input signals rather than transitions via a global clock. Asynchronous computation (Spars and Furber 2002), just like its synchronous counterpart, is Turing complete (Manohar and Martin 1996), meaning that any bounded-tape Turing machine can be implemented with an asynchronous circuit, providing that the implementation of that circuit has isochronous forks. An isochronous fork is the propagation of a signal from a Fig. 4 Simulation of the motif given in Example 2 under different semantics. The input is initially X hi and the output is Y lo , both starting at a concentration of 10 molecules. In a, b, c we respectively show the deterministic solution, stochastic simulation and LNA plot with variance, where Y lo is seen in blue and Y hi is seen in yellow. We can observe that in all cases Y hi is present after 0.4 s. In d we show the state transitions of the induced CTMC that correspond to the output switching from Y lo to Y hi . The computation reaches the correct output state and stabilises, with no more transitions enabled. (Color figure online) single source to multiple receivers, with the important constraint that the signal must reach the receivers at precisely the same time. In classical digital circuitry this corresponds to the propagation of a signal down wires of exactly the same length from one component to another. In CRNs this could be seen as two species reaching a threshold N at precisely the same time.
Asynchronous circuits, which are designed based upon the theoretical principles of asynchronous computation, are widely used for low-power microprocessor designs, e.g., by ARM, and are increasing in popularity with the increase in distributed computing (Myers 2004). Asynchronous designs offer a number of advantages, the main one being correctness independent of timing, although they require a greater overhead in terms of silicon area.
We illustrate the principles of asynchronous circuit design by describing its key component called the Muller C-element and showing how it is used to construct a pipeline that propagates signals.
Muller C-element The cornerstone of asynchronous computation is the Muller C-element. A Muller C-element has two Boolean inputs, X and Y, and one output Z. By definition these inputs can either be low or high (represented by 0 or 1). When both inputs are low the output is low. Similarly, when both inputs are high the output is high. The variation from a classical logic gate, however, is that if the inputs are high, or low, and one of them changes, it 'remembers' the last state. In other words, it retains the last 0 or 1 state. This is summarised in Fig. 5a. An important property of the C-element is that it allows an observer to conclude on seeing output change from 0 to 1 that both inputs are now 1, and similarly for input change from 1 to 0.
The table specification indicates that asynchronous circuits exhibit concurrency and causality, and hence their specifications need to reflect these characteristics. A common way is as a timing diagram, seen in Fig. 5d for the C-element, which represents a set of signals and their interactions over time. Each row of a timing diagram represents one signal and how it switches from low to high over time. If a signal displays a change before another signal on another line then this signal must precede the other. An arrow represents that one signal change triggers the change of another. In the C-element diagram, note that X and Y have to precede Z, both in the transition to 1 and down to zero. However, there is no causal dependency between X and Y.
Asynchronous diagrams, and in particular the C-element, are accurately described using Petri nets (Spars and Furber 2002, p. 86) or process algebras (Wang and Kwiatkowska 2007). We present a (1-bounded) Petri net Fig. 5 Four specifications of a C-element with inputs X, Y and an output Z. a Conventional logic table, where 'unchanged' means that the state of the output is the last stable configuration of 1 or 0. b State graph, in which each state denotes a configuration and a transition is caused by the presence of a signal, where * indicates that the signal is excited. c 1-bounded Petri net specifying the Muller C-element. d Timing diagram for the C-element Chemical reaction network designs for asynchronous logic circuits 115 for the C-element in Fig. 5c, in which transitions are interpreted as signal transitions and places and arcs capture the causal relations between the signal transitions. Following the usual convention, the Petri net is drawn in simpler form where most places have been omitted. We can observe that both tokens are needed in order to excite the transitions that cause the event Z hi , which in turn will require the events X hi and Y hi to be triggered. The same is true for Z lo . When considering circuit synthesis, we typically employ a state graph specification, which can be obtained from the Petri net representation (Myers 2004) and is given for the C-element in Fig. 5b. The values in each state correspond to the values of inputs X, Y and output Z, respectively. A * symbol indicates that the corresponding variable is excited by the outgoing transition (and will be changed in the following state). Observe how we can only transition to a state 1Ã1Ã1 from a state 110Ã requiring X as 1 and Y as 1. Because this is derived from a 1-bounded Petri net, we can assume that the transitions 0Ã10 ! 110Ã and 10Ã0 ! 110Ã do not conflict.
Muller C-pipeline In order to replace the need for a global clock, asynchronous computation relies on 'local cooperation' in the form of handshaking protocols. These protocols exchange completion signals in order to establish when a computation has terminated. These handshaking protocols rely heavily on the C-element described above.
The Muller pipeline, shown in Fig. 6, is constructed by the composition of Muller C-elements (depicted by the gate symbol labelled with C) and classical NOT-gates, which receive and send data to/from the environment (Left, Right in the figure). Its function is to propagate a high and low signal along the pipeline, emulating the 'wave' of high and low signals of a classical synchronous clock. Initially, all C-elements are set to a value of 0. The ith C-element C[i] will propagate a 1 from its predecessor, C½i À 1, only if its successor, C½i þ 1, is 0. Similarly, it will propagate a 0 from its predecessor only if its successor is 1. Eventually the first request initialised on the left hand side of our pipeline is propagated to the final request on the right.
The protocol enacted upon this pipeline uses request and acknowledge rails that can be set to high or low. The Muller pipeline implements a basic four phase protocol, which is as follows. Firstly, the sender sends data and sets request to high, viewed in Fig. 6 as the signal Req being high. The receiver then records this data and sets acknowledge to high (Ack). Then the sender responds by setting request to low (Req), and finally the receiver acknowledges this by setting acknowledgement to low (Ack). If at any point a handshake along the pipeline is slower than another, the pipeline will behave like a FIFO queue with data preserved. Herein lies the important purpose of the pipeline: it allows for the delay-insensitive transfer of information from one place to another. In combination with a latch we can create the propagation of information across latches using the pipeline as a control structure.
The construction of data storage and control structures such as queues and adders is similar to the Muller pipeline.
Timing properties of asynchronous circuits Asynchronous circuits can be classified as being self-timed, speed independent or delay-insensitive, depending upon delay assumptions that are made. Assume a circuit is composed of gates and wires. A self-timed circuit operates correctly if both gates and wires experience measurable and fixed delays. A speed independent circuit operates correctly if gates exhibit some unknown time delay within gates but exhibits no time delay on wires. A delay-insensitive circuit operates correctly if there is both unknown delay within gates but also unknown delay within wires. The set of delay insensitive circuits is small, essentially those built from the Muller C-element and NOT gates, and so a broader class of quasi-delay insensitive circuits are identified. Quasi-delay insensitive circuits, which can be composed of purely C-elements, NOT-gates and forks, are Turing complete Manohar and Martin (1996). They are not possible to achieve without an isochronous fork.

Asynchronous circuit designs as chemical reaction networks
In this section we present our dual-rail designs for an asynchronous computing device in CRNs. The key component is the Muller C-element, whose design is inspired by the well known Approximate Majority (AM) CRN (Angluin et al. 2008). We begin by providing a detailed justification for our C-element design, and then describe the remaining simple components, including latches, logic gates and control flow. Finally, we present complex circuits such as the pipeline, queue and adder.
To justify the designs, we demonstrate that each component we design exhibits correct behaviour. Considering as an example the C-element, this amounts to working with an informal specification of the C-element in terms of high/ Fig. 6 Muller pipeline. Signals are propagated from left to right using request and acknowledge signals. The pipeline effectively queues data, only allowing a transition to occur when a further signal has been acknowledged low signals as in Fig. 5, and then showing that our (continuous) CRN empirically satisfies that specification according to appropriate thresholding for high/low signals, as normally done in electronics for transistor logic. To this end, we explore the time evolution of the components under the deterministic and stochastic semantics, including also LNA for scalability. We additionally employ probabilistic model checking with PRISM, where temporal logic is used to express the temporal ordering of events. We remark that, although we validate the components for all possible input configurations, this does not amount to full verification of the correctness of the designs. We discuss the challenges of achieving full verification in Sect. 6.

Muller C-element as a CRN
The C-element design is based on the AM CRN, which computes the majority of two finite populations by converting the minority population into the majority population, so that a single population emerges as output. It uses a third 'undecided' state of the population, from where catalysis can drive the individuals into either of the final states. Interestingly, since approximate majority cannot be exactly computed by a bi-molecular CRN with less than 4 reactions (Mertzios et al. 2014), below we present the bimolecular AM CRN with exactly four reactions: where X, Y are both the input and output species and k is the aforementioned catalytic driver. The intuition behind this reaction network is that we have two competing initial populations of X and Y, both of which try to eliminate the other by transforming their counterpart into the intermediary k. If k then interacts with X, it transforms itself into X, else into Y. Presented below in Fig. 7a is the same AM CRN in our diagrammatic notation.
Deconstructed, if we consider the left-hand side of the diagram, the reaction X ! k is catalysed by the species Y and so yields the bi-molecular reaction Y þ X ! Y þ k. Similarly, X catalyses the reaction k ! X in the other direction and yields the reaction X þ k ! X þ X. On the right-hand side of the diagram, k is again catalysed by X and Y to produce Y, yielding the other two reactions from the AM CRN. The CRN in Fig. 7b, with inputs X and Y and outputs W, Z, is similar to the AM CRN, except that we produce new arbitrary outputs W, Z instead of producing greater quantities of X and Y.
Since we wish to apply the dual-rail methodology, we represent each signal as a pair of species and encode the value 1 (0) as a molecular population of at least M for some population threshold M (population 0). We thus present in Fig. 8a a dual-rail CRN which computes approximate majority with four inputs X hi ; X lo ; Y lo ; Y hi and outputs Z hi ; Z lo . Like before, our input catalyses an intermediary species k, but this is split over two reactions. In addition, Z lo and Z hi catalyse with reactions Z lo k and Z hi k. This has the effect that, if there are larger numbers of either Z lo or Z hi , the network enlarges its majority by converting the other into k. Two rounded arrows ( ) over a reaction (") indicates that either species can act as a catalyst to that reaction. We demonstrate this AM effect in Fig. 8b, where, given the initial configuration of inputs X lo ; Y lo of 10 molecules and output Z hi , we observe under deterministic semantics that, because jX lo j; jY lo j [ jX hi j; jY hi j, an output of Z lo where jZ lo j ¼ 10 molecules and jZ hi j ¼ 0 is produced after 0.5 s.
We now need to justify the correctness of the design against the C-element specification in Fig. 8a. Given a starting configuration of input X, Y both 1, and any starting output, the C-element should eventually output 1. Similarly, if X, Y are both 0, on any initial output our final output should be 0. For any other configuration of X, Y, the output signal should remain the same. Thus, given an input X hi ; Y hi , and crucially any starting output configuration Z, we wish to reach a state where Z hi has at least M molecules where M is a population threshold, and Z lo has 0 molecules. Similarly, for X lo ; Y lo and any Z we should to see a presence of Z lo after some time t. For all other configurations of the inputs we wish the output species to remain the same.
When validated against this informal specification with a starting configuration of X lo ; Y hi at 10 molecules and Z hi at 10 molecules, our CRN unfortunately fails, seen in Fig. 9. More specifically, we observe that species Z hi is at a concentration of 6 molecules and k has a concentration of 4 molecules after 0.3 s. This simulation was produced using LNA, see Sect. 3.1, which outputs standard deviation of the mean concentrations of species, seen in the shaded regions. This output configuration is incorrect since Z hi is below the required threshold, see Sect. 3.1, of 8 molecules to represent an output of Z ¼ 1. Fig. 7 Two diagrammatic CRNs which are capable of computing Approximate Majority Angluin et al. (2008). In a we present the original in which the inputs X, Y, dependent on which species has the majority, influence outputs X, Y. b Shows a similar AM circuit, but now the input species are catalysed to arbitrary output species W, Z We present an amended CRN that resolves this issue in Fig. 11a. This CRN is composed of two approximate majority circuits connected to each other, with the outcome of the first AM amplifying the outcome of the second. As we can see from the plot in Fig. 9, we need to amplify the Z hi species and suppress the k species. The second AM corrects exactly this issue. Figure 11b is a simulation with inputs X lo ; Y hi at 10 molecules and output Z hi at 10 molecules. Here we can see that all three species are now at 10 molecules throughout the duration of the simulation.
To strengthen the validation of the final C-element design, we provide two further plots for selected initial configurations. We include in Fig. 11c a deterministic plot with starting configuration of input X lo ; Y lo at 10 molecules and output Z hi at 10 molecules. After 1 s the system converges to output Z lo at 10 molecules. Figure 11d shows an LNA simulation with starting configuration of input X hi ; Y hi at 10 molecules and output Z lo at 10 molecules. After 1 s we reach output Z hi with probability % 1. Both these simulations show that our output changes if both inputs change. From Fig. 11b we observe that the output does not change if both inputs are different.
An issue to address is the use of dual-rail systems in which our chemical output is precisely a value of 0 molecules or K molecules, where K is the largest number of molecules achievable by a population in the system. In reality, a species may not reach its maximum population and, due to variance, we may have a situation, as seen in Fig. 9, where one species has moderate probability of being higher than the output species we wish to present. Fortunately, we can use our approximate majority circuit to boost species. Figure 10 shows an example where species X hi and Y hi are weak, but the output is boosted by the C-element back to the maximum output of 10 molecules. The gates are also reusable: specifically, in Fig. 11d we can observe the inherent reusability of the C-element because the output reacts to the change in input.

Latch design
A latch is a device used in electronics to store a logical 0 or 1, which therefore needs to have at least two stable states that are cycled between. Latches are used in asynchronous computing both for storage and for synchronisation purposes. When an input of 1 is received a latch will ideally display an output of 1, and likewise for an input of 0. We present two latch designs in Fig. 12, each intended to interface in a specific way when used within a larger Fig. 8 A dual-rail Approximate Majority CRN. a Circuit diagram. b Deterministic simulation of the CRN in a. Given inputs X lo ; Y lo at 10 molecules and output Z hi , we observe that, because jX lo j; jY lo j [ jX hi j; jY hi j, an output of Z lo where jZ lo j ¼ 10 molecules (seen in blue) and jZ hi j ¼ 0 (seen in red) is expected to be produced after 0.5 s. (Color figure online) system. The first simple latch, shown in Fig. 12a(i), has input X lo ; X hi and output species Y lo ; Y hi , the intuition being that either X lo catalyses Y hi to Y lo or X hi catalyses Y lo to Y hi . There are also two additional reactions that catalyse Y lo to itself and Y hi to itself, creating a feedback loop. These additional reactions ensure that, if there is a drop in the molecular concentrations of input species, the latch retains its state. For some larger systems we may need the output state of a latch to be neither Y lo nor Y hi , to signify that no value is stored within the latch, known as a neutral state in electronics. A reset wire, to reset a latch to neutral state, is also commonly used in circuits. To this end, the latch in Fig. 12a(ii) has an input R hi ; R lo used to reset the latch to a central state k, as well as the standard inputs X hi ; X lo and outputs Y hi ; Y lo . The advantage of this central state, k, is that the latch can be in a state where neither Y hi nor Y lo are present, which is useful if these reactions are catalytic to any other component, in which case they will not be triggered directly. A comparison of the behaviour of the two latches is displayed in Fig. 12b, c. With the same initial conditions X hi and output Y lo at 10 molecules, the latch in Fig. 12a(ii) produces an output Y hi at 10 molecules in 0.2 s. We contrast this with the latch in Fig. 12c, which outputs Y hi at 10 molecules in the slower time of 0.5 s.

AM as a control flow element
An arbiter is used to decide an output signal based on which signal arrived first or if one signal is dominant over another. They are used in error correction where a signal may have degraded. Essentially, an arbiter computes the well known function maxðjX 1 j; jX 2 jÞ for two inputs. In terms of CRNs, this can be seen as either one species arriving before another or having higher molecular concentration. Since the AM circuit computes the maxðjX 1 j; jX 2 jÞ function, as one population is biased over another depending upon which has the majority, it therefore serves as an appropriate candidate for an arbiter. The Fig. 9 LNA plot for the candidate CRN for the dual-rail Muller C-element design in Fig. 8. With a starting configuration of inputs X lo ; Y hi at a concentration of 10 molecules and output Z hi at a concentration of 10 molecules, we can observe that after 0.3 s Z hi (seen in red) is at a concentration of 6 molecules and k is at a concentration of 4 molecules (seen in green). The shaded regions represent standard deviation. As we can see, with a non-zero probability we cannot distinguish the signals. (Color figure online) proposed arbiter design, seen in Fig. 13a, is the same as our AM CRN presented in Fig. 8a, except that there are two inputs, X hi and X lo , and two outputs, Y hi and Y lo , instead of four inputs. This works as desired since the output Y hi ; Y lo begins to be converted from k as soon as either of X hi ; X lo arrives, therefore automatically biasing whichever species is present first. The ability for approximate majority to reach a consensus means that this circuit can deal with stochastic fluctuations in input. Although, within electronic circuits, an arbiter outputs which signal arrived first, we assume that this information is revealed through the promotion of an output species linked to an input species.
In Fig. 13b we demonstrate the operation of this arbiter by LNA simulation on inputs of X hi at a concentration of 5 molecules and X lo with a concentration of 0 molecules. After 0.4 s we see Y hi (in red) at a concentration of 10 molecules, representing the majority.

Other control flow circuits
Control flow is used to mediate or propagate the flow of information throughout the computing device. In digital circuitry forks and joins, both control flow elements, are naturally implemented using wires. Unfortunately, there is no natural fork or join within CRNs and consequently we present designs for them. The fork, shown Fig. 14a, is used to split signals. It has two input species X hi ; X lo and four output species, Y hi ð1Þ; Y hi ð2Þ to represent the splitting of X hi and Y lo ð1Þ; Y lo ð2Þ to represent the splitting of X lo . The join, see Fig. 14b, joins two input signals to create one output signal. There are 4 inputs X hi ð1Þ; X hi ð2Þ with output Y hi to represent the merging of the two input signals and X lo ð1Þ; X lo ð2Þ to merge to an output Y lo .

Dual rail asynchronous logic gate designs
Although gate designs for Boolean operators have been proposed in CRNs (Soloveichik et al. 2008), we present Fig. 10 LNA simulation demonstrating two weak input signals, X hi and Y hi at 6 and 4 molecules respectively, boosted by the C-element with an output Z hi at 10 molecules Fig. 11 Simulation of the final dual-rail Muller C-element design on selected input configurations. a Dual-rail AM circuit which is our final C-element design. b Deterministic simulation with inputs X lo ; Y hi at 10 molecules and Z hi at 10 molecules, which does not exhibit any change over time. c Deterministic simulation resulting in Z lo (seen in light blue) based on the initial configuration X lo ; Y lo at 10 molecules and initial output Z hi at 10 molecules. d LNA simulation resulting in Z hi (seen in yellow) based on the initial configuration of input X hi ; Y hi at 10 molecules and initial output Z lo at 10 molecules. Note that some plots are overlayed but are either set to 0 or 10 molecules. (Color figure online) Chemical reaction network designs for asynchronous logic circuits 121 dual-rail implementations of logic gates in line with other designs proposed within this paper. In contrast to the gates in Soloveichik et al. (2008), our gates account for all inputs X hi ; X lo , Y hi ; Y lo and outputs Z hi ; Z lo . They are also reusable and respond to changes in input. The simplest gate, NOT, in Fig. 15a(i), inverts the inputs X hi ; X lo to outputs Y lo ; Y hi . The AND-gate, shown in Fig. 15a(iii), has initial concentrations of k 1 ; k 2 as well as input concentrations. With a presence of species Y hi we can catalyse k 2 into the state k 1 , and with the species X hi we can catalyse k 1 to Z hi ; thus both species are needed for the gate to output the signal Z hi . The state Z hi catalyses a reaction between species Z lo and k 3 , therefore showing that only one output signal can be present at any time. Conversely, with either X lo ; Y lo we can convert k 3 to Z lo , which in turn can convert Z hi back to k 2 and k 2 to k 1 . Using a similar trail of thought we can see how the other gates are devised, with OR ( Fig. 15a(ii)) having initial concentrations of k 2 ; k 3 , Fig. 12 Two latch designs and their comparison. a(i) Simple latch with two feedback loops. a(ii) Latch with reset to a neutral state. b Deterministic plot of simple latch in a(i). c Deterministic plot of the latch in a(ii) which shows faster convergence Fig. 13 Arbiter circuit design and its simulation. a The arbiter CRN with inputs X hi and X lo and outputs Y hi and Y lo . The output reflects the input with the higher concentration of molecules (or which ever species appeared first). b LNA simulation of the arbiter, demonstrated with an input of X hi at a concentration of 5 molecules and X lo with a concentration of 0 molecules. After 0.4 s we see Y hi (in red) at a concentration of 10 molecules, representing the majority input. (Color figure online) NOR ( Fig. 15a(vi)) having initial concentrations of k 2 ; k 3 and NAND (Fig. 15a(iv)) having initial concentrations of k 1 ; k 2 . We provide an example, showing a deterministic simulation of the AND gate, seen in Fig. 15b, in which inputs X hi ; Y hi at 10 molecules and initial output Z lo is converted to Z hi after 0.8 s.
XOR is slightly different. XOR, traditionally a gate that requires a composition of many other logic gates, has to be constructed with all combinations of inputs considered. The XOR gate (Fig. 15a(v)) has initial concentrations of k 1 ; k 2 ; k 3 ; k 4 . In Fig. 16 we show an example validation of the XOR gate for all four input configurations using LNA.

Muller C-pipeline
We construct a CRN to emulate the behaviour of the C-pipeline outlined in Sect. 3.2. The pipeline is a mechanism to relay handshakes between components, for example latches to store data. We construct the pipeline by placing three of our C-element CRNs, shown Fig. 11a, in series. At each intermediate stage between the C-elements we add a fork. One path of the fork is negated and fed back into the previous C-element, and the other path is fed into the new C-element. The key interaction between the components of the C-pipeline is that no C-element can Fig. 14 CRNs for control flow components. a The fork used to split a signal into two. b The join used to merge signals Fig. 15 a Dual-rail logic gate designs: we present a NOT (i), OR (ii), AND (iii), NAND (iv), XOR (v) and NOR (vi) over inputs X, Y and output Z. b Deterministic simulation for AND on inputs X hi ; Y hi at 10 molecules and initial output Z lo , which is converted to Z hi after 0.8 s (seen in grey) Chemical reaction network designs for asynchronous logic circuits 123 output a positive species or negative species without the previous displaying a positive or negative one.
The inputs to the C-pipeline CRN are Req hi ; Req lo ; Acc hi and Acc lo . The only output is C hi ; C lo , which is the output species of the third C-element. However, for the sake of Fig. 16 XOR gate validation demonstrated using LNA for all input combinations and k 1 ; k 2 ; k 3 ; k 4 having initial concentrations of 10 molecules. In a we demonstrate that XOR on an input configuration X lo ; Y lo at 10 molecules produces an output Z lo (seen in grey) at 10 molecules after 0.4 s. d A similar plot on input X hi ; Y hi , which results in output Z lo (seen in grey). In b, c we show that both X hi ; Y lo and X lo ; Y hi demonstrate the correct output of Z hi (seen in pink). (Color figure online) clarity, we also include four other species A hi ; A lo corresponding to the output of the first C-element and B hi ; B lo corresponding to the second. On an input of Req hi at 10 molecules we would expect to see that A hi , B hi and C hi are all at 10 molecules after a staggered amount of time. If we then changed the input to Req lo we would expect to see A hi , B hi , C hi diminish with A lo , B lo , C lo , reaching 10 molecules to reflect the change in input. This is seen as a 'wave' through the pipeline propagating a high signal and then a low signal.
We design an experiment, see Fig. 17, where we initialise the pipeline with the input species Req hi at 10 molecules, and all C-elements are initialized with the intermediary species k at 10 molecules such that no Celement yet outputs a species. From both the deterministic and LNA simulations of this we can observe how the species A hi ; B hi ; C hi approach 10 molecules one after another, showing that indeed the value of Req hi is being propagated along the pipeline. In order to show that our pipeline design is continuously reactive, we convert all of the species Req hi to Req lo via the introduction of a reaction Req hi ! Req lo . This effect occurs at around 1 s and we can observe that the pipeline responds by reducing the species A hi ; B hi ; C hi to a molecular count of 0. We also observe (not shown on the simplified plot) that the species A lo ; B lo ; C lo reach 10 molecules at the same time as A hi ; B hi ; C hi reach 0 molecules (2 s). The LNA plot reveals that all output species are separated by a significant time difference such that no two species can be conflated.
In addition to simulation experiments which plot expected concentrations of species over time, we also check temporal properties concerning the interactions between species and components. Using the PRISM model checker we interrogate the CTMC models of the pipeline with specific queries. We give some important examples of such queries in Table 1, which are verified by PRISM as being true with very high probability by checking 20 paths against the property. ''Th'' refers to the required population threshold which can be set by the user. PRISM also has the ability to track and plot concentrations of species over specific time intervals and number of samples. For example, for the C-pipeline we may wish to focus specifically on species concentrations of the second C-element whilst ignoring the others. We can isolate the species in question starting at a time [ 0 and simulating only the species of the second C-element. We demonstrate this property on the pipeline with initial input species Req hi at 10 molecules in Fig. 18.

Queue
We have also designed and validated a queue, the schematic for which is shown in Fig. 19, built by the addition of latches at each C-element block to the Muller pipeline. The queue is used in electronics to regulate and store the flow of information. The asynchronous queue uses the pipeline as a control mechanism to propagate signals between the Fig. 17 Validation of the Muller C-pipeline. The input request signal, encoded by the species Req hi , is propagated to the end of the pipeline (represented by the species A hi ; B hi ; C hi ); we then set the request signal to low. The pipeline then responds by the presence of A hi ; B hi and C hi diminishing to zero. In b we show that the variance is low, even for low molecular counts Chemical reaction network designs for asynchronous logic circuits 125 latches. We use the latch with reset seen in Fig. 12 for this purpose. As the species Req hi is propagated along the pipeline, it sends a signal to the queue to read and store the value in the next latch along. Each latch represents some computation that could be completed within each time interval.
For the latches of our queue we have input species Am hi and Am lo and output species Ams hi ; Bms hi ; Cms hi representing the output of each latch. Deterministic simulation the queue pipeline is shown in Fig. 19. In this experiment we propagate a 1 (represented as Am hi at 10 molecules) followed by 0 (Am lo at 10 molecules). We can observe the species Ams hi ; Bms hi , which represent the outputs of the first and second latches, noting an oscillatory pattern of cycling between 1 and 0.

Adder
We have also designed a three bit ripple-carry adder, which operates in a similar fashion to the queue but instead of latches we compose adders in series, seen in Fig. 20b. An adder, the circuit design for which is seen in Fig. 20a, is a composition of two XOR gates, two OR gates and an AND gate. The adder produces two outputs, the value of the summation and the carry. Within the ripple-carry adder, the carry output of each adder is fed into the next adder, which outputs the sum and a carry. In this way, with three adders, we can add three two-bit numbers together.
Our ripple-carry CRN has four input species per adder representing the two inputs, and two output species representing the output. In Fig. 21 we show that the adder exhibits correct behaviour by producing the desired output species for a specific input, where each sum is calculated only in the next stage in the pipeline. If we view each Celement and adder as one stage in the pipeline, labelled A, B and C, then we can view the output species of each adder as a bridge to the next adder. The six output species, representing the carry-bit output, are denoted by A abridgeOneOut ; B bbridgeOneOut ; C cbridgeOneOut and A abridgeZeroOut ; B bbridgeZeroOut ; C cbridgeZeroOut . In order to show correct operation the output of the adders (represented in this case by 10 molecules) should be interleaved with the control species of the C-pipeline (A hi ; B hi ; C hi ), allowing time for the carry species to catalyse with the input of the following adder. Whilst a simulation provides the intuition behind the ripple-carry adder, using the PRISM model checker we can query the output of all three adders after 10 s to confirm if the correct output is present. We have four input species per adder, excluding the carry, which represent two numbers. We expect one species from each adder as output, representing the addition of two inputs, plus the carry. The third adder relies on the previous adder's carry being correct. We therefore only need to look at the desired output of each adder plus the carry of the final adder. We summarise this with the following example PRISM property: With this example we can only satisfy these three outputs and the carry, by seeing each of their molecular concentrations rise above the threshold Th, based on a specific input configuration. This particular input are the species representing 0 and 1 for the first adder (for example A aOneZeroIn and A aTwoOneIn ), the species representing inputs 1,1 for the second adder and 1,1 for the third adder. 0 þ 1 in the first adder should give us an outcome of 1 carry 0 and so satisfies A abridgeOneOut [ Th. 1 þ 1 plus the 0 carry from the first adder gives an output of 0 carry 1, and so satisfies B bbridgeZeroOut [ Th. An input of 1 ? 1 plus 1 carry from the second adder means that our output should be 1 as well as the final carry should be 1, represented by C cbridgeOneOut [ Th and C cbridgeCarryOneOut [ Th. Our adder satisfies this property based upon the inputs given and therefore shows correct operation for an adder.

Experimentation
All designs 3 presented in this paper have been validated using both Microsoft's Visual GEC tool (Pedersen and Phillips 2009) and the PRISM model checker Fig. 19 Deterministic simulation of the queue pipeline. We propagate a value of 1 through the queue. The species Ams hi ; Bms hi represent the outputs of the first and second latches. Note that through oscillatory patterns generated by the pipeline we can mimick properties of a synchronous system Fig. 20 a Circuit diagram for a ripple-carry adder. b Three adders in series controlled by the C-pipeline. A carry-bit output from one adder is fed into the next as part of the input (Kwiatkowska et al. 2011), both for the deterministic and stochastic semantics of the CRNs. Visual GEC provides a programming language, LBS, for designing and simulating any given CRN. We systematically tested each component in isolation by simulating its behaviour against all input and output configurations. Next, we examined how a component might behave in a larger system, where it will be exposed to a change in input. To this end, we introduced new reactions to emulate a signal change. For instance, if we wished to change a carrier signal from high to low, we would introduce an additional reaction X hi ! k X lo , which converts all of the signal X hi into a signal X lo while the component is operating. Since deterministic semantics is not accurate for low molecular populations, we additionally explored stochastic semantics. Visual GEC exports models to the probabilistic model checker PRISM, which then enables verification of the induced continuous-time Markov chain against temporal logic properties. This allows one to check that the circuits ensure the correct temporal ordering of the events, for example, for the Muller pipeline seen in Fig. 6, that the species in the first stage of the pipeline is present before the species in the second, i.e. with probability 1, and that the signal is eventually propagated to the end of the pipeline. PRISM implements numerical solution of the CME, which is exponential in the initial number of molecules and hence not scalable, and analysis based on stochastic simulation, which is time consuming. We thus additionally used an experimental implementation of the LNA within Visual GEC, based on Cardelli et al. (2015). As well as being capable of checking temporal logic properties (Cardelli et al. 2015;Bortolussi et al. 2016), the LNA can plot the species concentration over time together with standard deviation, and is fast and reasonably accurate even for low molecule counts. Moreover, compared to the deterministic semantics, LNA provides important information about stochasticity that may affect the robustness of the circuits, and which can be explored further with CME, stochastic simulation, or verifying that the circuit converges with probability 1 to a single value.

Discussion
When modelling asynchronous circuits as a chemical system, the wires are chemical species from the output set of one gate component to the input set of another. We cannot bound the time for a gate to transform an input species to an output species. This excludes the class of self-timed circuits. Under deterministic semantics, we could guarantee an isochronous fork since two chemical species, either high or low, could theoretically reach the threshold M at precisely the same time given equal rates and initial Fig. 21 Deterministic simulation of the ripple-carry adder circuit responding to various inputs. We plot the output species from each section of the pipeline used to coordinate the output from each adder. The carry-bit output from each adder is represented by A abridgeOneOut ; B bbridgeOneOut and C cbridgeOneOut . The output of the Celement (A hi ; B hi ; C hi ) arrives strictly before the output from the adder. The logical input for a is 1 and 0 for the first adder, 1 and 0 for the second adder, and 1 and 0 for the third adder. In b we show the computation on different inputs, namely 1 and 0 for the first adder, 0 and 0 for the second adder, and 1 and 0 for the third adder. The crossover in the concentrations of output species of the C-element and the logical output of 1 (resulting from inputs 0, 1 and carry of 1) in the third adder (plots C hi and C csOneOut ) indicates faster convergence but does not affect the results in further stages of the pipeline concentrations, and therefore under deterministic semantics we have a Turing-complete method of computation. We cannot guarantee this under stochastic semantics (Cook et al. 2009). This is because there is a non-zero probability that one species could reach M before the other. We thus conclude that our circuits, at worst, can be classified as speed independent. We can calculate approximately the delay on wires based upon rates and concentrations for each semantics.
Direct chemical implementations of CRNs have been theorised and realised, but involve complicated reaction mechanisms (Shin 2011). The most common substrate for chemical kinetics is DNA strand displacement (DSD), which involves the displacement of DNA strands in solution. These strands are labelled with the chemical species and, once the reaction has taken place, an outputting strand represents an output from the CRN that the strand displacement system is trying to emulate. DNA strand displacement has been shown to be a universal substrate for chemical kinetics, specifically for bi-molecular reactions used here (Soloveichik et al. 2010). Most importantly, the AM circuit seen in Fig. 8b has been implemented as a strand displacement device (Lakin et al. 2012). However, a potential difficulty with this approach is scalability: as the number of components increases, the number of chemical species representing them also increases. Large numbers of chemical species result in large numbers of DSD complexes in solution, and consequently crosstalk needs to be accounted for. Fortunately, a recent experiment with the implementation of a square-root circuit in solution provided a new experimental ceiling on the number of species that can be used .
Another challenge is to provide formal verification of the correctness of the designs. We remark that, although we have validated the behaviour of the components for all possible input configurations and verified that correct outputs are produced and in the correct order using simulation and simulation-based probabilistic model checking with PRISM, this does not amount to full verification. Asynchronous diagrams are represented using a variety of notations, see Fig. 5, and correspond to certain classes of (safe) Petri nets (Myers 2004) known as Signal Transition Graphs (STGs), representing the rise and fall of signals. Our designs are systems containing many molecules which exhibit stochastic behaviour. Firstly, one would need to show that our CRN designs meet the specification given as a Petri net, which would involve relating the two formally via a refinement relation, where one needs to relate structures with many molecules of a given a species to structures with at most one. This presents us with two major challenges: scalability and stochasticity. Scalability can be addressed using compositional verification, which has been developed for (non-probabilistic) process algebraic specifications of asynchronous circuits (Wang and Kwiatkowska 2007) (equivalent to STGs) and it would be interesting to see if they can be applied in this setting. However, no probabilistic extension of this approach is known. Another possibility is to capture stochasticity by employing stochastic Petri nets to model the designs as done for molecular walkers in Barbot and Kwiatkowska (2015), and then perform temporal logic verification using the tool Cosmos. Cosmos relies on an implementation of statistical model checking that exploits parallelism of the Petri net specifications and achieves greater scalability than PRISM.

Conclusion
We have proposed a novel design for an asynchronous computing device based on Chemical Reaction Networks. CRNs are inherently asynchronous, and thus particularly well suited to this computational paradigm. Our designs are based on a simple, bi-molecular reaction motif inspired by Approximate Majority (Angluin et al. 2008;Cardelli and Csikász-Nagy 2012), employ catalytic reactions and assume well-mixed solution and constant, uniform rates. Moreover, they do not rely on the universal clock which is difficult to realise. Since an arbitrary CRN can be physically realised using DNA strand displacement (Soloveichik et al. 2010), as recently demonstrated experimentally in , the proposed designs are in principle implementable, and we have confirmed this in theory by modelling them in the two-domain setting (Cardelli 2010) using Visual DSD (Phillips and Cardelli 2009;Lakin et al. 2011). Our designs are the first feasible implementation of an asynchronous computing device in chemical kinetics and are relevant for a multitude of applications in nanotechnology and synthetic biology. As future work we would like to investigate alternative experimental settings.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.