Error-free stable computation with polymer-supplementedchemical reaction networksbyAllison Yeu Yang TaiB.Sc., University of British Columbia, 2016A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Bioinformatics)The University of British Columbia(Vancouver)November 2019c© Allison Yeu Yang Tai, 2019The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Error-free stable computation with polymer-supplemented chemical re-action networkssubmitted by Allison Yeu Yang Tai in partial fulfillment of the requirements forthe degree of Master of Science in Bioinformatics.Examining Committee:Anne Condon, Computer ScienceSupervisorDavid Kirkpatrick, Computer ScienceSupervisory Committee MemberRaymond Ng, Computer ScienceSupervisory Committee MemberSteven Jones, Medical Genetics and Molecular Biology BiochemistryAdditional ExamineriiAbstractWhen disallowing error, traditional chemical reaction networks are very limited incomputational power: Angluin et al. and Chen et al. showed that only semilinearpredicates and functions are stably computable by CRNs. Qian et al. and othershave shown that polymer-supplemented chemical reaction networks are capableof Turing-universal computation. However, their model requires that inputs areloaded onto the polymers at protocol initialization, in contrast with the traditionalconvention that inputs are represented by counts of molecules in solution. Here, weshow that polymer-supplemented chemical reaction networks can stably simulateTuring-universal computations even with solution-based inputs. However, suchsimulations use a unique ”leader” polymer per input type and thus involve manyslow bottleneck reactions. We further refine the polymer-supplemented chemicalreaction network model to allow for clone polymers, that is, multiple functionally-identical copies of a polymer, and provide an illustrative example of how bottleneckreactions can be reduced in this new model.iiiLay summaryIn this paper, we design algorithms for the chemical reaction networks that lever-ages the additional power of polymers to compute functions and simulate powerfulcomputational models that cannot be done traditionally. In particular, we demon-strate that we can simulate these functions and models in an error-free manner,meaning that our algorithm arrives at the correct answer every single time. This isnot possible with the traditional chemical reaction network model, which is capa-ble of computing only a very limited class of functions and predicates. Chemicalreaction networks use molecules and their reactions to perform computation, byusing a finite set of molecular species and a specified input to generate an outputthat corresponds to the answer. We also suggest methods to speed up some of thealgorithms we design, by allowing multiple copies of each polymer species to existsimultaneously.ivPrefaceThis thesis is the original, independent work of the author, A. Tai, with input onideas and writing from A. Condon. A version of the document body was publishedas conference proceedings for DNA 25 [Tai A, Condon A. Error-free stable com-putation withpolymer-supplemented chemical reaction networks. DNA Computingand Molecular Programming, 197-218, 2019].vTable of contentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Contributions and highlights . . . . . . . . . . . . . . . . . . . . 32 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Polymer-supplemented chemical reaction networks . . . . . . . . . 63 Stable, Turing-universal computation by sequential psCRNs withleader polymers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Faster computation of Square by threaded psCRNs with clone poly-mers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215 psCRN time complexity analysis and simulation . . . . . . . . . . . 30vi6 Conclusions and future work . . . . . . . . . . . . . . . . . . . . . . 36Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38viiList of tablesTable 3.1 Instruction abstractions of psCRN reactions. The decrementdec(σ ) instruction can complete only if the σ -polymer has lengthis at least 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 13viiiList of figuresFigure 5.1 (a) Simulation of Algorithm 1: Sequential-n2-psCRN, withfast-forwarding and input n on X-polymer, n = 100. (b) Sim-ulation of Algorithm 1: Sequential-n2-psCRN, with input de-tection and loading, n = 100. Each line in the plots shows thecount of the outputs as a function of the number of interac-tions, with the blue line being the count of output species Yfinally released into solution, while the red line shows the sizeof the Yint polymer. By interaction, we mean a collision of twomolecular species in the system, which may or may not resultin a reaction. . . . . . . . . . . . . . . . . . . . . . . . . . . 34Figure 5.2 Plotting simulation runtimes as a function of initial speciescount n for (a) Algorithm 1: Sequential-n2-psCRN startingwith pre-loaded input (b) Algorithm 1: Sequential-n2-psCRNwith input detection and loading. Each blue marker plots a sin-gle simulation, showing the number of interactions required forprotocol completion as a function of the initial count of inputspecies X . In both cases, data was fitted with an equation ofform cx4, where (a) c = 78.4, and (b) c = 99.2. . . . . . . . . 34ixFigure 5.3 Comparing runtimes of Algorithm 1: Sequential-n2-psCRN(red) to Algorithm 4: Threaded-n2blognc-psCRN (blue) withfast-forwarding, where (a) n = 128, and (b) n = 2i, i ∈ [3,7].We plot the average number of interactions over five differentsimulations for each input count n. Sequential data was fittedwith an equation of form cx4, c = 50.3. Threaded data wasfitted with an equation of form cx3 logx, c = 267. . . . . . . . 35xAcknowledgmentsThank you to my supervisor, Anne Condon, for her support with ideas and feed-back throughout my degree, and for helping me with the technical aspects of thewriting. I also thank my committee for their valuable feedback and insight whilepreparing this thesis.xiChapter 1IntroductionThe logical, cause-and-effect nature of chemical reactions has long been recog-nized for its potential to carry information and make decisions. Indeed, biologicalsystems exploit interacting digital molecules to perform many important processes.Examples include inheritance with DNA replication, passing information from thenucleus to the cytoplasm using messenger RNA, or activating different cellularstates through signal transduction cascades. Chemical reaction networks (CRNs)exploit these capabilities to perform molecular computations, using a finite set ofmolecular species (including designated input and output species) and reactions.Reactions among molecules in a well-mixed solution correspond to computationsteps. CRN models may use mass-action kinetics, where the dynamics of the reac-tions are governed by ordinary differential equations, or stochastic kinetics, wherethe choice of reaction and length of time between reactions depends on counts ofmolecular species. We focus on stochastic CRNs in this work.Stochastic CRNs using unbounded molecular counts are Turing-universal [19],but have a non-zero chance of failure. One challenge is that CRNs are unable todetect the absence of a molecule, and therefore when all molecules of a particu-lar species have been processed. For example, when trying to simulate a countermachine, if their count of a species corresponds to a counter value, then a test-for-zero instruction needs to detect when all molecules have been depleted, i.e., theircount is zero. Indeed, while the error of a CRN simulation of Turing-universalcomputation can be made arbitrarily small, it can never reach zero [9].1Error-free CRNs include those that exhibit stable computation: the output canchange as long as it eventually converges to the correct answer; and committingcomputation: the presence of a designated “commit” species indicates that the out-put is correct and does not change subsequently. The class of predicates stablycomputable by CRNs is limited to semilinear predicates [3], and functions com-putable by committing CRNs are just the constant functions [9].Cummings et al. [11] introduced the notion of limit-stable computation, whichrelaxes the stability requirement. In a computation of a limit-stable CRN, the out-put may change repeatedly, but the probability of changing the output from its cor-rect value goes to zero in the limit. Cummings et al. show that any halting countermachine can be simulated by a limit-stable CRN. Their construction involves re-peated simulations of a counter machine, resetting the counter values each time,along with slowing down any error-prone reactions each time they occur. Theyshow that the computational power then becomes equivalent to that of a Turingmachine with the ability to change its output a finite number of times, capable ofdeciding predicates in the class ∆02 of limit-computable functions.From these insights, we can see that CRN computations that produce the cor-rect answer with probability 1 are still severely limited. We ask, ”Is there any wayto extend CRNs to work around the lack of ability to detect absence?” Qian et al.[17] gave a promising answer to this question by introducing a CRN model thatis supplemented by polymers that behave as stacks, onto which monomers can bepushed and popped. Most importantly, this extended model allows for the stackbase unit ⊥ to be used as a reactant in implementing a ”stack empty” operation.Indeed, Qian et al. use this operation in an error-free simulation of a stack ma-chine. The resulting protocol, however, requires that the entire input is pre-loadedonto one of the stacks, a large change from traditional CRNs which assumes theinputs are well-mixed in a solution.Motivated by the work of Qian et al., we wish to go one step further: Is Turing-universal stable computation by polymer-supplemented CRNs (psCRNs) possiblewhen the input is represented by counts of monomers in a well-mixed solution? In-tuitively, if input monomers can be loaded on to a polymer, absence of that speciesfrom the system could be detected by emptiness of the polymer, and thus circum-vent a significant barrier to error-free, Turing-universal computation. The obvious2obstacle is that it seems impossible to guarantee that all inputs are loaded on thepolymer before computation can begin, if we can’t reliably check for absence ofinputs in the environment. At first glance, the logic appears circular, but we showthat indeed stable Turing-universal computation is possible, and also present ideasfor speeding up such computations.In the rest of this chapter we describe our four main contributions and reviewrelated work. Chapter 2 introduces our polymer CRN models, Chapters 3, 4, and5 describe our results, and Chapter 6 concludes with a summary and directions forfuture work.1 Contributions and highlightsStable counter machine simulation using CRNs with leader polymers. We designa polymer-supplemented CRN (psCRN) that unlike Qian et al., performs error-free Turing universal computation without requiring all inputs to be pre-loaded onpolymers. To this end, instead of a binary stack machine, we simulate a schemathat depends on unary counters. We do this by designing a protocol simulating acounter machine that performs computations assuming perfect loading of inputs,then add a mechanism that forces the simulation to continuously restart as long asinputs still exist in solution, only stopping once all inputs have been loaded ontotheir respective polymers. Therefore the protocol then completes successfully ifand only if the last simulation began with all inputs properly loaded. Our schemeis similar to the error correction scheme of Cummings et al. in how we restartthe simulation multiple times, but leverages polymers to ensure stable computationwith 0 probability of error. Our polymer simulation of counter machines, and thusTuring-universal computation, uses a unique polymer per input species, as well asa “program counter” to ensure that execution of reactions follows the proper order.These molecules each have a count of one and so in the parlance of traditionalCRNs, serve as “leaders”. As a consequence, the simulation has a high numberof so-called bottleneck reactions, which involve two leader reactants. Bottleneckreactions are undesirable because they are slow.3Clone polymers can help reduce bottleneck reactions. To reduce the number bot-tleneck reactions, we propose a CRN polymer model with no limit on the numberof polymers of a given species, other than the limit posed by the volume of thesystem. In addition to type-specific increment, decrement and jump-if-empty oper-ations (which are applied to one polymer of a given species), polymer stubs can becreated or destroyed. We call such polymers “clones.” We illustrate the potential ofpsCRNs with clone polymers to reduce bottleneck reactions, by describing psCRNto compute f (n) = n2blgnc (which is the same as n2 when n is a power of 2).Abstractions for expressing CRN multi-threading and synchronization. Our CRNfor f (n) = n2blgnc uses threading to ensure that polymer reactions can happen asyn-chronously, and uses a ”leader” polymer for periodic synchronization. To describeour psCRN, we develop threading abstractions for psCRNs with clone polymers.Time complexity and a simulator for CRNs with clone polymers. To test the cor-rectness of our psCRNs and evaluate their running times, we developed a customCRN simulator designed to support clone polymers and their associated reactions.Underlying our simulator is a stochastic model of psCRN kinetics that is a naturalextension of traditional stochastic CRNs and population protocols. We also usethis model to analyze the expected time complexities of our psCRNs examples inthis paper, showing how speedups are possible with clone polymers.2 Related workSoloveichik et al. [18] demonstrated how Turing-universal computation is possiblewith traditional stochastic CRNs, achieving arbitrarily small (but non-zero) errorprobability. For the CRN model without polymers Cummings et al. [11] showedhow to reset computations so as to correct error and achieve limit-stable computa-tion (which is weaker than stable computation).In order to understand the inherent energetic cost of computation, Bennett[4, 5] envisioned a polymer-based chemical computer, capable of simulating Tur-ing machines in a logically reversible manner. Qian et al. [17] introduced a stack-supplemented CRN model in which inputs are pre-loaded on stacks, and showed4how the model can stably simulate stack machines. Johnson et al.[15] introduce aquite general linear polymer reaction network (PRN) model, allowing for arbitrarygrowth from both ends of a polymer as well as for polymers to have monomersof multiple types. With this model, they can simulate not only the original DNAstack machine model from Qian et al., which we also focus on, but copy-tolerantTuring Machines and microtubule dynamics. However, as the focus of their re-search is on simulation and verification, and does not address the issues raised inour work, namely achieving stable computation without pre-loaded inputs and re-ducing bottleneck reactions. Cardelli et al. [6] also demonstrated Turing-universalcomputation using polymers, using process algebra systems, but these systems arenot stochastic. Jiang et al. [14] also worked on simulating computations withmass-action chemical reactions, using a chemical clock to synchronize reactionsand minimize errors.Lakin et al. [16] described polymerizing DNA strand displacement (DSD)systems, and showed how to model and verify stack machines at the DSD level.They also simulated their stochastic systems using a “just-in-time” extension ofGillespie’s algorithm. Their model has a single complex to represent a stack. Rec-ognizing limitations of this, they noted that “it would be desirable to invent analternative stack machine design in which there are many copies of each stackcomplex...”, which is what we do in this paper. They propose that updates to stackscould perhaps be synchronized using a clock signal such as that proposed by Jianget al. [14]. In contrast, our synchronization mechanism is based on detection ofempty polymers.The population protocol (PP) model introduced by Angluin et al. [1], whichis closely related to the CRN model, focuses on pairwise-interacting agents thatcan change state. In Angluin et al.’s model, agents in a PP are finite-state. Aninput to a computation is encoded in the agents’ initial states; the number of agentsequals the input size. Any traditional CRN can be transformed into a PP and viceversa. Chatzigiannakis et al. [7] expand the n agents to be Turing machines, thenexamine what set of predicates such protocols can stably compute using O(logn)memory. Although the memory capacity of our polymers can surpass O(logn),polymer storage access is constrained to be that of a counter or stack, unlike themodel of Chatzigiannakis et al.5Chapter 2Polymer-supplemented chemicalreaction networksA polymer-supplemented stochastic chemical reaction network (psCRN) modelsthe evolution of interacting molecules in a well-mixed volume, when monomerscan form polymers. We aim for simplicity in our definitions here, providing justenough capability to communicate our key ideas. Many aspects of our definitionscan be generalized, for example by allowing multiple monomer species in a poly-mer, or double-end polymer extensibility, as is done in the work of Johnson et al.[15], Lakin et al. [16], Qian et al. [17], and others.Reactions. A traditional CRN describes reactions involving molecules whose speciesare given by a finite set Σ. A reactionr+ r′ −−→ p+ p′ (2.1)describes what happens when at most one each of species r,r′ ∈ Σ are picked to bereactant molecules: they produce molecules of species p ∈ Σ and p′ ∈ Σ. In ourmodel, exactly one of r, r′ and one of p, p′ can be null, so each reaction has atleast one, but at most two, reactants and products. When both r, r′ are not null, asingle r molecule and single r′ are assumed to have interacted or collided, and thereaction is bimolecular. Otherwise we say the non-null molecule decomposes into6p, p′, one of which may be null, and the reaction is unimolecular. We assume thatfor all valid reactions, the multi-sets {r,r′} and {p, p′} are not equal, that for any rand r′, there is at most one reaction with reactants of species r and r′. For now wedo not ascribe a rate constant to a reaction; we will do that in Chapter 5.Polymer-supplemented chemical reaction networks (psCRNs) also have reac-tions pertaining to polymers. A designated subset Σ(m) of Σ is a set of monomers. Apolymer of type σ ∈ Σ(m), which we also call a σ -polymer, is a string ⊥σσ i, i≥ 0;its length is i and we say that the polymer is empty if its length is 0. We call ⊥σa stub and let ⊥= {⊥σ | σ ∈ Σ(m)} ⊂ Σ. [16, 17]. Polymer reactions also involvemolecules in a setA = {Aσpush | σ ∈ Σ(m)}∪{Aσpop | σ ∈ Σ(m)} of guard molecules,whereA ⊆ Σ−Σ(m). Generally, for each σ ∈ Σ(m) there are two polymer reactions,corresponding to σ -push and σ -pop, respectively:[⊥σ . . . ]+σ −−→ [⊥σ . . .σ ]+Aσpush[⊥σ . . .σ ]+Aσpop −−→ [⊥σ . . . ]+σA σ -push consumes a single monomer σ and a σ -polymer of length i to generatea guard Aσpush and a σ -polymer of length i+ 1. A σ -pop instead does almost thereverse, consuming a single guard Aσpop and a σ -polymer of length i to generate amonomer σ along with a σ -polymer of length i−1. We cover the details of thesereactions in the section discussing configurations.Configurations. A configuration specifies how many molecules of each speciesare in the system at a specified time, keeping track also of the lengths of all σ -polymers. Formally, a configuration is a mapping c : Σ∪ {⊥σσ i | i ≥ 1} → N,where N is the set of nonnegative integers. We let c([⊥σ . . .]) denote total numberof σ -polymers in the system at a given time (including stubs) and let c([⊥σ . . .σ ])denote the total number of σ -polymers in the system that have length at least 1.With respect to configuration c, we say that a molecule of species σ ∈ Σ is a leaderif c(σ) = 1.A reaction of type (2.1) is applicable to configuration c if, when r 6= r′, c(r)≥ 1and c(r′) ≥ 1, and when r = r′, c(r) ≥ 2. If the reaction is applied to c, a newconfiguration c′ is reached, in which the counts of r and r′ decrease by 1 (when7r = r′ the count of r decreases by 2), the counts of p and p′ increase by 1 (whenp = p′ the count of p increases by 2), and all other counts remain unchanged.We also for now add a special restriction to all σ -polymer molecules: the onlyreactions they can be involved in are σ -pushes, σ -pops, and bimolecular reactionswhere they are not consumed or generated, such as Li+⊥σ −−→ Lk +⊥σ . At thesame time, we enforce that c0([⊥σ . . .]) = 1 for all [⊥σ . . .], restricting the totalnumber of any σ -polymer to be 1, meaning every polymer in the system is a leaderpolymer. We defer the case of clone polymers, where ∃[⊥σ . . .], c([⊥σ . . .])> 1, toChapter 4.Under the above, a σ -push is applicable if and only if c(σ)> 0. The result ofapplying a σ -push reaction is that c(σ) decreases by 1, c(Aσpush) increases by 1 andalso for some i ≥ 0 such that c(⊥σσ i) = 1, c(⊥σσ i) is set to 0 and c(⊥σσ i+1) isset to 1.On the other hand, a σ -pop is applicable if and only if c([⊥σ ])= 0 and c(Aσpop)>0. This differs significantly from the traditional pops seen when computing withstacks, as unlike there, where only the stack itself is required for the pop, the pres-ence of at least one Aσpop is absolutely necessary. Similarly, the result of applyinga σ -pop reaction is that c(σ) increases by 1, c(Aσpop) decreases by 1, and for somei≥ 1 such that c(⊥σσ i) = 1, c(⊥σσ i) is set to 0 while c(⊥σσ i−1) is set to 1.Intuitively, after these reactions, the length of the σ -polymer in the system ei-ther grows or shrinks by 1, and correspondingly, the count of Aσpush either increasesby 1, or Aσpop decreases by 1.If c′ results from the application of some reaction to c, we write c→ c′ and saythat c′ is directly reachable from c. We say that c′ is reachable from c if for somek ≥ 0 and configurations c1,c2, . . . ,ck,c→ c1→ c2 . . .→ ck→ c′.Computations and stable computations. We’re interested in CRNs that compute,starting from some initial configuration c0 that contains an input. For simplicity,we focus on CRNs that compute functions f : Nk→ N. For example, the functionmay be Square, namely f (n) = n2.8In a function-computing psCRN, the input n = (n1, . . . ,nk) ∈Nk is representedby counts of species in a designated setI = {X1,X2, . . . ,Xk}⊆ Σ(m) and the outputis represented by the count of a different designated species Y ∈ Σ(m). In the initialconfiguration c0 = c0(n), the initial counts of the input species Xi is ni, 1 ≤ i ≤ k,and the counts of all species other than the input species, including polymers andguard molecules, is 0, with the exception that there may be some leader moleculesor polymers present. A leader molecule can function as a program counter, forexample.A computation of a psCRN is a sequence of configurations starting with aninitial configuration c0, such that each configuration (other than the first) is directlyreachable from its predecessor.Let C be a psCRN, and let c be a configuration of C . We say that c is stableif for all configurations c′ reachable from c, c(Y ) = c′(Y ), where Y is the outputspecies. The psCRN stably computes a given function f : Nk→ N if on any inputn ∈ Nk, for any configuration c reachable from c0(n), a stable configuration c′ isreachable from c and moreover, c′(Y ) = f (n). Finally if psCRN C stably computesa given predicate, we say that C is committing if C has a special “commit” speciesLH such that for all n ∈ Nk, for any configuration c reachable from c0(n) thatcontains species LH , if c′ is reachable from c then c′ also contains LH and c′(Y ) =f (n).To summarize, a polymer-supplemented chemical reaction network (psCRN)function computer is a 9-tuple C = (Σ,Σ(m),I ,⊥,A ,R,L ), where− Σ(m) ⊆ Σ is the set of monomer types,− I ⊆ Σ is the set of input types,− ⊥= {⊥σ | σ ∈ Σ(m)} is the set of stub types,− A = {Aσpush | σ ∈ Σ(m)∪{Aσpop | σ ∈ Σ(m)} is the set of guard types,− the sets Σ(m), ⊥, and A are disjoint,−L ⊆ Σ is the set of initial leader species,− Σ contains the sets Σ(m),I ,⊥,A andL ,− Σ also contains Y , and LH , the output, and commit species types, and−R is a set of reactions, including σ -push and σ -pop for all σ ∈ Σ(m).9Bottleneck reactions. In our CRN algorithms of Chapter 3, many reactions in-volve a sole leader ”state” molecule, or program counter, that reacts with a soleleader polymer. Such reactions, in which the count of both reactants are 1, are of-ten described as bottleneck reactions [8]. As explained in Chapter 5, in a stochasticCRN that executes in a well-mixed system with volume V , the expected time forsuch a reaction isΘ(V ) [18]. Our motivation for the clone polymer model in Chap-ter 4 is to explore how to compute with polymers in a way that significantly reducesbottleneck reactions.10Chapter 3Stable, Turing-universalcomputation by sequentialpsCRNs with leader polymersHere we describe how psCRNs with leader polymers can stably simulate countermachines, thereby achieving Turing-universal computation. Before doing so, wefirst introduce psCRN “pseudocode” which is convenient for describing psCRNalgorithms. Then, as an illustration, we describe a psCRN to compute the Squarefunction f (n) = n2, under a special condition: we “fast-forward” from initial con-figuration c0 to some configuration cs where all input-loading onto polymers hasfinished, meaning ∀σ ∈ Σ(m), c(σ) = 0, and only run the protocol starting fromthat configuration.We then describe afterwards how to guarantee that we can reach this specialconfiguration cs, by coupling it with a simple counter machine simulation thatbuilds strongly on Cummings et al. [11]’s register machine illustration. Notethat although our method, like theirs, relies on re-simulating the protocol multipletimes, we only re-simulate one final time after we reach cs, upon which the prob-ability of error drops to 0, thus achieving stable computation, while Cummings etal. must keep repsimulating indefinitely, with the error of the simulations fallingover time, but never reaching 0.11Sequential psCRN pseudocode. Following earlier work [2, 11, 17, 18], we de-scribe a psCRN program as a set of instructions. Because one instruction mustfinish before moving on to the next, we call these sequential psCRNs. Correspond-ing to each instruction number i is a molecular “program counter” species Li ∈ Σ.One copy of L1 is initially present, and no other Li′ for i′ 6= i is initially present.The instructions inc(σ ) and dec(σ ) of Table 3.1 increase and decrease thelength of a σ -polymer by 1, respectively, making it possible to use the polymersas counters. The inc instruction produces a σ which is then pushed on to theσ -polymer. In order to ensure that the push and pop reactions happen only withinthe inc() and dec() reactions, the inc(σ ) operation generates the guard Aσpush ,which we convert into a separate molecule Iσ ∈ Σ−Σ(m) that cannot be used to poppolymers before the instruction execution completes, while the dec(σ ) instructionchanges Iσ into Aσpop in order to reduce the length of a σ -polymer by 1. Whena psCRN executes instructions of Table 3.1, starting from an initial configurationin which there is no polymer of length greater than 0, then we have the followinginvariant:Invariant: Upon completion of any instruction, the sum of the countof Iσ , Aσpush , and Aσpop equals the length of the σ -polymer.Even more specifically, within this chapter, by using only dec(σ ) and inc(σ )to pop and push σ onto polymers, upon the completion of any instruction, c(Aσpush)and c(Aσpop) are always 0, and therefore the sum of the count of Iσ equals thelength of the σ -polymer. Note that a σ -polymer is empty (has a length of 0) if andonly if a stub ⊥σ is present in the system. Therefore assuming that our invariantholds, the leader σ -polymer is not empty if and only if at least one Iσ moleculeis in the system. Therefore the jump-if-empty instruction is very useful inensuring that the program counter advances properly: a danger is that when theσ -polymer is empty, any dec(σ ) instructions cannot proceed and may cause analgorithm to stall. The jump-if-empty(σ ,k) instruction provides a way to firstcheck whether the σ -polymer is empty, and if not, dec(σ ) can safely be used. Thecreate and destroy instructions provide a way to create and destroy copiesof a species. While often more than one reaction is needed to implement oneinstruction, all will have completed when the instruction has completed and the12program counter is set to the number of the next instruction to be executed in thepseudocode.i: inc(σ) Li −−→ L∗i +σσ + [⊥σ . . .]−−→ Aσpush + [⊥σ . . .σ ]L∗i +Aσpush −−→ Li+1+ Iσi: dec(σ) Li+ Iσ −−→ L∗i +AσpopAσpop + [⊥σ . . .σ ]−−→ σ + [⊥σ . . .]L∗i + σ −−→ Li+1i: jump-if Li +⊥σ −−→ Lk +⊥σ-empty(σ ,k) Li + Iσ −−→ Li+1 + Iσi: goto(k) Li −−→ Lki: create(σ) Li −−→ Li+1+σi: destroy(σ) Li + σ −−→ Li+1i: halt Li −−→ LHTable 3.1: Instruction abstractions of psCRN reactions. The decrementdec(σ ) instruction can complete only if the σ -polymer has length is atleast 1.Pseudocode instructions may also be function calls, where a function is it-self a sequence of instructions expressed as pseudocode. Suppose again that thereis a leader σ -polymer and also a leader σ ′-polymer in the system. Then theadd-to(σ ,σ ′) function (using a temporary τ-polymer) extends the length of theσ ′-polymer by the length of the σ -polymer. Another useful function is flush(σ )which decrements the (leader) σ -polymer until its length is 0. A third function,release-output(σ ), is useful to “release” molecules on a (leader) σ -polymeras Y molecules into the solution. This function uses an additional special leaderY ′-polymer which is empty in the initial configuration, and whose length at the endof the function equals the number of released Y molecules. The Y ′ molecule willbe useful later, when we address how a psCRN can be restarted (and should not beused elsewhere in the code).13i: add-to(σ ,σ ′) i : goto(i.1)i.1: jump-if-empty(σ , i.6)i.2: dec(σ)i.3: inc(σ ′)i.4: inc(τ)i.5: goto(i.1)i.6: jump-if-empty(τ, i.10)i.7: dec(τ)i.8: inc(σ)i.9: goto(i.6)i.10: goto(i+1)i: flush(σ) i: goto(i.1)i.1: jump-if-empty(σ , i+1)i.2: dec(σ)i.3: goto(i.1)i: release-output(σ) i: goto(i.1)i.1 jump-if-empty(σ , i+1)i.2 dec(σ )i.3 inc(Y ′)i.4 create(Y )i.5 goto(i.1)Numbering of function instructions. For clarity, we use i.1, i.2, and so on to la-bel the lines of a function called from line i of the main program. Upon such afunction call, the CRN’s program counter first changes from Li to Li.1. The pro-gram counter is restored to Li+1 upon completion of the function’s instructions,e.g., via a goto(i+1) instruction or a jump-if-empty(σ , i+1) instruction. Ifone function fB is called from line a.b of another function fA, the program counterlabels would be a.b.1, a.b.2 and so on, and so the label “i” in the function descrip-tion should be interpreted as “a.b”. In this case, when the function fB completes,control is passed back to line a.(b+ 1) of function fA; that is, the “goto(i+ 1)”14statement should be interpreted as “goto(a.(b+ 1))”. Also for clarity, we usespecial labeling of instructions in a few special places, such as the restart functionbelow, in which instructions are labeled S1, S2 and so on.psCRNs with fast-forwarding. As noted in the introduction, a challenge in achiev-ing stable computation with psCRNs is detecting the absence of inputs. To buildup to our methods for addressing this challenge, we first ignore this issue, fast-forwarding to a new configuration cs directly reachable from c0. By this we meanthat if the input contains ni molecules of a given species Xi, then in configurationcs there is a unique Xi-polymer of length ni. Furthermore, there are ni copies of themolecule IXi in the system. Intuitively, the configuration cs is one that would bereached if ni inc(Xi) operations were performed from an initial configuration withno inputs and an empty Xi-polymer, for every input species Xi.A committing, sequential psCRN with fast-forwarding for Square. Our psCRNfor the Square function f (n) = n2 has one input species X and one output speciesY , and Σ(m) = {X ,X ′,X ′′,Yint ,Y ′,τ}. In the fast-forwarded configuration, the inputis represented as the length n of a leader X-polymer, and the count of IX is n. Theonly other molecules in the initial configuration are the leader program counterL1, and stubs ⊥X ′ , ⊥X ′′ , ⊥Yint , and ⊥Y ′ , and ⊥τ . Note that instead of using the X-polymer directly, we first copy it to the X ′-polymer in a line not numbered as partof the main program, line P. The distinction will become useful when we removethe assumption of fast-forwarding later. The psCRN has a loop (implemented usingjump-if-empty and goto) that executes n times, adding n to an intermediateYint-polymer each time. When the loop completes, the output is released from theYint-polymer in the form of Y , so that the number of Y ’s in solution is n2, and thepsCRN halts. The halting state is in effect a committing state, since no transitionis possible from LH .Committing Turing-universal computation by psCRNs with fast-forwarding. Turing-universal computation is possible with a two counters machine (2CM). To simulatea halting 2CM that computes function f :Nk→N using unary counters, we create a15Algorithm 1 Sequential-n2-psCRN, with input n on X-polymer.P: add-to(X ,X ′)1: add-to(X ′,X ′′)2: jump-if-empty(X ′,6)3: dec(X ′)4: add-to(X ′′,Yint)5: goto(2)6: release-output(Yint)7: haltpsCRN has 2 unary counters R′1,R′2, where R′1 contains the input n in unary, whileR′2 is initially 0. Throughout the simulation of the counter machine, the psCRNhas exactly one R′1-polymer and one R′2 polymer, to represent the aforementionedunary counters. In addition, there is one additional polymer, a Y ′-polymer, whichis initially empty and is used by the release-output function. A two countersmachine program is a sequence of instructions, where instructions can incrementa counter; decrement a non-empty counter; test if a counter is empty (0) and jumpto a new instruction if so, or halt. Table 3.1 already shows how all four of theseinstructions can be implemented using a psCRN. We assume in what follows thatthese are the only instructions used by the psCRN simulator. At the end, the outputis stored on counter R′1, then released into solution, using release-output(R′1),once the machine being simulated reaches its halt state.We assume that release-output(R′1) is the only function call of the 2CMsimulator.Stable, Turing-universal computation by psCRNs. We now handle the case wherewe start from initial configuration c0, where the input is represented as countsof molecules instead. That is, in the initial configuration of the psCRN all poly-mers are empty, including the R′1 and R′2-polymers; instead, we have n copies ofmolecule R1 in solution, and a new R1-polymer that is not directly part of the 2CMsimulation for them to be loaded on. Our scheme uses the R1-push reaction to loadinputs. Note that here, R1 counter, although not used by the 2CM simulator forcomputation, serves a very important role as a back-up counter: it is used to restorethe value of R′1 in case the simulation began too early, before input loading finishes.16We will demonstrate how, by adding CRNs to detect input-loading and to restartthe simulator in this case. Once all inputs are loaded, and the protocol starts onefinal time, the system is never subsequently restarted. Overall our simulation hasfour components:• Input loading: This is done as R1-push, which can happen at any time untilall inputs are loaded. Recall that the R1-push reaction is[⊥R1 . . . ]+R1 −−→ [⊥R1 . . .R1]+AR1 pushEach such reaction generates a guard molecule AR1 which, as explained be-low, triggers input detection.• Two counters machine (2CM) simulation: Algorithm 2 shows the simu-lator. This psCRN program has a “prelude” line P, that copies the inputs onthe R1-polymer to the R′1-polymer. Then starting from line numbered 1, the sim-ulation does the computation as described before using the R′1 and R′2-polymers.Upon completion of the computation, the output is released from counter R′2, andthe simulator halts (produces theLH species).Algorithm 2 Sequential-2CM-psCRN.P: add-to(R1,R′1)1: // Rest of psCRN simulation pseudocode here, using2: // R′1 and R′2: . . .: // ending with release-output(R′2) function and halt instruction.• Input detection: This is triggered by the presence of an active guard moleculeAR′1 push . For each value i of the main program counter, after the prelude add-to,as well as for LH , we have the following reactions, where S1 is the first numberof the restart pseudocode (see below). The reactions convert the guard AR1 pushmolecule into a molecule, IR1 , since the input molecule is now loaded, and alsochanges the program counter to LS1, which triggers restart.Li+AR1 push −−→ LS1 + IR1 ,LH +AR1 push −−→ LS1 + IR117LH is no longer a committing species, since it may change to LS1.Line P is a function call to add-to(R1,R′1), which executes intructions numberedP.1 through P.11 of add-to. Input detection is only done at line P.2, the firstjump-if-empty instruction:LP.2+AR1 push −−→ LP.2+ IR1It does not trigger a restart, but simply converts the guard molecule AR1 push into IR1 .• Restart: Restart happens a number of times that is exactly the total input lengthn1 + n2 + . . .nk, since each input molecule is detected by the system exactly once,generating one guard molecule. The counters R′1 and R′2 are flushed, and any out-puts that have been released in solution are destroyed, assuming that the number ofoutputs released into the solution was tracked by some Y ′ counter, as before. Thenthe program counter is set to line P of the simulator (leader molecule LP). Algorithm3 shows the restart pseudocode.Algorithm 3 RestartS1: flush(R′1)S2: flush(R′2)S3: destroy-output()S4: goto(P)i: destroy-output() i: goto(i.1)i.1: jump-if-empty(Y ′,i.5)i.2: dec(Y ′)i.3: destroy(Y )i.4: goto(i.1)i.5: goto(i+1)Correctness of Algorithm 2: Sequential-RM-psCRN. We claim that our complete2CM simulator stably computes the same function as the counter machine. (Wenote that no “fairness” assumption regarding the order in which reactions happenis necessary to show stability, since stability is a “reachability” requirement.)18The delicate part of the simulation lies in the instruction on line P, which copiesinput counter R1 to R′1.Unlike the previous chapter, the using the R1-push reaction for input loadingmeans that the count of IR1 is no longer guaranteed to equal the length of the(leader) R1-polymer. Instead, we have that IR1 is less than or equal to the lengthof the R1-polymer. This is because when input monomer R1 is pushed onto itspolymer, the AR1 push guard molecule is produced, not the IR1 molecule, and the AR1molecule is not converted to IR1 until input detection happens.This delay in converting AR1 push and generating IR1 can cause the jump-if-emptyinstruction numbered P4.2 in the add-to function to stall, when there is no IR1 andalso no ⊥R1 (since the R1-polymer is not empty due to input loading). In this case,the input detection reaction (introduced above)LP.2+AR1 push −−→ LP.2+ IR1will convert AR1 push to IR1 . This averts stalling, since jump-if-empty can pro-ceed using IR1 . The subsequent lines of the add-to(R1,R′1) code can then proceed.Once line P has completed, the correctness of the psCRN simulation of the2CM, using the copies R′1 and R′2, is not affected by input loading. Input loadingand input detection can also proceed. These are the only viable reactions from the“halting” state LH , and so eventually (since the RM machine being simulated isa halting machine), on any sufficiently long computation path, all inputs must beloaded and detected. Input detection after the prelude phase produces a “missing”Iσ molecule and triggers a restart. The restart flushes all counters used by thesimulator, and also, using the Y -polymer, destroys any outputs in solution. (Sincerestart is triggered only when the program counter is at a line of the main program,restart does not interrupt execution of the release-output function.) A newsimulation is then started at line P, ensuring that any inputs that have been loadedsince the last detect are copied to the simulator’s input counters R′1.Once all inputs have been detected, the invariant is restored and the simula-tor proceeds correctly, producing the correct output. This correct output is neversubsequently changed, and so the computation is stable.19Bottleneck reactions. In our sequential psCRNs, both inc(σ ) and dec(σ ) con-tain bottleneck reactions, and so add-to(σ ,σ ′) has Θ(|σ |) bottleneck reactions.Thus the psCRN for Square has Θ(n2) bottleneck reactions. In the next chapter weshow how to compute a close variant of the Square function with fewer bottleneckreactions, using clone polymers rather than leader polymers.20Chapter 4Faster computation of Square bythreaded psCRNs with clonepolymersConfigurations when considering clone polymers. Previously we enforced that∀[⊥σ ...], c([⊥σ ...]) = 1 by enforcing that the only reactions involving polymerswere σ -pushes and σ -pops. Now, we will remove this requirement, and let c([⊥σ ...])be any non-negative integer. We do this by introducing two new instructions:create-polymer(σ ) and destroy-polymer(σ ).i: create-polymer(σ ) Li −−→ Li+1 +⊥σi: destroy-polymer(σ ) Li +⊥σ −−→ Li+1create-polymer(σ ) generates an empty polymer while destroy-polymer(σ )consumes a empty polymer. For simplicity, we also enforce that the initial con-figuration c0 of any psCRN contain no polymers, that is, ∀[⊥σ ...], c0([⊥σ ...]) =0, and so any polymers that are needed must be explicitly created during com-putation: for the polymer-dependent functions we introduced in Chapter 3, weadd create-polymer and destroy-polymer as needed. As an example,add-to(σ ,σ ′) now requirescreate-polymer(τ) before line i.1 and create-polymer(τ) before line21i.10.With this, the mechanics of the σ -pushes and σ -pops also change. Now, theresult of a σ -push is the same except that for exactly one i≥ 1 such that c[⊥σσ i]>0, c[⊥σσ i+1] increases by 1, while c[⊥σσ i] decreases by 1. In the same way,the result of a σ -pop now has that for exactly one i ≥ 1 such that c[⊥σσ i] > 0,c[⊥σσ i−1] increases by 1, while c[⊥σσ i] decreases by 1. Intuitively, the length ofsome arbitrary σ -polymer in the system either grows or shrinks by one, and manyσ -polymers may be present in the system simultaneously.As we will see in Chapter 5, the time complexity of σ -pops and σ -pushesare now different as they are no longer always simple bottleneck reactions; thisis due to the fact that all clone polymers are indistinguishable to the system: inboth cases, the affected polymer is chosen nondeterministically from among theclones. Exactly how the polymer is chosen is not important in the context of stablecomputation, as long as every σ -polymer with positive count in c has positiveprobability of being chosen. For example, the polymer could be chosen uniformlyat random, consistent with the model of Lakin and Phillips [16], and in this work,this is what we assume.These modifications mean that inc(σ) and dec(σ) instructions can nowoperate on any of the many functionally-identical clone σ -polymers that may be inthis system, rather than just a single leader polymer, and thus reduce the number ofmandatory bottleneck reactions. Now we demonstrate the increased power arisingfrom these changes, using the function f (n) = n2blgnc as an example, and focusonly on a system that has been fast-forwarded. Input detection and loading can belayered on, in a manner similar to Chapter 3.Under these conditions, we start with a single Y -polymer of length n, and wewish to create a total of 2blgnc Y -polymers, whose lengths sum to n2blgnc. Algo-rithm 4 proceeds in blgnc rounds (lines 6-9), doubling the number of Y ’s on eachround. To keep track of the Y molecules, we introduce a distributed σ -counter datastructure, and use it with σ = Y . The data structure consists of σ -polymers thatwe call σ -thread-polymers, plus a thread-polymer counter Tσ , which is a leaderpolymer whose length, |Tσ |, is the number of σ -thread-polymers. The value of thisdistributed counter is the total length of all σ -thread-polymers. We explain belowhow operations on this distributed counter work.22Algorithm 4 Threaded-n2blognc-psCRN.P1: create-polymer(H)P2: add-to(X ,H)P3: create-distributed-counter(Y )P4: add-thread-polymer(Y ,1)P5: add-to(X ,Y)1: halve(H)2: jump-if-empty(H,5)3: double(Y )4: goto(1)5: haltAlgorithm 4 counts the number of rounds using a leader H-polymer, whoselength is halved on each round. The halve function is fairly straightforward toimplement, using instructions and functions already introduced in Chapter 3.i: halve(H) i: goto(i.1)i.1: create-polymer(H ′)i.2: add-to(H,H ′)i.3: jump-if-empty(H ′, i.9)i.4: dec(H ′)i.5: dec(H)i.6: jump-if-empty(H ′, i.9)i.7: dec(H ′)i.8: goto(i.3)i.9: destroy-polymer(H ′)i.10: goto(i+1)The double(σ ) function of Algorithm 4 is where we leverage our distributedY -counter (with σ = Y ). Recall that a distributed σ -counter data structure consistsof a set of clone σ -polymers, which we call σ -thread-polymers, plus a thread-polymer counter Tσ , which is a leader polymer whose length is the number of σ -thread-polymers. The double function first creates two other distributed countersτ and τ ′ (lines i.1 and i.2), and gives each the same number of thread-polymersas σ , namely |Tσ | thread-polymers (lines i.3 and i.4), all of which are empty. The23heart of double (line i.5) transfers the contents of the distributed σ -counter toτ and τ ′, emptying and destroying all σ -thread-polymers in the process. It thencreates double the original number of (empty) σ -thread-polymers (lines i.6 andi.7; note that the number of threads of τ is the original value of |Tσ |). It finallytransfers the τ and τ ′ polymers back to σ (lines i.8 and i.9), thereby doubling σ .i: double(σ ) i.1 create-distributed-counter(τ)i.2 create-distributed-counter(τ ′)i.3 add-thread-polymers(τ , Tσ )i.4 add-thread-polymers(τ ′, Tσ )i.5 transfer(σ , τ , τ ′)i.6 add-thread-polymers(σ , Tτ )i.7 add-thread-polymers(σ , Tτ ′)i.8 transfer(τ , σ )i.9 transfer(τ ′, σ )i.10 destroy-distributed-counter(τ)i.11 destroy-distributed-counter(τ ′)i.12 goto(i+1)Next are details of instructions used to create an empty distributed counter, andto add empty threads to the counter. Again, these are all straightforward sequentialimplementations (no threads), using leader polymers to keep track of counts.24i: create-distributed-counter(σ ) i: goto(i.1)// Creates an empty counter i.1 create-polymer(Tσ )// with zero polymers i.2 goto(i+1)i add-thread-polymers(σ ,T ) i goto(i.1)// Adds |T | empty i.1: create-polymer(Temp)// thread-polymers to the i.2: add-to(T ,Temp)// distributed σ -counter, i.3: jump-if-empty(Temp, i.8)// where T is a counter i.4: dec(Temp)i.5: create-polymer(σ )i.6: inc(Tσ )i.7: goto(i.3)i.8: destroy-polymer(Temp)i.9: goto(i+1)i add-thread-polymer(σ ,1) i goto(i.1)// Adds one empty i.1: create-polymer(σ )// thread-polymer to the i.2: inc(Tσ )// distributed σ -counter i.3: goto(i+1)The transfer function transfers the value of a distributed σ -counter to twoother distributed counters called τ and τ ′. In line i.2 of transfer, functioncreate-threads creates Tσ identical ”thread” program counters, Lt . Onceagain this is straightforward, using a leader polymer to keep track of counts. All ofthe thread program counters execute the thread-transfer function in line i.4of transfer, thereby reducing bottleneck reactions (details below). The “main”program counter, now at line i.4 of the transfer function, can detect when allthreads have completed, because each decrements Thread-Count exactly once, andso Thread-Count has length zero exactly when all threads have completed. Atthat point, the main program counter progresses to line i.5, destroying the threadprogram counters using the destroy-threads function (not shown, but usesthe destroy function to destroy each single thread).The function transfer(σ ,τ), not shown but used in double, is the same astransfer(σ ,τ ,τ ′), except the call to thread-transfer does not include τ ′and the ”inc(τ ′)” line is removed in the implementation of thread-transfer.25i: transfer(σ ,τ,τ ′) i: goto(i.1)// transfer σ to i.1: create-polymer(Thread-Count)// both τ and τ ′ i.2: create-threads(Tσ , Lt)i.3: add-to(Tσ , Thread-Count)i.4: loop-until-empty(Thread-Count, i.5)thread-transfer(σ ,τ,τ ′,Thread-Count)i.5: destroy-threads(Lt)i.6: destroy-polymer(Thread-Count)i.7: goto(i+1)i: create-threads(Tσ ,Lt) i: goto(i.1)// create Tσ thread i.1: create-polymer(Temp)// program counters, Lt i.2: add-to(Tσ , Temp)i.3: jump-if-empty(Temp, i.7)i.4 dec(Temp)i.5 create(Lt)i.6 goto(i.3)i.7: destroy-polymer(Temp)i.8: goto(i+1)i : loop-until-empty(σ ,k) Li+⊥σ −−→ Lk +⊥σFinally, we describe how threads work in thread-transfer. The threadon()instruction executes |Tσ | times, one per copy of Lt , thereby creating |Tσ | Lt1 pro-gram counters that execute computation “threads”. Using the instructiondec-until-destroy-polymer, each thread repeatedly (zero or more times)decrements one of the σ -thread-polymers and then increments both τ and τ ′. Thiscontinues until the thread finds an empty σ -thread-polymer, i.e., the stub ⊥σ , inwhich case it destroys the stub and moves to line t.5. The dec(σ ) and inc(σ )functions of Chapter 3 work exactly as specified, even when applied to distributedcounters. A key point is that the threads work “clonely” with the thread-polymers;it is not the case that each thread “owns” a single thread-polymer. Accordingly,one thread may do more work than another, but in the end all thread-polymers areempty.26A thread exits the dec-until-destroy-polymer loop by destroying ex-actly one σ -polymer. Since at the start of thread-transfer the number ofσ -thread-polymers equals the number of thread program counters, all thread pro-gram counters eventually reach line t.5, and there are no σ -thread-polymers onceall threads have reached line t.5 of the code. At line t.5, each thread decrementsThread-Count, and then stalls at line t.6. Moreover, once all threads have reachedline t.6, polymer ThreadCount is empty. At this point, the program counter fortransfer changes from line i.4 to line i.5, and all thread program counters aredestroyed.i: thread-transfer(σ ,τ ,τ ′,Thread-Count) i: threadon()t.1: dec-until-destroy-polymer(σ , t.5)t.2: inc(τ)t.3: inc(τ ′)t.4: goto(t.1)t.5: dec(Thread-Count)t.6:i: threadon(): Li+Lt −−→ Li+Lt.1i: dec-until-destroy-polymer(σ ,k) Li+ Iσ −−→ L∗i +AσAσ + [⊥σ . . .σ ]−−⇀↽− σ + [⊥σ . . .]L∗i +σ −−→ Li+1Li+⊥σ −−→ L∗∗iL∗∗i + ITσ −−→ L∗∗∗i +ATσATσ + [⊥Tσ . . .Tσ ]−−⇀↽− Tσ + [⊥Tσ . . .]L∗∗∗i +Tσ −−→ LkCorrectness. We claim that on any input n ≥ 0, as long as it’s been loaded be-forehand on a leader X-polymer using inc(X), Algorithm 4: Threaded-n2blgnc-27psCRN eventually halts with the value of the distributed-Y -counter being f (n) =n2blognc.The algorithm creates and initializes H to be a polymer of length n (lines 1-2),and the Y -distributed-counter to have a single polymer-thread of length n (lines3-5). When n = 0, H is empty, so from line 6 the algorithm jumps to line 10 andhalts, with the value of Y being f (0) = 0 as claimed.Suppose that n > 0. Reasoning about the halve function is straightforward,since it is fully sequential. We claim that in each round of the algorithm (lines 6-9),lines 7 and 8 complete successfully, with |H| halving (that is, |H| → b|H|/2c) inline 7, and with both the value of Y and |TY |, the number of Y -thread-polymers,doubling in line 8. As a result, |H|= 0 after blgnc rounds and the algorithm haltswith value(Y ) = f (n).Correctness of the double function is also straightforward to show, if weshow that the transfer(σ ,τ,τ ′) (and the transfer(σ ,τ) variant) works cor-rectly.Line i.4 is the core of transfer. We show that line i.4 does complete, that is,Thread-Count does become empty, that execution of line i.4 increases the valuesof distributed counters τ and τ ′ by the value of σ (while leaving the number of τ-and τ ′-thread-polymers unchanged), and also changes value(σ ) and the number ofσ -thread-polymers to 0.The loop-if-empty instruction ensures that the main program counter muststay at line i.4 of function transfer until Thread-Count is empty. Meanwhile,this main program counter can also activate threads using the threadon() func-tion, that is, change the thread program counters from Lt to Lt.1. From line i.2 oftransfer, the number of such thread program counters is |Tσ |.Each of these program counters independently executes thread-transfer.At line t.1, either (i) a dec(σ ) is performed (first three reactions ofdec-until-destroy-polymer), or (ii) a σ -polymer-thread is destroyed andthe polymer-thread-count Tσ is decremented (last four reactions). In case (i),both τ and τ ′ are incremented (lines t.2 and t.3), and the thread goes back to thedec-until-destroy-polymer instruction. In case (ii), the thread moves toline t.4, decrements Thread-Count exactly once, and moves to line t.5.Because the number of threads equals the value of Thread-Count at the start28of the loop-until-empty (line i.4), and because the main program counter can’tproceed beyind line i.4 of the transfer function until Thread-Count is zero, allthreads must eventually be turned on each of these threads must reach line t.4 andmust decrement Thread-Count. Only then can the main program counter proceedto line i.5 of transfer. This in turn means that each thread must destroy aσ -polymer-thread. Since the number of σ -polymer-threads, |Tσ |, equals Thread-Count, all threads are destroyed (and the Tσ -polymer is empty) upon completionof thread-transfer.Bottleneck reactions. In each round, the halve(σ) function decreases the lengthof the H-polymer by a factor of 2, starting from n initially. Each decrement orincrement of the H-polymer includes a bottleneck reaction, so there are Θ(n) bot-tleneck reactions in total, over all rounds. The double function creates 2l thread-polymers in round l, for a total of Θ(n) thread-polymers over all rounds. Thetransfer function creates 2l threads in round l and similarly destroys 2l threads,and copies a polymer of length 2l , so again has Θ(n) bottleneck reactions over allrounds. The reactions in thread-transfer are not bottleneck reactions (ex-cept in round 1); we analyze these in the next chapter.29Chapter 5psCRN time complexity analysisand simulationWe follow the stochastic model of Soloveichik et al [18] for well-mixed, closedsystems with fixed volume V . To achieve a fixed volume, we make two new as-sumptions. First, that all our reactions are bimolecular, of the form r+ r′→ p+ p′,by adding a blank molecule B as a second reactant r′ to each reaction with onlya single reactant r, and as a second product p′ to each reaction with only a singleproduct p. Second, that every psCRN initializes with sufficient numbers of B suchthat the psCRN can reach its final committing state. In addition, we assume that allreactions have rate constant 1.Then when in configuration c, the propensity of reaction R : r+ r′→ p+ p′ isc(r)c(r′)/V if r 6= r′, and is (c(r)2 )/V if r = r′. Let ∆(c) be the sum of all reactionpropensities, when in configuration c. When a reaction occurs in configuration c,the probability that it is reaction R is the propensity of R divided by ∆(c), and theexpected time for a reaction is 1/∆(c). When the only applicable reaction is a bot-tleneck reaction, and the volume V is Θ(n2), the expected time for this bottleneckreaction is Θ(n2). Soloveichik et al [18] consider CRNs without polymers, but thesame stochastic model is used by Lakin et al. [16] and Qian et al. [17], where thereactants r or r′ (as well as the products) may be polymers.30Expected time complexity of Algorithm 1: Sequential-n2-psCRN. This psCRNhas n rounds, with Θ(n) instructions per round; for example, the copy of length nin each round has n inc instructions. So the total number of instructions executed,over all rounds is Θ(n2); moreover, there are Θ(n2) inc(σ ) instruction overall.The program’s instructions execute sequentially, that is, the ith instruction com-pletes before the (i+ 1)st instruction starts, so the total expected time is the sumof the expected times of the individual instructions. Each instruction involves aconstant number of reactions. Some instructions involve bottleneck reactions; forexample, the push reaction of the inc instruction is a bottleneck reaction. So anexecution of the program involves Θ(n2) bottleneck reactions. Each of these takesΘ(n2) time, so the overall expected time is Θ(n4).Expected time complexity of Algorithm 4: Threaded-n2blognc-psCRN. We notedearlier that Algorithm 4 has Θ(n) non-threaded instructions, and in fact Θ(n) bot-tleneck instructions. These take expected time Θ(n3) overall, since the time foreach is Θ(V ) =Θ(n2).Now, consider the threaded function, thread-transfer. In round l,1 ≤l ≤ blgnc, thread-transfer has 2l threads, and pushes n2l Y monomers onto 2l clone Y -polymers. Since each Y -push reaction is independent and is equallylikely to increment each of the 2l Y -polymers, the expected number of moleculesper polymer is n. Using Chernoff tail bounds, we can show that all polymers havelength in the range [n/2,2n] with all but negligibly small probability for large n.Chernoff Tail Bounds [10] Consider N independent trials, each with the samesuccess probability, such that the total number of successes is X and the expectednumber of successes is µ . Then for 0 < δ ≤ 1,(a) P[X ≤ (1−δ )µ]≤ exp(− δ 2µ2 ), and(b) P[X ≥ (1+δ )µ]≤ exp(− δ 2µ3 ).For a given polymer, let X be the number of increments to that polymer. Forthe upper tail bound, and setting δ = 1,µ = n we have:P(X ≥ (1+δ )µ)≤ e− δ22+δ µ = P(X ≥ 2n)≤ e− n3 . (5.1)31Since this bound holds for any given polymer, and since there are at most n poly-mers, the probability that some polymer has length greater than 2n is at most ne−n3 ,which is exponentially small in n.For the lower tail bound, setting δ = 1/2,µ = n we have:P(X ≤ (1−δ )µ)≤ eµδ 2/2 = P(X ≤ n/2)≤ e− n8 , (5.2)and so the probability that all polymers have length at least n/2 is at least 1−ne− n8 .Assume that all polymers have length in the range n/2 and 2n. During thefirst ≥ n/2 of the thread-transfer decrements in round l, the count of eachof the reactants is 2l: one program counter per thread and 2l polymers in total.So the expected time for each of these decrements is Θ(V/22l). Pessimistically,all decrements would happen to the same polymer, whose length could be as lit-tle as n/2 by our assumption above, and so there are 2l − 1 polymers and threadsavailable for the next decrements, 2l − 2 polymers and threads available for thenext Θ(n) decrements after that once a second polymer is depleted, and so on.On top of that, the polymers would be decremented in increasing order of length,as the later in a round a polymer is decremented, the longer each decrement isexpected to take. Therefore the expected time for a polymer to be decrementedwhen there are j polymers left is Θ(V n/ j2). Then the total expected time forround l is O(V n∑2lj=1(1/ j2)) = Θ(nV ). Multiplying by blgnc, the number ofrounds, and noting that V = Θ(n2), we have that the total expected time for thethread-transfer over all rounds is Θ(n3 logn), given that all polymers havelength in the range n/2 and 2n.When our assumption does not hold, in the worst case, all Y -increments andY -decrements happen to a single Y -polymer. In this case each Y -decrement be-comes a bottleneck reaction, taking time Θ(n2). As the number of Y -decrementsduring round l is Θ(n2l), the total number of Y -pops across all rounds must beΘ(n∑lgn−1i=0 2i) = Θ((2lgn − 1)n) = Θ(n2). This means the worst-case time forthread-transfer over all rounds becomes Θ(n4), the same as Algorithm 1.Finally, since the probability that our assumption does not hold is exponentiallysmall in n, we have that the total expected time for the thread-transfer overall rounds is Θ(n3 logn).32Simulator. To test the correctness of our protocols, we developed a custom CRNsimulator, implemented in Python, designed to support clone polymers. The simu-lator uses a slightly modified version of Gibson and Bruck’s next reaction method[12], which itself is an extension of Gillespie’s algorithm [13]. We redefine what asingle “species” is from the algorithm’s point of view, classifying all σ -polymersas one species, and track polymer lengths separately.Interestingly, simulation of our stable, sequential psCRN with input detectionfor Square appears to take only constant factor extra time compared to the com-mitting, fast-forwarded sequential psCRN (see both Figures 5.1 and 5.2). This isbecause each of the n error detection steps, and subsequent restart, is expected tohappen in O(n2) time, which is negligible compared to theΘ(n4) expected runningtime of the psCRN when fast-forwarding is assumed. On the other hand, compar-ing the runtimes of our sequential versus threaded algorithms shows that we gainsignificant speed from the parallelization (Figure 5.3). For all simulations, curvefitting was performed using the scipy package.33(a) (b)Figure 5.1: (a) Simulation of Algorithm 1: Sequential-n2-psCRN, with fast-forwarding and input n on X-polymer, n = 100. (b) Simulation ofAlgorithm 1: Sequential-n2-psCRN, with input detection and loading,n= 100. Each line in the plots shows the count of the outputs as a func-tion of the number of interactions, with the blue line being the count ofoutput species Y finally released into solution, while the red line showsthe size of the Yint polymer. By interaction, we mean a collision of twomolecular species in the system, which may or may not result in a reac-tion.(a) (b)Figure 5.2: Plotting simulation runtimes as a function of initial species countn for (a) Algorithm 1: Sequential-n2-psCRN starting with pre-loadedinput (b) Algorithm 1: Sequential-n2-psCRN with input detection andloading. Each blue marker plots a single simulation, showing the num-ber of interactions required for protocol completion as a function of theinitial count of input species X . In both cases, data was fitted with anequation of form cx4, where (a) c = 78.4, and (b) c = 99.2.34(a)(b)Figure 5.3: Comparing runtimes of Algorithm 1: Sequential-n2-psCRN (red)to Algorithm 4: Threaded-n2blognc-psCRN (blue) with fast-forwarding,where (a) n = 128, and (b) n = 2i, i ∈ [3,7]. We plot the average num-ber of interactions over five different simulations for each input countn. Sequential data was fitted with an equation of form cx4, c = 50.3.Threaded data was fitted with an equation of form cx3 logx, c = 267.35Chapter 6Conclusions and future workIn this work, we’ve expanded the computing model of stochastic chemical reactionnetworks with polymers, by considering inputs that are represented as monomersin solution, as well as clone polymers that facilitate distributed data structures andthreaded computation. We’ve shown that stable, error-free Turing-universal com-putation is possible in the monomer input model, by introducing an error-correctionscheme that takes advantage of the ability to check for empty polymers. We’ve il-lustrated how programming with clone polymers can provide speed-ups, comparedwith using leader polymers only, and how leader polymers can be used for syn-chronization purposes by CRNs with clone polymers.There are many interesting directions for future work. First, we have shownhow to use clone polymers to get a speed-up for the Square problem, but we havenot shown that such a speed-up is not possible without the use of clone polymers.Is it possible to show lower bounds on the time complexity of problems whenonly leader polymers are available? Or, could bottleneck reactions be reduced bya psCRN computing Square? Second, our faster psCRN for Square with clonepolymers still uses leader polymers for synchronization. Is the speed-up possibleeven without the use of leader polymers? More generally, how can synchroniza-tion be achieved in leaderless psCRNs? Are there faster psCRNs, with or withoutleader polymers? It would be very interesting to know what problems have sta-ble psCRNS that use no leaders, but can use clone polymers. Finally, it would bevaluable to have more realistic models of reaction propensities for psCRN mod-36els. Underlying the model in this paper is the assumption that the extensible endsof polymers are well-mixed in solution, along with other monomeric moleculesin the system. This seems implausible, in practice, since steric hindrance couldisolate a long polymer’s extensible end from other reactants, for example. In gen-eral, subunits that are close by on the polymer would also be close by in solution,which could increase the chances of certain interactions from happening, whiledecreasing the chances of other ones. Research into this area would at least en-tail building running simple simulations of how polymers interact with monomersand each other in solution, to build a new model that accounts for these aspectsof polymer behaviour. We expect that our contributions of correctness and thread-ing would remain relevant even with an updated model, with the focus being moresophisticated runtime analysis for psCRN programs.37Bibliography[1] D. Angluin, J. Aspnes, Z. Diamadi, M. J. Fischer, and R. Peralta.Computation in networks of passively mobile finite-state sensors.Distributed Computing, pages 235–253, Mar. 2006. → page 5[2] D. Angluin, J. Aspnes, and D. Eisenstat. Fast computation by populationprotocols with a leader. In Dolev S. (eds) Distributed Computing (DISC),Lecture Notes in Computer Science, volume 4167, pages 61–75. Springer,Berlin, Heidelberg, 2006. → page 12[3] D. Angluin, J. Aspnes, and D. Eisenstat. Stably computable predicates aresemilinear. In PODC ’06: Proceedings of the twenty-fifth annual ACMsymposium on Principles of distributed computing, pages 292–299, NewYork, NY, USA, 2006. ACM Press. → page 2[4] C. Bennett. Logical reversibility of computation. IBM journal of Researchand Development, 17(6):525–532, 1973. → page 4[5] C. Bennett. The thermodynamics of computation - a review. InternationalJournal of Theoretical Physics, 21(12):905–940, 1981. → page 4[6] L. Cardelli and G. Zavattaro. Turing universality of the biochemical groundform. Mathematical. Structures in Comp. Sci., 20(1):45–73, Feb. 2010.ISSN 0960-1295. → page 5[7] I. Chatzigiannakis, O. Michail, S. Nikolaou, A. Pavlogiannis, and P. G.Spirakis. Passively mobile communicating machines that use restrictedspace. In Proceedings of the 7th ACM ACM SIGACT/SIGMOBILEInternational Workshop on Foundations of Mobile Computing, FOMC ’11,pages 6–15, New York, NY, USA, 2011. ACM. ISBN 978-1-4503-0779-6.→ page 5[8] H.-L. Chen, R. Cummings, D. Doty, and D. Soloveichik. Speed faults incomputation by chemical reaction networks. In F. Kuhn, editor, Distributed38Computing, pages 16–30, Berlin, Heidelberg, 2014. Springer BerlinHeidelberg. ISBN 978-3-662-45174-8. → page 10[9] H.-L. Chen, D. Doty, and D. Soloveichik. Deterministic functioncomputation with chemical reaction networks. Natural Computing, 13(4):517–534, Dec 2014. → pages 1, 2[10] H. Chernoff. A measure of asymptotic efficiency for tests of a hypothesisbased on the sum of observations. The Annals of Mathematical Statistics,pages 493–507, 1952. URL http://www.jstor.org/stable/2236576. → page 31[11] R. Cummings, D. Doty, and D. Soloveichik. Probability 1 computation withchemical reaction networks. Natural Computing, 15(2):245–261, 2014. →pages 2, 4, 11, 12[12] M. A. Gibson and J. Bruck. Efficient exact stochastic simulation of chemicalsystems with many species and many channels. The Journal of PhysicalChemistry A, 104(9):1876–1889, 2000. → page 33[13] D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions.The Journal of Physical Chemistry, 81(25):2340–2361, 1977. → page 33[14] H. Jiang, M. Riedel, and K. Parhi. Synchronous sequential computation withmolecular reactions. In Proceedings of the 48th Design AutomationConference, DAC ’11, pages 836–841, New York, NY, USA, 2011. ACM.→ page 5[15] R. Johnson and E. Winfree. Verifying polymer reaction networks usingbisimulation. 2014. URLhttp://www.dna.caltech.edu/Papers/Polymers2014-VEMDP.pdf. → pages5, 6[16] M. R. Lakin and A. Phillips. Modelling, simulating and verifyingTuring-powerful strand displacement systems. In L. Cardelli and W. Shih,editors, DNA Computing and Molecular Programming, pages 130–144,Berlin, Heidelberg, 2011. Springer Berlin Heidelberg. → pages5, 6, 7, 22, 30[17] L. Qian, D. Soloveichik, and E. Winfree. Efficient Turing-universalcomputation with DNA polymers. In Proceedings of the 16th internationalconference on DNA computing and molecular programming, pages 123–140,Berlin, Heidelberg, 2010. Springer-Verlag. → pages 2, 4, 6, 7, 12, 3039[18] D. Soloveichik, M. Cook, E. Winfree, and J. Bruck. Computation with finitestochastic chemical reaction networks. Natural Computing, 7(4):615–633,Dec 2008. → pages 4, 10, 12, 30[19] D. Soloveichik, M. Cook, E. Winfree, and J. Bruck. Computation with finitestochastic chemical reaction networks. Natural Computing, 7(4):615–633,2008. → page 140
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Error-free stable computation with polymer-supplemented...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Error-free stable computation with polymer-supplemented chemical reaction networks Tai, Allison Yeu Yang 2019
pdf
Page Metadata
Item Metadata
Title | Error-free stable computation with polymer-supplemented chemical reaction networks |
Creator |
Tai, Allison Yeu Yang |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | When disallowing error, traditional chemical reaction networks are very limited in computational power: Angluin et al. and Chen et al. showed that only semilinear predicates and functions are stably computable by CRNs. Qian et al. and others have shown that polymer-supplemented chemical reaction networks are capable of Turing-universal computation. However, their model requires that inputs are loaded onto the polymers at protocol initialization, in contrast with the traditional convention that inputs are represented by counts of molecules in solution. Here, we show that polymer-supplemented chemical reaction networks can stably simulate Turing-universal computations even with solution-based inputs. However, such simulations use a unique ''leader'' polymer per input type and thus involve many slow bottleneck reactions. We further refine the polymer-supplemented chemical reaction network model to allow for clone polymers, that is, multiple functionally-identical copies of a polymer, and provide an illustrative example of how bottleneck reactions can be reduced in this new model. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-11-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0385114 |
URI | http://hdl.handle.net/2429/72211 |
Degree |
Master of Science - MSc |
Program |
Bioinformatics |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2020-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2020_may_tai_allisonyeuyang.pdf [ 410.37kB ]
- Metadata
- JSON: 24-1.0385114.json
- JSON-LD: 24-1.0385114-ld.json
- RDF/XML (Pretty): 24-1.0385114-rdf.xml
- RDF/JSON: 24-1.0385114-rdf.json
- Turtle: 24-1.0385114-turtle.txt
- N-Triples: 24-1.0385114-rdf-ntriples.txt
- Original Record: 24-1.0385114-source.json
- Full Text
- 24-1.0385114-fulltext.txt
- Citation
- 24-1.0385114.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0385114/manifest