A Molecular Implementation of the Least Mean Squares Estimator Christoph Zechner∗1 and Mustafa Khammash†2 1Department of Biosystems Science and Engineering, ETH Zurich, Basel, Switzerland Abstract – In order to function reliably, synthetic molec- Inthepresentwork,weconsideranotherpowerfulclass ular circuits require mechanisms that allow them to adapt of estimation schemes that are frequently used in adap- to environmental disturbances. Least mean squares (LMS) tive signal processing and control theory. These schemes schemes,suchascommonlyencounteredinsignalprocessing –termedleastmeanssquares(LMS)estimators[13,14]– and control, provide a powerful means to accomplish that iterativelycomputethesolutionofageneralleastsquares 7 goal. In this paper we show how the traditional LMS algo- problemthroughagradient-basedparametersearch. This 1 rithm can be implemented at the molecular level using only iterative structure allows a circuit to estimate unknown 0 2 a few elementary biomolecular reactions. We demonstrate quantities in an adaptive fashion by processing measure- our approach using several simulation studies and discuss its ments in realtime. In this paper we demonstrate how n a relevance to synthetic biology. LMS-type estimators can be realized using elementary J biomolecularreactions. Ourworkisrelatedto[10],where 3 theauthorshaveproposedaDNA-basedgradient-descent 1 Introduction scheme, which is able to learn static linear functions. In ] this work, however, we focus on dynamical and possibly N Engineered circuits in living cells often exhibit poor ro- stochasticbiochemicalsystemsthatshallbe identifiedby M bustness and substantial variations from one cell to the a molecular LMS estimator. . next [1, 2]. In extreme cases, they are found functional o in only a small fraction of cells in an isogenic popula- The remainder of the paper is structured as follows. i b tion, while others act unpredictably. A major cause for InSection2weintroduce the mathematicalnotationand q- such behavior is that the biochemical components that models required to describe molecular circuits. In Sec- [ constitute a circuit depend on factors in their molecular tion 3 we introduce the concept of LMS estimation and environment(or context)of the cell[3]. For instance, the present possible molecular realizations. In Section 4 we 1 rate at which a protein is expressed depends on the gene test the performance of the proposed circuits using sev- v 2 dosage or the number of available ribosomes and so on. eralsimulationstudiesanddiscusshowtheymaybeused 0 Recently, progress has been made in taking into ac- in practical applications. 6 count such environmental factors into the modeling and 0 design of molecular circuits [4, 5, 6, 7]. This can tremen- 2 Biochemical Reaction Networks 0 . dously improve the faithfulness of computational models 1 and in turn the predictability of rationally designed cir- Weconsiderwell-mixedmolecularreactionnetworkscom- 0 7 cuits. However, in practical scenarios, the origins and prising K molecular species Z = (Z1,...,ZK)T that in- 1 propertiesofpotentialdisturbancesarebarelyknownand teractwitheachotherthroughLreactionchannelsofthe : hardto anticipate duringdesigntime. Fromthis pointof form v i view, it seems barely realistic to tune a circuit in silico Z−h−i−(Z−(−t)⇀) Z+νi, (1) X suchthatitactsrobustlyunderallpossibleperturbations with i as the reaction index, Z(t) as the abundance of Z r that it may encounter in the real environment of a cell. a at time t, hi as a rate function determined by the law of A viable alternative is to employ adaptive design prin- mass-action and ν as the stoichiometric change associ- i ciples in which a circuit continuously senses and adjusts ated with reaction i. Throughout this paper, we follow itselftochangingenvironmentalconditions. Thisrequires the convention to denote molecular species as boldface molecular circuits that learn and make inference about symbols. Note that we use the same symbol also to refer their surroundings. A few attempts have been made re- to the circuit that those species constitute. cently to devise such circuits in the form of chemical re- We describe the time-evolution of Z as a continuous- action networks, for example to perform neural network time Markov chain (CTMC) that can take into account computations [8], to realize message passing inference [9] the inherent randomnessof biochemical reactions [15]. It or supervised learning [10]. Along these lines, we have canbeshown[16]thatthemolecularabundanceZ(t)sat- recentlyproposedamolecularimplementationofanopti- isfies a stochastic integral equation of the form mal filter that allowsone to estimate dynamically chang- ing noise signals [11]. This estimator was derived under L t Z(t)=Z + P h (Z(s))ds ν , (2) a Bayesianoptimality criterion by employing a Kushner- 0 i(cid:18)Z i (cid:19) i Xi=1 0 Stratonovich differential equation [12]. where P is an independent unit Poisson processes de- i ∗[email protected] scribing the firings of reaction i. Eq. (2) is commonly †[email protected] known as the random time change model. 1 Assuming the chemical species to be highly abundant, withα(n)asatuneablestep-size. Thechoiceofthelatter molecular fluctuations become negligible and eq. (2) can usuallyinvolvesatradeoffbetweentherateofconvergence be approximated by a deterministic rate equation of the andthesteadysateerrorofthescheme(i.e.,excesserror). form Since we consider continuous-time dynamical systems, L d we seekfor aninfinitesimal variantofthe LMSalgorithm Z˜(t)= h (Z˜(t))ν (3) dt i i [17], i.e., Xi=1 with Z˜(t)≈Z(t). We will make use of equations (2) and dθˆ(t)=−α(t) ∂ J(θˆ) , (7) (3) at a later in this manuscript to model stochastic and dt ∂θˆ (cid:12)(cid:12) (cid:12)θˆ(t) deterministic reaction networks, respectively. (cid:12) (cid:12) in which case α(t) can be understood as a rate at which the scheme adapts. 3 A Continuous-Time LMS Algorithm Suppose a circuit requires knowledge about certain envi- Target system ronmentalfactorsθ. Forexample, θ couldbe the number of phosphotases available to the circuit. However, these factors are typically not accessible directly by the circuit butonlyindirectlythroughavailableintermediatesY. In the example above, Y couldbe a proteinthat is targeted Y(t) e(t) + by that phosphatase, for instance. - The idea is now to use a second molecular circuit X X(t) Adaptive system thatisable toidentify the dynamicsofY throughasuit- able adaptation scheme. We assume here that X and Y areequivalentintheir structurebuthavedistinctparam- eters θˆ and θ, respectively. The goal of the adaptation scheme is to adjust the parametersθˆsuchas to minimize the discrepancy between the measured output Y(t) and the output of X (termed X(t)). The resulting optimal LMS scheme parameters θˆ∗ then represent an estimate of θ. Figure 1: Schematic illustration of an adaptive molecular cir- A suitable and analytically convenient metric to as- cuit. The goal is to construct a biomolecular circuit X that ad- sess the discrepancy between Y(t) and X(t) is the mean justsitsparameterstomimickthoseofatargetsystemY. While squared error those parameters are inaccessible by the circuit, it can measure the output Y(t) of the target system. This output is compared J(θˆ)=E (X(t)−Y(t))2 . (4) to the output of X to construct an error signal e(t), which is in (cid:2) (cid:3) turn used tofind theoptimal parametersof X. According to this measure, we seek for the set of param- eters that minimizes J(θˆ), i.e., θˆ∗ =argmin J(θˆ), (5) For the sake of simplicity, we make a few simplifica- θˆ tions related to X and Y before deriving the scheme. A moregeneralscenario,however,willbesubjectofafuture in which case θˆ∗ is referred to as the least mean squared manuscript. First, we assume that θ (and correspond- (LMS) estimator. The closed-form solution of this opti- inglyθˆ)containsonlyasingleparameterthatneedstobe mizationproblemcanbe foundin certainspecific scenar- identified from Y(t). Furthermore, we allow only Y(t) to ios,forinstanceinthecaseoflinearsystemdynamics[14]. becorruptedbymolecularnoise,andreactionsassociated Inmostscenarios,however,(5) isanalyticallyintractable with X(t) are assumed to evolve deterministically (e.g., and one has to minimize J(θˆ∗) numerically. Note that through appropriate rescaling of the associated reaction iterative schemes may be beneficial even when (5) is an- rates). alytically tractable, because it gives X the flexibility to Under these assumptions, the gradientof J(θˆ) is given readapt to changes in θ as will be shown later in this by manuscript. A common strategy to minimize J(θˆ) is to employ a ∂ ∂ J(θˆ)= E (X(t)−Y(t))2 gradient-basedmethod that, at eachiteration,moves the ∂θˆ ∂θˆ parameters θˆalong the direction of the steepest descent, (cid:2) (cid:3) ∂ giving rise to the well-known LMS algorithm. This algo- =E (X(t)−Y(t))2 (cid:20)∂θˆ (cid:21) rithm is usually used in a discrete-time scenario, for in- (8) ∂ stance when operated on a digital signal processing unit. =E 2(X(t)−Y(t)) X(t) Insuchcase,ateachtimeiterationn,thealgorithmwould (cid:20) ∂θˆ (cid:21) update the parameters using the relation ∂ =2E (X(t)−Y(t)) X(t) . (cid:20) ∂θˆ (cid:21) ∂ θˆn+1 =θˆn−α(n) J(θˆ) (6) ∂θˆ (cid:12)(cid:12) Since we assume that X(t) is deterministic, (8) further (cid:12)θˆn (cid:12) (cid:12) 2 simplifies to can be represented by ∂ ∂ f J(θˆ)=2E (X(t)−Y(t)) X(t) A+B−↽⇀−O ∂θˆ (cid:20) ∂θˆ (cid:21) b =2X(t) ∂ X(t)−2E[Y(t)] ∂ X(t) (9) O+C−λ−b−/⇀f D, ∂θˆ ∂θˆ =2X(t)S(t)−2E[Y(t)]S(t), assuming b,f >>λ. Specific implementations of the derived LMS adapta- with S(t) := ∂ X(t) as the sensitivity of X(t) with re- tion scheme will be given in the subsequent section. ∂θˆ spect to θˆ. Note that the particular formof this sensitiv- itydependsonXandhowθˆentersitsdynamics. Specific 4 Case Studies examples will be given later in Section 4. Inthissectionweprovideseveralnumericalandanalytical Now, plugging eq. (9) into (7) yields a dynamic equa- examples to demonstrate our molecular LMS estimation tion for the LMS estimator θˆ(t), i.e., framework. In all of the examples, we will consider a target process d θˆ(t)=−2α(t)X(t)S(t)+2α(t)E[Y(t)]S(t) ρ φ dt (10) ∅−⇀Y −⇀∅, (13) =−α˜(t)X(t)S(t)+α˜(t)E[Y(t)]S(t), with ρ and φ as the process parameters that we aim to identify using an LMS scheme. As indicated earlier, we with α˜(t)=2α(t). restrict ourselves to the case of a single unknown param- eter, meaning that we either have θ ={ρ} or θ ={φ}. 3.1.OnlineAdaptation. Eq. (10)providesthedesired continuous-time solution of the LMS problem. However, 4.1.Self-adjusting birth-rate. We first consider the in its current form it is not adaptive, meaning that it as- case where the birth-rate ρ is unknown to the circuit. sumes known (and fixed) statistics of Y (i.e., the output The goal is to construct a corresponding adaptive circuit mean E[Y(t)]). In practice, however, such statistics are X oftenunknownandtheymightalsovaryovertime. Using ∅−⇀θˆ X−⇀φ ∅, (14) LMS estimation, this problem can be bypassed by esti- matingE[Y(t)]online fromavailablemeasurementsY(t). whose birth-rate θˆ adapts to that of Y. To accomplish This way, the required statistics are extracted directly this, we require a molecular implementation of relation from data, which in turn allows the scheme to readapt (11). The firststepis to derivethe particularformofthe when θ changes. A common and simple approach is to sensitivity function S(t) from the dynamics of X. The approximate E[Y(t)] by the current value of Y(t) such latter is given by the rate equation (3) which in this case that reads d X(t)=θˆ−φX(t). (15) dθˆ(t)=−α˜(t)X(t)S(t)+α˜(t)Y(t)S(t) (11) dt dt Differentiating both sides of eq. (15) with respect to θˆ and we adopt this strategy also in the present work. A yields graphicaldepiction of the online LMS scheme is depicted d ∂ d X(t)= S(t)=1−φS(t). (16) in Fig. 1. dt∂θˆ dt Fortunately, this equation is already in the form of a 3.2.Molecular Implementations. The goalisnowto valid rate equation. In particular, it describes the time- synthesize eq. (11) using biochemical reactions. How- evolution of a birth-death process ever, in its present form eq. (11) is incompatible with mass-actionratelawsbecause it containsa negative(i.e., ∅−⇀1 S−⇀φ ∅. (17) degradation) flux that does not depend on the current value of θˆ(t). To account for this, we choose the adap- Inconjunctionwith(11),the overalladaptivesystemcan tation rate to be proportional to θˆ(t), i.e., α˜(t) := λθˆ(t) be implemented through reactions and thus, ∅−⇀θ Y dθˆ(t)=−λθˆ(t)X(t)S(t)+λθˆ(t)Y(t)S(t). (12) Y −⇀φ ∅ dt θˆ −⇀1 θˆ+X While (12) is now in principle compatible with mass- φ actionkinetics,itinvolvestrimolecularreactions,thatare X−⇀∅ (18) hardormaybeimpossibletorealizeinpractice. However, ∅−⇀1 S the trimolecular reaction can be composed from two bi- φ molecular reactions with appropriately chosen rate con- S−⇀∅ stants. For example, the reaction S+θˆ+X−⇀λ S+X A+B+C−⇀λ D S+θˆ+Y −⇀λ S+2θˆ+Y. 3 We first studied the adaptation performance of (18) 0.3 as a function of the tuning parameter λ in an idealized noise-free scenario. In this case, the network from (18) is 0.2 described by the rate equations 0.1 d Y(t)=θ−φY(t) (19) dt 0 d X(t)=θˆ(t)−φX(t) (20) 0 20 40 60 80 100 dt Time d θˆ(t)=−λθˆ(t)X(t)S(t)+λθˆ(t)Y(t)S(t) (21) true value dt d Figure 2: Convergence ofestimated birth-rateasafunction of S(t)=1−φS(t). (22) λ. Theadaptivecircuitfrom(18)wassimulatedusingparameters dt θ =0.2 and φ= 0.01. Small λ lead to slow convergence, while This equation was simulated for different values of λ as λ > φ3/(4θ) = 1.25e −6 cause overshooting and oscillatory depicted in Fig. 2. The results indicate a tradeoff that is behavior. associated with the choice of λ: too small λ lead to slow convergence of the scheme, while too large λ cause the 40 0.4 adaptation scheme to “overshoot” the target value and exhibit oscillations. 30 0.3 In order to analyze the convergence properties of (18), X(t)20 0.2 we performed a local stability analysis of (22) based 10 0.1 on linearization. Noting that Y(t) and S(t) evolve au- 0 0 tonomously,wecanreplacethosevariablesbytheirsteady 0 200 400 600 0 200 400 600 Time Time state values Y∞ =θ/φ and S∞ =1/φ, respectively. This True Estimated yields the reduced system Figure 3: Algorithm performance in the presence of molecular d X(t)=θˆ(t)−φX(t) (23) noise. The adaptive circuit from (18) was simulated using the dt stochastic simulation algorithm to account for molecular noise. dθˆ(t)=−λθˆ(t)X(t)+ λθθˆ(t) (24) We assume that the target value θ changes spontaneously at dt φ φ2 certain time points tocheck thecircuit’s ability to readapt. The parameters used for the simulations were θ =0.2, φ=0.01 and which has equilibrium points X∞1 = θˆ∞1 = 0 and X∞2 = λ=1e−6. θ/φ and θˆ2 =θ. The respective Jaccobians are given by ∞ −φ 1 −φ 1 A1 =(cid:18) 0 λφθ2(cid:19) A2 =(cid:18)−θφλ 0(cid:19). (25) vioarluoefSθˆ(∞t),=m1e/aφn)indgoetshantotθcwhiallngreemthaienaassymthpetootnilcybsethaabvle- The eigenvalue λθ/φ2 of A1 is always positive meaning equilibriumpoint. However,theadaptationlawsimplifies that for initial conditions θˆ(0) > 0, the system will not to converge to that equilibrium point. From A2 we find d λ λ θˆ(t)=− θˆ(t)X(t)+ θˆ(t)Y(t) eigenvalues dt φ φ (27) λ φ4−4θλφ+φ2 − φ4−4θλφ+φ2 = θˆ(t)(Y(t)−X(t)), µ1 =−p 2φ µ2 =− p 2φ . φ (26) whichisstructurallyequivalenttoaspecificcontrolmotif For any λ > 0, the real part of both eigenvalues is that has been studied previously [19, 20]. In particular, negative, meaning that θ is the only stable equilibrium it was shown to act as an integral controlcircuit exhibit- of the adaptation scheme. However, we find that for ing robust perfect adaptation [21]. This points out the λ > φ3/(4θ), the system will have complex eigenvalues, potential use of our adaptive estimation framework for indicatingoscillatorybehavior. Thiscanalsobeseenfrom studying robustness in biological networks. in Fig. 2, especially for the case λ=3e−6. Wenextanalyzedthecircuit’sperformanceinthepres- 4.2.Self-adjusting death-rate. We further show how ence of molecular fluctuations and spontaneous changes the LMS adaption can be used to identify the target cir- inθ. We usedthe modelfromeq. (2)andstochasticsim- cuit’s death rate θ = {φ}. The corresponding sensitivity ulations [18] to simulate the reaction network from (18). S(t) can be shown to satisfy The results from Fig. 3 show that the birth-rate θ and in turn the system output Y(t) is accurately trackedby the d molecularLMSscheme. Adetailedandquantitativeerror S(t)=−X(t)−θˆS(t). (28) dt analysis in the presence of molecular fluctuations will be subject of future work. There are two issues associated with the above equation. Remark: We want to point out another interesting First, it depends on the tuning parameter θˆ, which will propertyofthis particularLMSestimator. Replacingthe change over time due to the adaptation scheme. Corre- sensitivityS(t)byapositiveconstant(e.g.,itsstationary spondingly, the value of S(t) will be different from the 4 actual sensitivity ∂ X(t). While this will have an im- found that for any λ > 0, only the non-zero equilibrium ∂θˆ pact on the convergence rate of the LMS scheme, it does point is stable. For any λ > θ4/(4ρ2), the adaptation not affect its steady state behavior (see analyticalresults exhibit oscillatory behavior, which should be taken into below). The second problem is that S(t) is incompati- consideration when designing this circuit. These results ble with mass-action rate laws due to the negative de- are confirmed in Fig. 4, for which we simulated (35) for pendency on X(t). In order to address this problem, we three different values of λ. − consider the equation for S (t) = −S(t), which is given by 0.02 d S−(t)=X(t)−θˆS−(t) (29) dt and correspondingly use the LMS update rule 0.01 d θˆ(t)=−λθˆ(t)X(t)S(t)+λθˆ(t)Y(t)S(t) dt (30) =λθˆ(t)X(t)S−(t)−λθˆ(t)Y(t)S−(t). 0 0 50 100 Time Overall, the adaptive circuit is given by the reactions true value ρ ∅−⇀Y Figure 4: Convergenceofestimateddeath-rateasafunctionof Y−⇀θ ∅ λ. Theadaptivecircuit from (31)wassimulated forθ=0.1and φ=0.01anddifferentvaluesofλ. Forλ>θ4/(4ρ2)=2.5e−7, ρ ∅−⇀X theadaptation scheme exhibits oscillatory behavior. θˆ+X−⇀1 θˆ (31) X−⇀1 X+S θˆ+S−⇀1 ∅ 5 CONCLUSIONS S+θˆ+Y−⇀λ S+Y S+θˆ+X−⇀λ S+2θˆ+X. LMSschemesprovideapowerfulandversatileframework for adaptive estimation. In this work we have shown Similar to the previous section, we performed simula- how simple LMS-type estimators that usually run on a tions to check the adaptation performance of the circuit computer can be implemented biochemically for the pur- as a function of λ under idealized noise-free conditions. pose of synthetic biology. Such algorithms would allow This allowsus to describe the adaptivecircuitby the dif- a circuit to make inference about its environment and ferential equations facilitate adaptive behavior. We have shown by simula- tion that the LMS circuit is able to accurately estimate d Y(t)=ρ−θY(t) (32) and track unknown parameters of a birth-death process. dt We arecurrently extending the proposedschemeto more d X(t)=ρ−θˆ(t)X(t) (33) general scenarios when multiple parameters are to be es- dt timated simultaneously and when both X and Y are ar- d θˆ(t)=−λθˆ(t)X(t)S(t)+λθˆ(t)Y(t)S(t) (34) bitrary, possibly multivariate stochastic circuits. dt We anticipate several important potential applications d S(t)=−X(t)−θˆ(t)S(t). (35) of the presented framework. In [11], optimal filters (such dt as the Kalman filter [22]) were employed to design noise- cancelling synthetic circuits. To this end, the LMS ap- We again performed a local stability analysis of the proachprovidesanattractivealternativetooptimalfilter- differentialequationstoinvestigatetheconvergenceofthe ingdue toits genericandsimple structure. Inthefuture, circuit. Inthiscase,thesensitivityS(t)iscoupledtoX(t) and θˆ(t) and thus, has to be included in the dynamic wewillextendtheapproachtononlinearandmultivariate modelsandprovideanin-depthanalysisofitsproperties. analysis. Aftereliminatingtheequationcorrespondingto Y(t), we obtain a three-dimensional system d ACKNOWLEDGMENTS. Thisprojectwasfinanced X(t)=ρ−θˆ(t)X(t) (36) dt with a grantfrom the Swiss SystemsX.ch initiative, eval- d λρ uated by the Swiss National Science Foundation. θˆ(t)=−λθˆ(t)X(t)S(t)+ θˆ(t)S(t) (37) dt θ d S(t)=−X(t)−θˆ(t)S(t), (38) COPYRIGHT INFORMATION. This article was dt presented and published at the 2016 IEEE 55th Confer- whichhas equilibria at the originandat the point X∞3 = ence on Decision and Control (CDC) in Las Vegas. In ρ/θ,θˆ3 =θ and S3 =−1/θ2. For compactness, we skip the original version, Figure 1 was unreferenced and θ in ∞ ∞ explicit expressions of the respective Jaccobian matrices eq. (28) shouldhave been θˆ. Inthe presentversion,both andeigenvalues. However,asinthepreviousexample,we issues have been corrected. 5 References [11] Christoph Zechner, Georg Seelig, Marc Rullan, and Mustafa Khammash.Molecularcircuitsfordynamicnoisefiltering.Pro- [1] Timothy S Gardner, Charles R Cantor, and James J Collins. ceedings of the National Academy of Sciences, 113(17):4729– ConstructionofagenetictoggleswitchinEscherichiacoli.Na- 4734,2016. ture,403(6767):339–342, 2000. [12] Harold J Kushner. On the differential equations satisfied by [2] Michael B Elowitz and Stanislas Leibler. A synthetic conditional probablitity densities of Markov processes, with oscillatory network of transcriptional regulators. Nature, applications. SIAM Journal on Control and Optimization, 403(6767):335–338, 2000. 2(1):106, 1964. [3] StefanoCardinaleandAdamPaulArkin.Contextualizingcon- [13] Simon Haykin. Adaptive Filter Theory (3rd Ed.). Prentice- textforsyntheticbiology–identifyingcausesoffailureofsyn- Hall,Inc.,UpperSaddleRiver,NJ,USA,1996. theticbiologicalsystems.BiotechnologyJournal,7(7):856–866, [14] S.M.Kay. FundamentalsofStatisticalSignalProcessing: Es- 2012. timation Theory. PrenticeHall,EnglewoodCliffs,NJ,1993. [4] Artemis Llamosi, Andres M. Gonzalez-Vargas, Cristian Ver- [15] Nicolaas Godfried Van Kampen. Stochastic processes in sari, Eugenio Cinquemani, Giancarlo Ferrari-Trecate, Pascal physics and chemistry,volume1. Elsevier,1992. Hersen,andGregoryBatt. Whatpopulationrevealsaboutin- dividualcellidentity: Single-cellparameterestimationofmod- [16] D. F. Anderson and T. G. Kurtz. Continuous time Markov elsofgeneexpressioninyeast. PLoSComputBiol,12(2):1–18, chain models for chemical reaction networks. In Design and 022016. AnalysisofBiomolecularCircuits,pages3–42.Springer,2011. [5] Christoph Zechner and Heinz Koeppl. Uncoupled analysis of [17] S. Karni and G. Zeng. The analysis of the continuous-time stochasticreactionnetworksinfluctuatingenvironments. Plos lmsalgorithm. IEEE Transactions on Acoustics, Speech, and Computational Biology,10(12):e1003942, 2014. Signal Processing,37(4):595–597, 1989. [6] Tina Toni and Bruce Tidor. Combined model of intrin- [18] DanielT.Gillespie. Stochasticsimulationofchemicalkinetics. sic and extrinsic variabilityfor computational network design Annual Review of Physical Chemistry,58(1):35–55, 2007. with application to synthetic biology. PLoS Comput Biol, [19] Oren Shoval, Lea Goentoro, Yuval Hart, Avi Mayo, Eduardo 9(3):e1002960, 2013. Sontag, andUri Alon. Fold-change detection and scalar sym- [7] J. Hasenauer, S. Waldherr, M. Doszczak, N. Radde, metry of sensory input fields. Proceedings of the National P. Scheurich, and F. Allgower. Identification of models of Academy of Sciences,107(36):15995–16000, 2010. heterogeneouscellpopulationsfrompopulationsnapshotdata. [20] C. Briat, C. Zechner, and M. Khammash. Design of a syn- BMC Bioinformatics,12(1):125, 2011. thetic integral feedback circuit: dynamic analysis and DNA [8] LuluQian,ErikWinfree,andJehoshuaBruck.Neuralnetwork implementation. ACSSynthetic Biology (in press),2016. computationwithDNAstranddisplacementcascades. Nature, [21] Tau-Mu Yi, Yun Huang, Melvin I. Simon, and John Doyle. 475(7356):368–372, 2011. Robustperfectadaptationinbacterialchemotaxisthroughin- [9] Nils E Napp and Ryan P Adams. Message passing inference tegralfeedback control. Proceedings of the National Academy withchemicalreactionnetworks.InC.J.C.Burges,L.Bottou, of Sciences,97(9):4649–4653, 2000. M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, [22] RudolphEKalmanandRichardSBucy. Newresultsinlinear AdvancesinNeuralInformationProcessingSystems26,pages filteringandpredictiontheory. JournalofFluidsEngineering, 2247–2255. CurranAssociates,Inc.,2013. 83(1):95–108, 1961. [10] Matthew R.LakinandDarkoStefanovic. Supervisedlearning in adaptive DNA strand displacement networks. ACS Syn- thetic Biology,5(8):885–897, 2016. 6