Computing equilibrium concentrations of molecular complexes is generally analytically intractable and requires numerical approaches. In this work we focus on the polymer-monomer level, where indivisible molecules (monomers) combine to form complexes (polymers). Rather than employing free-energy parameters for each polymer, we focus on the athermic setting where all interactions preserve enthalpy. This setting aligns with the strongly bonded (domain-based) regime in DNA nanotechnology when strands can bind in different ways, but always with maximum overall bonding - and is consistent with the saturated configurations in the Thermodynamic Binding Networks (TBNs) model. Within this context, we develop an iterative algorithm for assigning polymer concentrations to satisfy detailed-balance, where on-target (desired) polymers are in high concentrations and off-target (undesired) polymers are in low. Even if not directly executed, our algorithm provides effective insights into upper bounds on concentration of off-target polymers, connecting combinatorial arguments about discrete configurations such as those in the TBN model to real-valued concentrations. We conclude with an application of our method to decreasing leak in DNA logic and signal propagation. Our results offer a new framework for design and verification of equilibrium concentrations when configurations are distinguished by entropic forces.

Discrete chemical reaction networks formalize the interactions of molecular species in a well-mixed solution as stochastic events. Given their basic mathematical and physical role, the computational power of chemical reaction networks has been widely studied in the molecular programming and distributed computing communities (e.g., in the language of population protocols). While for Turing-universal systems there is a universal measure of optimal information encoding based on Kolmogorov complexity, chemical reaction networks are not Turing universal unless error and unbounded molecular counts are permitted. Nonetheless, here we show that the optimal number of reactions to generate a specific count x ∈ ℕ with probability 1 is asymptotically equal to a "space-aware" version of the Kolmogorov complexity of x, defined as Ks(x) = minp {|p|/log|p| + log(space(U(p))) :U(p) = x}, where p is a program for universal Turing machine U. This version of Kolmogorov complexity incorporates not just the length of the shortest program for generating x, but also the space usage of that program. Probability 1 computation is captured by the standard notion of stable computation from distributed computing, but we limit our consideration to chemical reaction networks obeying a stronger constraint: they "know when they are done" in the sense that they produce a special species to indicate completion. As part of our results, we develop a module for encoding and unpacking any b bits of information via O(b/log b) reactions, which is information-theoretically optimal for incompressible information. Our work provides one answer to the question of how succinctly chemical self-organization can be encoded---in the sense of generating precise molecular counts of species as the desired state.

The field of chemical computation attempts to model computational behavior that arises when molecules, typically nucleic acids, are mixed together. By modeling this physical phenomenon at different levels of specificity, different operative computational behavior is observed. Thermodynamic binding networks (TBNs) is a highly abstracted model that focuses on which molecules are bound to each other in a "thermodynamically stable" sense. Stability is measured based only on how many bonds are formed and how many total complexes are in a configuration, without focusing on how molecules are binding or how they became bound. By defocusing on kinetic processes, TBNs attempt to naturally model the long-term behavior of a mixture (i.e., its thermodynamic equilibrium). We study the problem of signal amplification: detecting a small quantity of some molecule and amplifying its signal to something more easily detectable. This problem has natural applications such as disease diagnosis. By focusing on thermodynamically favored outcomes, we seek to design chemical systems that perform the task of signal amplification robustly without relying on kinetic pathways that can be error prone and require highly controlled conditions (e.g., PCR amplification). It might appear that a small change in concentrations can result in only small changes to the thermodynamic equilibrium of a molecular system. However, we show that it is possible to design a TBN that can "exponentially amplify" a signal represented by a single copy of a monomer called the analyte: this TBN has exactly one stable state before adding the analyte and exactly one stable state afterward, and those two states "look very different" from each other. In particular, their difference is exponential in the number of types of molecules and their sizes. The system can be programmed to any desired level of resilience to false positives and false negatives. To prove these results, we introduce new concepts to the TBN model, particularly the notions of a TBN’s entropy gap to describe how unlikely it is to be observed in an undesirable state, and feed-forward TBNs that have a strong upper bound on the number of polymers in a stable configuration. We also show a corresponding negative result: a doubly exponential upper bound, meaning that there is no TBN that can amplify a signal by an amount more than doubly exponential in the number and sizes of different molecules that comprise it. We leave as an open question to close this gap by either proving an exponential upper bound, or giving a construction with a doubly-exponential difference between the stable configurations before and after the analyte is added. Our work informs the fundamental question of how a thermodynamic equilibrium can change as a result of a small change to the system (adding a single molecule copy). While exponential amplification is traditionally viewed as inherently a non-equilibrium phenomenon, we find that in a strong sense exponential amplification can occur at thermodynamic equilibrium as well - where the "effect" (e.g., fluorescence) is exponential in types and complexity of the chemical components.

Understanding the algorithmic behaviors that are in principle realizable in a chemical system is necessary for a rigorous understanding of the design principles of biological regulatory networks. Further, advances in synthetic biology herald the time when we will be able to rationally engineer complex chemical systems and when idealized formal models will become blueprints for engineering. Coupled chemical interactions in a well-mixed solution are commonly formalized as chemical reaction networks (CRNs). However, despite the widespread use of CRNs in the natural sciences, the range of computational behaviors exhibited by CRNs is not well understood. Here, we study the following problem: What functions f : ℝ k → ℝ can be computed by a CRN, in which the CRN eventually produces the correct amount of the “output” molecule, no matter the rate at which reactions proceed? This captures a previously unexplored but very natural class of computations: For example, the reaction X 1 + X 2 → Y can be thought to compute the function y = min ( x 1 , x 2 ). Such a CRN is robust in the sense that it is correct whether its evolution is governed by the standard model of mass-action kinetics, alternatives such as Hill-function or Michaelis-Menten kinetics, or other arbitrary models of chemistry that respect the (fundamentally digital) stoichiometric constraints (what are the reactants and products?). We develop a reachability relation based on a broad notion of “what could happen” if reaction rates can vary arbitrarily over time. Using reachability, we define stable computation analogously to probability 1 computation in distributed computing and connect it with a seemingly stronger notion of rate-independent computation based on convergence in the limit t → ∞ under a wide class of generalized rate laws. Besides the direct mapping of a concentration to a nonnegative analog value, we also consider the “dual-rail representation” that can represent negative values as the difference of two concentrations and allows the composition of CRN modules. We prove that a function is rate-independently computable if and only if it is piecewise linear (with rational coefficients) and continuous (dual-rail representation), or non-negative with discontinuities occurring only when some inputs switch from zero to positive (direct representation). The many contexts where continuous piecewise linear functions are powerful targets for implementation, combined with the systematic construction we develop for computing these functions, demonstrate the potential of rate-independent chemical computation.

Life is built upon amazingly sophisticated molecular machines whose behavior combines mechanical and chemical action. Engineering of similarly complex nanoscale devices from first principles remains an as yet unrealized goal of bioengineering. In this paper we formalize a simple model of mechanical motion (mechanical linkages) combined with chemical bonding. The model has a natural implementation using DNA with double-stranded rigid links, and single-stranded flexible joints and binding sites. Surprisingly, we show that much of the complex behavior is preserved in an idealized topological model which considers solely the graph connectivity of the linkages. We show a number of artifacts including Boolean logic, catalysts, a fueled motor, and chemo-mechanical coupling, all of which can be understood and reasoned about in the topological model. The variety of achieved behaviors supports the use of topological chemical linkages in understanding and engineering complex molecular behaviors.

Engineering molecular systems that exhibit complex behavior requires the design of kinetic barriers. For example, an effective catalytic pathway must have a large barrier when the catalyst is absent. While programming such energy barriers seems to require knowledge of the specific molecular substrate, we develop a novel substrate-independent approach. We extend the recently-developed model known as thermodynamic binding networks, demonstrating programmable kinetic barriers that arise solely from the thermodynamic driving forces of bond formation and the configurational entropy of forming separate complexes. Our kinetic model makes relatively weak assumptions, which implies that energy barriers predicted by our model would exist in a wide variety of systems and conditions. We demonstrate that our model is robust by showing that several variations in its definition result in equivalent energy barriers. We apply this model to design catalytic systems with an arbitrarily large energy barrier to uncatalyzed reactions. Our results could yield robust amplifiers using DNA strand displacement, a popular technology for engineering synthetic reaction pathways, and suggest design strategies for preventing undesired kinetic behavior in a variety of molecular systems.

Biological regulatory networks depend upon chemical interactions to process information. Engineering such molecular computing systems is a major challenge for synthetic biology and related fields. The chemical reaction network (CRN) model idealizes chemical interactions, allowing rigorous reasoning about the computational power of chemical kinetics. Here we focus on function computation with CRNs, where we think of the initial concentrations of some species as the input and the equilibrium concentration of another species as the output. Specifically, we are concerned with CRNs that are rate-independent (the computation must be correct independent of the reaction rate law) and composable ( f°g can be computed by concatenating the CRNs computing f and g). Rate independence and composability are important engineering desiderata, permitting implementations that violate mass-action kinetics, or even "well-mixedness", and allowing the systematic construction of complex computation via modular design. We show that to construct composable rate-independent CRNs, it is necessary and sufficient to ensure that the output species of a module is not a reactant in any reaction within the module. We then exactly characterize the functions computable by such CRNs as superadditive, positive-continuous, and piecewise rational linear. Thus composability severely limits rate-independent computation unless more sophisticated input/output encodings are used.

Formal methods have enabled breakthroughs in many fields, such as in hardware verification, machine learning and biological systems. The key object of interest in systems biology, synthetic biology, and molecular programming is chemical reaction networks (CRNs) which formalizes coupled chemical reactions in a well-mixed solution. CRNs are pivotal for our understanding of biological regulatory and metabolic networks, as well as for programming engineered molecular behavior. Although it is clear that small CRNs are capable of complex dynamics and computational behavior, it remains difficult to explore the space of CRNs in search for desired functionality. We use Alloy, a tool for expressing structural constraints and behavior in software systems, to enumerate CRNs with declaratively specified properties. We show how this framework can enumerate CRNs with a variety of structural constraints including biologically motivated catalytic networks and metabolic networks, and seesaw networks motivated by DNA nanotechnology. We also use the framework to explore analog function computation in rate-independent CRNs. By computing the desired output value with stoichiometry rather than with reaction rates (in the sense that X → Y+Y computes multiplication by 2), such CRNs are completely robust to the choice of reaction rates or rate law. We find the smallest CRNs computing the max, minmax, abs and ReLU (rectified linear unit) functions in a natural subclass of rate-independent CRNs where rate-independence follows from structural network properties.

Embedding computation in molecular contexts incompatible with traditional electronics is expected to have wide ranging impact in synthetic biology, medicine, nanofabrication and other fields. A key remaining challenge lies in developing programming paradigms for molecular computation that are well-aligned with the underlying chemical hardware and do not attempt to shoehorn ill-fitting electronics paradigms. We discover a surprisingly tight connection between a popular class of neural networks (binary-weight ReLU aka BinaryConnect) and a class of coupled chemical reactions that are absolutely robust to reaction rates. The robustness of rate-independent chemical computation makes it a promising target for bioengineering implementation. We show how a BinaryConnect neural network trained in silico using well-founded deep learning optimization techniques, can be compiled to an equivalent chemical reaction network, providing a novel molecular programming paradigm. We illustrate such translation on the paradigmatic IRIS and MNIST datasets. Toward intended applications of chemical computation, we further use our method to generate a chemical reaction network that can discriminate between different virus types based on gene expression levels. Our work sets the stage for rich knowledge transfer between neural network and molecular programming communities.

Population protocols are a distributed computing model appropriate for describing massive numbers of agents with very limited computational power (finite automata in this paper), such as sensor networks or programmable chemical reaction networks in synthetic biology. A population protocol is said to require a leader if every valid initial configuration contains a single agent in a special "leader" state that helps to coordinate the computation. Although the class of predicates and functions computable with probability 1 (stable computation) is the same whether a leader is required or not (semilinear functions and predicates), it is not known whether a leader is necessary for fast computation. Due to the large number of agents n (synthetic molecular systems routinely have trillions of molecules), efficient population protocols are generally defined as those computing in polylogarithmic in n (parallel) time. We consider population protocols that start in leaderless initial configurations, and the computation is regarded finished when the population protocol reaches a configuration from which a different output is no longer reachable. In this setting we show that a wide class of functions and predicates computable by population protocols are not efficiently computable (they require at least linear time), nor are some linear functions even efficiently approximable. It requires at least linear time for a population protocol even to approximate division by a constant or subtraction (or any linear function with a coefficient outside of N), in the sense that for sufficiently small gamma > 0, the output of a sublinear time protocol can stabilize outside the interval f(m) (1 +/- gamma) on infinitely many inputs m. In a complementary positive result, we show that with a sufficiently large value of gamma, a population protocol can approximate any linear f with nonnegative rational coefficients, within approximation factor gamma, in O(log n) time. We also show that it requires linear time to exactly compute a wide range of semilinear functions (e.g., f(m)=m if m is even and 2m if m is odd) and predicates (e.g., parity, equality).

Understanding the algorithmic behaviors that are in principle realizable in a chemical system is necessary for a rigorous understanding of the design principles of biological regulatory networks. Further, advances in synthetic biology herald the time when we'll be able to rationally engineer complex chemical systems, and when idealized formal models will become blueprints for engineering.

The connection between self‐assembly and computation suggests that a shape can be considered the output of a self‐assembly “program,” a set of tiles that fit together to create a shape. It seems plausible that the size of the smallest self‐assembly program that builds a shape and the shape’s descriptional (Kolmogorov) complexity should be related. We show that when using a notion of a shape that is independent of scale, this is indeed so: in the tile assembly model, the minimal number of distinct tile types necessary to self‐assemble a shape, at some scale, can be bounded both above and below in terms of the shape’s Kolmogorov complexity. As part of the proof, we develop a universal constructor for this model of self‐assembly that can execute an arbitrary Turing machine program specifying how to grow a shape. Our result implies, somewhat counterintuitively, that self‐assembly of a scaled‐up version of a shape often requires fewer tile types. Furthermore, the independence of scale in self‐assembly theory appears to play the same crucial role as the independence of running time in the theory of computability. This leads to an elegant formulation of languages of shapes generated by self‐assembly. Considering functions from bit strings to shapes, we show that the running‐time complexity, with respect to Turing machines, is polynomially equivalent to the scale complexity of the same function implemented via self‐assembly by a finite set of tile types. Our results also hold for shapes defined by Wang tiling—where there is no sense of a self‐assembly process—except that here time complexity must be measured with respect to nondeterministic Turing machines.