{'fraction_non_alphanumeric': 0.03214294326155627, 'fraction_numerical': 0.00863253840032794, 'mean_word_length': 4.507569721115538, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chemical Boltzmann Machines

How smart can a micron-sized bag of chemicals be? How can an artificial or real cell make inferences about its environment? From which kinds of probability distributions can chemical reaction networks sample? We begin tackling these questions by showing four ways in which a stochastic chemical reaction network can implement a Boltzmann machine, a stochastic neural network model that can generate a wide range of probability distributions and compute conditional probabilities. The resulting models, and the associated theorems, provide a road map for constructing chemical reaction networks that exploit their native stochasticity as a computational resource. Finally, to show the potential of our models, we simulate a chemical Boltzmann machine to classify and generate MNIST digits in-silico.


Introduction
To carry out complex tasks such as finding and exploiting food sources, avoiding toxins and predators, and transitioning through critical life-cycle stages, singlecelled organisms and future cell-like artificial systems must make sensible decisions based on information about their environment [1,2]. The small volumes of cells makes this enterprise inherently probabilistic: environmental signals and the biochemical networks within the cell are noisy, due to the stochasticity inherent in the interactions of small, diffusing molecules [3,4,5]. The small volumes of cells raises questions not only about how stochasticity influences circuit function, but also about how much computational sophistication can be packed into the limited available space.
Perhaps surprisingly, neural network models provide an attractive architecture for the types of computation, inference, and information processing that cells must do. Neural networks can perform deterministic computation using circuits that are smaller and faster than boolean circuits composed of AND, OR, and NOT gates [6], can robustly perform tasks such as associative recall [7], and naturally perform Bayesian inference [8]. Furthermore, the structure of biochemical networks, such as signal transduction cascades [9,1,10] and genetic regulatory networks [11,12,13,14,15], can map surprisingly well onto neural network architectures. Chemical implementations of neural networks and related machine In a micron-scale environment, molecular counts are low and a real (or synthetic) cell will have to respond to internal and environmental cues. Probabilistic inference using chemical Boltzmann machines provides a framework for how this may be achieved.
learning models have also been proposed [16,17,18,19,20], and limited examples demonstrated [21,22,23,24], for synthetic biochemical systems. Previous work on biochemical neural networks and biochemical inference invoked models based on continuous concentrations of species representing neural activities. Such models are limited in their ability to address questions of biochemical computation in small volumes, where discrete stochastic chemical reaction network models must be used to account for the low molecular counts. The nature of biochemical computation changes qualitatively in this context. In particular, stochasticity has been widely studied in genetic regulatory networks [25], signaling cascades [26], population level bet hedging in bacteria [27], and other areas [28,29] -where the stochasticity is usually seen as a challenge limiting correct function, but is occasionally also viewed as a useful resource [30]. Our work falls squarely in the latter camp: we attempt to exploit the intrinsic stochastic fluctuations of a formal chemical reaction network (CRN) to build natively stochastic samplers by implementing a stochastic neural network. This links to efforts to build natively stochastic hardware for Bayesian inference [31,32] and to the substantial literature attempting to model, and find evidence for, stochastic neural systems capable of Bayesian inference [33,34].
Specifically, we propose CRNs that implement Boltzmann machines (BMs), a flexible class of Markov random fields capable of generating diverse distributions and for which conditioning on data has straightforward physical interpretations [8,35]. BMs are an established model of probabilistic neural networks due to their analytic tractability and identity with spin systems in statistical physics [36] and Hopfield networks in computer science [7]. These networks have been studied extensively and used in a wide range of applications including image classification [37] and video generation [38]. We prove that CRNs can implement BMs and that this is possible using detailed balanced CRNs. Moreover, we show that many of the attractive features of BMs can be applied to our CRN constructions such as inference, a straightforward learning rule and scalability to real-world data sets. We thereby introduce the idea of a chemical Boltzmann machine (CBM), a chemical system capable of exactly or approximately performing inference using a stochastically sampled high-dimensional state space, and explore some of its possible forms.

Boltzmann Machines (BMs):
Boltzmann machines are a class of binary stochastic neural networks, meaning that each node randomly switches between the values 0 and 1 according to a specified distribution. They are widely used for unsupervised machine learning because they can compactly represent and manipulate high-dimensional probability distributions. Boltzmann machines provide a flexible machine learning architecture because, as generative models, they can be used for a diverse set of tasks including data classification, data generation, and data reconstruction. Additionally, the simplicity of the model makes them analytically tractable. The use of hidden units (described below) allows Boltzmann machines to represent high order correlations in data. Together, these features make Boltzmann machines an excellent starting point for implementing stochastic chemical computers.
Fix a positive integer N ∈ Z >0 . An N -node Boltzmann machine (BM) is specified by a quadratic energy function E : where θ i ∈ R is the bias of node i, and w ij = w ji ∈ R is the weight of the unordered pair (i, j) of nodes, with w ii = 0. One may specify a BM architecture, or graph topology, by choosing additional weights w ij that are to be set to 0. In this paper, we will use N (i) = {j s.t. w ij = 0} to denote the neighborhood of i. From a physical point of view, we are implicitly using temperature units k B T for energy, which we will continue to do throughout this paper. A BM describes a distribution P (x) over state vectors x = (x 1 , . . . , x N ) ∈ {0, 1} N , Nodes of a BM are often partitioned into sets V and H of visible and hidden, respectively. Nodes in V represent data, and auxiliary nodes in H allow more complex distributions to be represented in the visible nodes. An implementation of a BM is a stationary stochastic process that samples from this distribution in the steady state. A BM can be implemented in silico using the Gibbs sampling algorithm [39], which induces a discrete time Markov chain (DTMC) on the state space {0, 1} N in such a way that the stationary distribution of this Markov chain corresponds to the distribution P (x). In each round, one node i ∈ {1, . . . , n} is chosen at random for update. For any two adjacent configurations x and x which differ only at node i -i.e., x i = 1−x i and x j = x j for all j = i -we set the transition probabilities T x→x of the DTMC so that Any function T x→x can be chosen so long as (3) is satisfied. One common choice A Boltzmann machine is also an inference engine. One can do inference on P (x) by conditioning on the values of a subset of the nodes. Suppose nonempty node subsets U and Y form a partition of the nodes {1, 2, . . . , N }, and fix u ∈ {0, 1} U . To obtain samples from P (y | u) where y ∈ {0, 1} Y , one clamps every node i ∈ U to the state u i while running Gibbs sampling, i.e., one does not allow these nodes to update. Clamping nodes to an input state is the same as specifying the input data for a statistical model. Steady state samples y ∈ {0, 1} Y of this procedure are draws from the distribution P (y | u).
Boltzmann machines can be used to learn a generative model from unlabeled data. After specifying the architecture, one then proceeds to find the weights, w ij , and biases, θ i , that maximize the likelihood of the observed data according to the model, using gradient descent from a random initial parameterization. This reduces to a very simple two-phase learning rule where weights on active edges are strengthened in a "wake phase" during which the network is clamped to observed data and are weakened in a "sleep phase" during which the network runs free [8,35]. Given a target distribution Q(x), this gradient descent corresponds to calculating the gradient of the Kullback-Leibler divergence from P to , with respect to the parameters θ i and w ij : (4) where · P and · Q denote expected values with respect to the distributions P and Q respectively. When hidden units are present, the distribution Q (which is defined on visible units only) is extended to hidden units based on clamping the visible units according to Q and using the conditional distribution P (y|u) to determine the hidden units.

Chemical Reaction Networks (CRNs):
Fix a finite set S = (S 1 , S 2 , . . . , S M ) of M species. A reaction r is a formal chemical equation abbreviated as r = µ r → ν r where µ r , ν r ∈ N S are the stoichiometric coefficient vectors for the reactant and product species respectively, and N = Z ≥0 . A reaction rate constant, k r ∈ R >0 , is associated with each reaction. In this paper, we define a chemical reaction network (CRN) as a triple C = (S, R, k) where S is a finite set of species, and R is a set of reactions, and k is the associated set of reaction rate constants.
We will denote chemical species by capital letters, and their counts by lower case letters; e.g., s 1 denotes the number of species S 1 . Thus the state of a stochastic CRN is described by a vector on a discrete lattice, s = (s 1 , s 2 . . . s M ) ∈ N S . The dynamics of a stochastic CRN are as follows [40]. The probability that a given reaction occurs per unit time, called its propensity, is given by Each time a reaction fires, state s changes to state s + ∆ r , where ∆ r = ν r − µ r is called the reaction vector, and the propensity of each reaction may change.
Viewed from a state space perspective, stochastic CRNs are continuous time Markov chains (CTMCs) with transition rates and thus their dynamics follow dP (s, t) dt = s =s where P (s, t) is the probability of a state with counts s at time t. Equivalently, they are governed by the chemical master equation, A stationary distribution π(s) may be found by solving dP (s,t) dt = 0 simultaneously for all s; in general, it need not be unique, and even may not exist. Given an initial state s 0 , π(s) = P (s, ∞) is unique if it exists. For that initial state, the reachability class Ω s0 ⊆ N M is the maximal subset of the integer lattice accessible to the CRN via some sequence of reactions in R. We will specify a CRN and a reachability class given an initial state as a shorthand for specifying a CRN and a set of initial states with identical reachability classes.

Detailed Balanced Chemical Reaction Networks:
A CTMC is said to satisfy detailed balance if there exists a well-defined function of state s, E(s) ∈ R, such that for every pair of states s and s , the transition rates R s→s and R s →s are either both zero or If the state space Ω is connected and the partition function Z = s∈Ω e −E(s) is finite, then the steady state distribution π(s) = 1 Z e −E(s) is unique, and the net flux between all states is zero in that steady state.
There is a related but distinct notion of detailed balance for a CRN. An equilibrium chemical system is constrained by physics to obey detailed balance at the level of each reaction. In particular, for a dilute equilibrium system, each species S i ∈ S has an energy G[S i ] ∈ R that relates to its intrinsic stability, and where ∆ i r + is the ith component of ∆ r + = ν r + − µ r + , and we have defined Any CRN for which there exists a function G satisfying (11) is called a detailed balanced CRN. To see that the CTMC for a detailed balanced CRN also itself satisfies detailed balance in the sense of (10), let s = s + ∆ r+ and note that (6) and (11) imply that for all reactions r + . Here, G(s) is a well-defined function of state s (the free energy) that can play the role of E in (10). If there are multiple reactions that bring s to s , they all satisfy (12), and therefore the ratio R s→s /R s →s satisfies (10) and the CTMC satisfies detailed balance. It is possible to consider non-equilibrium CRNs that violate (11). Such systems must be coupled to implicit reservoirs of fuel molecules that can drive the species of interest into a non-equilibrium steady state [41,42,43]. Usually -but not always [44,45] -the resultant Markov chain violates detailed balance. In Section 3.1, we shall consider a system that exhibits detailed balance at the level of the Markov chain, but is necessarily non-equilibrium and violates detailed balance at the detailed chemical level.
Given an initial condition s 0 , a detailed balanced CRN will be confined to a single reachability class Ω s0 . Moreover, from the form of G(s), the stationary distribution π(s) on Ω s0 of any detailed balanced CRN exists, is unique, and is a product of Poisson distributions restricted to the reachability class [46], with the partition function Z = s ∈Ωs 0 e −G(s ) dependent on the reachability class. Note that this implies that the partition function is always finite, even for an infinite reachability class.

Clamping and Conditioning with Detailed Balanced CRNs:
In a Boltzmann machine that has been trained to generate a desired probability distribution when run, inference can be performed by freezing, also known as clamping, the values of known variables, and running the rest of the network to obtain a sample; this turns out to exactly generate the conditional probability. A similar result holds for a subclass of detailed balanced CRNs that generate a distribution, for an appropriate notion of clamping in a CRN. Imagine a "demon" that, whenever a reaction results in a change in the counts of one of the clamped species, will instantaneously change it back to its previous value. If every reaction is such that either no clamped species change, or else every species that changes is clamped, then the demon is effectively simply "turning off" those reactions. More precisely, consider a CRN, C = (S, R, k). We will partition the species into two disjoint groups Y = S f ree and U = S clamped , where S f ree will be allowed to vary and S clamped will be held fixed. We will define free reactions, R f ree , as reactions which result in neither a net production nor a net consumption of any clamped species. Similarly, clamped reactions, R clamped are reactions which change the counts of any clamped species. The clamped CRN will be denoted C| U =u to indicate the the species U i ∈ U have been clamped to the values u i . The clamped CRN is defined by C| U =u = (S, R f ree , k f ree ), that is, the entire set of species along with the reduced set of reactions and their rate constants. In the clamped CRN it is apparent that the clamped species will not change from their initial conditions because no reaction in R f ree can change their count. However, these clamped species may still affect the free species catalytically. If the removed reactions, R clamped , never change counts of non-clamped species, then C| U =u is equivalent to the action of the "demon" imagined above.
We use equation 13 to prove that clamping a detailed balanced CRN is equivalent to calculating a conditional distribution, and to show when the conditional distributions of a detailed balanced CRN will be independent. Together, these theorems provide guidelines for devising detailed balanced CRNs with interesting (non-independent) conditional distributions and for obtaining samples from these distributions via clamping.
We will need one more definition. Let C be a detailed balanced CRN with reachability class Ω s0 for some initial condition s 0 = (u 0 , y 0 ). Let Γ s0 be the reachability class of the clamped CRN C| U =u0 with species U clamped to u 0 and species Y free. We say clamping preserves reachability if where In other words, clamping preserves reachability if, whenever a state s = (u 0 , y) is reachable from s 0 by any path in C, then it is also reachable from s 0 by some path in C| U =u0 that never changes u. Theorem 1. Consider a detailed balanced CRN C = (S, R, k) with reachability class Ω s0 from initial state s 0 . Partition the species into two disjoint sets Let the projection of s 0 onto U and Y be u 0 and y 0 . The conditional distribution P (y | u) implied by the stationary distribution π of C is equivalent to the stationary distribution of a clamped CRN, C| U =u starting from initial state s 0 with u 0 = u, provided that clamping preserves reachability.
Then, the conditional probability distribution of the unclamped CRN is given by Simply removing pairs of forwards and backwards reactions will preserve detailed balance for unaffected transitions, and hence the clamped system remains a detailed balanced CRN with the same free energy function. We then readily see that the clamped CRN's stationary distribution, π c (y|u) is given by The original CRN and the clamped CRN do not need to have the same initial conditions as long as the initial conditions have the same reachability classes. However, even if the two CRNs have the same initial conditions, it is possible that the clamping process will make some part of Ω Y s0|U =u0 inaccessible to C| U =u , in which case this theorem will not hold.
Theorem 2. Assume the reachability class of a detailed balanced CRN can be expressed as the product of subspaces, Ω s0 = L j=1 Ω j s0 . Then the steady-state distributions of each subspace will be independent for each product space: π(s) = L j=1 π j (s j ), where s = (s 1 , . . . , s L ) and π j is the distribution over Ω j s0 . Proof.
If Ω s0 is decomposable into a product of subspaces Ω j s0 , with j = 1...L, then each subspace involves disjoint sets of species Y j = {S j 1 , . . . , S j Mj }. In this case the steady-state distribution of a detailed balanced CRN can be factorized due to the simple nature of G(s) given by eq. (12): where s j = (s j 1 , s j 2 , . . . , s j Mj ) is state of the set of species within the subspace j.
The product form π(s) means that species from separate subspaces Ω j s0 are statistically independent. To develop non-trivial conditional probabilities for the states of different species, therefore, it is necessary either to use a non-detailed balanced CRN by driving the system out of equilibrium, or to generate complex interdependencies through conservation laws that constrain reachability classes and "entangle" the state spaces for different species. We explore both of these possibilities in the following sections.

Direct Implementation of a Chemical Boltzmann Machine (DCBM):
We first consider the most direct way to implement an N -node Boltzmann machine with a chemical system. Recall that a BM has a state space Ω BM = {0, 1} N and an energy function  . Although there are 2 2N states in which where s is any state with x OF F i = 1, and s has x ON i = 1 but is otherwise the same. Such a choice would enforce detailed balance of the CTMC, with the desired steady-state distribution To implement such a CRN, we define a reaction set R that contains a distinct pair of reactions for each possible state of the neighbors of i for which w ij = 0. Let α i ∈ {ON, OF F } |N (i)| denote a state of neighboring species. Then, the necessary reactions and rate constants are for each i and every possible state α. In physical terms, the species representing the neighbors of node i collectively catalyze X OF F i X ON i , with a separate pair of reactions for every possible α i . While this entails a large number of reactions (2 |N (i)|+1 for each node i), it allows the rate constants for each configuration of neighbors to be distinct, and thus to satisfy the ratio of rate constants given in (19). For CRN states that satisfy the conservation laws x ON i + x OF F i = 1, there will be a unique reaction that can flip any given bit, and thus the CTMC detailed balance (17) also holds, yielding the correct π(s). The construction is illustrated by example in Figure 2 and compared to other constructions in Figure 3.
The distribution π(s) is identical to that of the BM, both with and without clamping. Reachability is preserved by clamping, as all states satisfying the conservation laws and clamping can be reached in the clamped CRN. All results derived for traditional BMs therefore apply, including conditional inference and the Hebbian learning rule. The construction can be generalized to any graphical model and indeed to any finite Markov chain defined on a positive integer lattice.
With the DCBM, we have shown that CRNs can express the same distributions as BMs, and are thus very expressive. However, since each possible state α i of N (i) is associated with two reactions, the number of reactions of the CRN is exponentially large in the typical node degree d in the original BM. Moreover, the scheme requires high molecularity reactions in which multiple catalysts effect a single transition (the molecularity grows linearly with d). Physical implementations are therefore likely to be challenging. We further note that as a consequence of Theorem 2, the DCBM cannot be detailed balanced at the level of the underlying chemistry, due to its simple conservation laws. Physically, this means that the DCBM must use a fuel species to drive each reaction. Details of this argument are given in the appendix (7.1).

The Edge Species CBM Construction (ECBM):
Can a detailed balanced CRN also implement a Boltzmann machine, or is it necessary to break detailed balance at the level of the CRN reactions, as in the DCBM? Here we show that it is not necessary by introducing a detailed balanced CRN that uses species to represent both the nodes and edges of a BM.
(20) This reaction scheme is visualized in Figure 3. Just like in the DCBM, there is a separate pair of reactions for each node i for each state of its neighbors α i . In this case, however, the backwards reaction in (20) does represent a transition that is a true chemical inversion of the forwards reaction. Further, given a valid initial state, clamping any subset of the X ON/OF F i species preserves reachability. The ECBM is detailed balanced and therefore its stationary distribution has the form (13). Substituting the conservation law w ij = x i x j and omitting species with 0 energy results in Comparing this expression to the distribution of a BM, equation (2), the above expressions are equivalent provided that their partition functions are equivalent. To see this is the case, notice that: 1) the partition function is just a sum over the Gibbs factors across the entire state space.
2) The Gibbs factors take the same form between the ECBM and BM (as shown above). And 3) the reachable state spaces spaces are equivalent. Thus a sum over all possible Gibbs factors will be equal. Therefore, Z BM = Z π and the theorem is proven.
Via the ECBM, we have shown that even detailed balanced CRNs can represent rich distributions and are able to calculate conditional distributions through clamping as proven in Theorem 1. Due to being detailed balanced, this construction requires no fuel molecules and performs sampling via the intrinsic equilibrium fluctuations of the CRN. Moreover, it is only necessary to tune molecular energies in this construction, since appropriate relative rate constants follow by definition. This construction is possible due to the complex set of conservation laws that ensure that the reachability classes of all the X ON/OF F i species are tightly coupled via the W ON/OF F ij species. One implication is that this construction does not generalize easily to non-binary species counts. Additionally, issues related to high molecularity reactions and large number of reactions remain.

Approximate Bimolecular Implementations
The DCBM and the ECBM both require reactions of high molecularity. High molecularity reactions and systems involving many species are physically challenging to implement and also potentially suffer from long mixing times. In this section, we discuss an approximation scheme to create CBMs with lower molecularity reactions and thus overcome these issues.

Taylor Series Chemical Boltzmann Machine (TCBM):
Here, we demonstrate a compact CBM that approximates a BM. It is not detailed balanced on either the Markov chain or the CRN level, but uses only 2N species and unimolecular and bimolecular reactions. The TCBM is a non-equilibrium CBM of the kind discussed in Section 3.1 that uses a dual-rail representation and single-node transitions to approximately implement a BM. The reactions are given by: which, with appropriate initial conditions, preserve the conservation law that This model's parameters can be taken directly from the weights of a BM, w ij . First, define a symmetric matrix W . Decompose this matrix into the difference of two positive matrices, W = A − B, where a ij ∈ A are all w ij > 0 and b ij ∈ B are the absolute values of all w ij < 0. Finally, k is an arbitrary overall rate. This construction can be understood as an approximation of equation (17), which dictates that for two states s and s that differ only in bit i with x ON i = 1 in state s , the CTMC transition rates must satisfy The bias θ i has been absorbed into w ij for notational clarity by assuming there is some x ON 0 = 1 whose weights act as biases. The TCBM is a bimolecular CRN obeying the same conservation laws as the DCBM in which each species j acts as an independent catalyst for transitions in i with reaction rates determined by a ij and b ij . The relative propensities of this network are exactly equal to the linear expansion of the relative propensities shown in the last term in (23). Specifically, the numerator is the sum of the reaction propensities for X OF F i → X ON i and the denominator is the sum of the reaction propensities for X ON i → X OF F i , in each case plus a constant term due to the unimolecular reactions. We thus propose the simple scheme in (22) as an approximate CBM; Figure 4A depicts this TCBM schematically. This model bears some resemblance to protein phosphorylation networks where adding or removing a phosphate group is analogous to turning a species on or off; both are driven, catalytic processes capable of diverse computation.

Approximate BCRN Inference:
Remarkably, this simple approximate CBM can reasonably approximate the inferential capabilities of a BM. We demonstrate this by using (22) to convert a BM trained on the MNIST dataset [47] to a TCBM ( Figure 4). We then compare the BM and the TCBM side by side. Digit classification is shown in Figure 4 panels C and D for a BM and in Figure 4 panels F and G for a CBM as confusion heatmaps. Classification is carried out by clamping the image nodes to MNIST images and averaging the values of the classification nodes. As is apparent from these plots, the BM does a fairly reasonable job classifying these digits, but struggles on the number 5. The CBM functions as a very noisy version of the BM with nodes in general much more likely to be on. The CBM has also faithfully inherited the capabilities and limitations of the BM and similarly struggles to classify the digit 5. Digit generation is shown in Figure 4E for a BM and 4H for a CBM. Generation was carried out by clamping a single class node to 1 and all other class nodes to 0, then averaging the output of the image nodes after the network had equilibrated. For each generated image, we show the raw output and the top 85th percentile of nodes, a thresholding which helps visualize the noisy output. As is apparent from the raw output, the CBM approximation scheme does not generate images nearly as distinctly as the BM. However, this approximation does faithfully reproduce plausible digits when filtered for the top 85th percentile. Additional training and simulation details can be found in the appendix (7.2) The overall performance of the CRN is reasonable, given the fact that weights were simply imported from a BM without re-optimization. The TCBM only approximates the distribution implied by these weights and, in the absence of detailed balance, does not have an established formal relationship between clamping and conditioning.

Detailed Balanced CRN Learning Rule
A broad class of detailed balanced chemical reaction networks can be trained with a Hebbian learning rule between a waking phase (clamped) and sleeping phase (free) that is reminiscent of the classic gradient descent learning algorithm for a BM [8,35]. A similar algorithm for detailed balanced CRNs follows.
First we state a simple case of Theorem 4 where we just want a CRN with stationary distribution π over Ω s0 to learn a target distribution Q also defined on Ω s0 . Then, the learning rule is given by Here, s i π and s i Q denote the expected count of the species S i with respect to the probability distributions π and Q, respectively, and g i = G[S i ] is the energy of species S i . Theorem 4 generalizes this procedure to cases where there are hidden species.
Theorem 4. Let C = (S, R, k) be a detailed balanced chemical reaction network with stationary distribution π(s) on Ω s0 . Consider a partition (V, H) of the set S of species into visible and hidden species such that π(s) = π(v, h). Require that for all visible states v, the clamped CRN C| V =v preserves reachability. Let be the weighted mixture of stationary distributions of the clamped CRNs C| V =v with v drawn from the distribution Q. Then, the gradient of the Kullback-Leibler divergence from π V to Q with respect to the energy, g i = G[S i ], of the species S i is given by Proof. Applying Theorem 1, the clamped CRN ensemble π Q (s) may be written Additionally we will need the partial derivative of a Gibbs factor and the partition function with respect to g i , Using these results, the partial derivative of any detailed balanced CRN's distribution at a particular state s, with respect to an energy g i , is ∂π(s) Noting that Q has no dependence on g i , the gradient of the Kullback-Leibler divergence can then be written, In the special case where there are no hidden species, which is to say the target distribution Q is defined over the whole reachability class Ω s0 , then π V (v) = π(s) and π Q (s) = Q(s) and the gradient has the simple form shown in equation (24). Applying gradient descent via dgi dt = − ∂D KL ∂gi , we thus have a simple in silico training algorithm to train any detailed balanced CRN such that it minimizes the Kullback-Leibler divergence from π V to Q. If H = ∅, simulate the CRN freely to estimate the average counts s i under π(s). Then compare to the average counts under the target Q(s) and update the species' energies accordingly. If H = ∅, clamp the visible species to some v ∈ Ω V s0 with probability Q(v) and simulate the dynamics of the hidden units. Repeat to sample an ensemble of clamped CRNs C| V =q . Because clamping v preserves reachability, Gillespie simulations of the CRN with the V species clamped to the data values v will sample appropriately. This gives the average counts under π Q .
This CBM learning rule is more general than the classical Boltzmann machine learning rule, as it applies to arbitrary detailed balanced CRNs, including those with arbitrary conservation laws and arbitrarily large species counts (but still subject to the constraint that reachability under clamping must be preserved). That said, at first glance the CBM learning rule appears weaker than the classical Boltzmann machine learning rule, as it depends exclusively on mean values s i , whereas the Boltzmann machine learning rule relies primarily on secondorder correlations x i x j . In fact, though, conservation laws within the CRN can effectively transform mean values into higher-order correlations. A case in point would be to apply the CBM learning rule to the ECBM network: For x ON j π , which exactly matches the classical Boltzmann machine learning rule if we assert that the energies of OF F species are fixed at zero.

Discussion
We have given one approximate and two exact constructions that allow CRNs to function as Boltzmann machines. BMs are a "gold standard" generative model capable of performing numerous computational tasks and approximating a wide range of distributions. Our constructions demonstrate that CRNs have the same computational power as a BM. In particular, CRNs can produce the same class of distributions and can compute conditional probabilities via the clamping process. Moreover, the TCBM construction appears similar in architecture to protein phosphorylation networks. Both models are non-equilibrium (i.e., require a fuel source) and make use of molecules that have an on/off (e.g., phosphorylated/unphosphorylated) state. Additionally, there are clear similarities between our exact schemes and combinatorial regulation of genetic networks by transcription factors. In this case, both models make use of combinatoric networks of detailed-balanced interactions (e.g., binding/unbinding) to catalyze a state change in a molecule (e.g., turning a gene on/off). We note that our constructions differ from some biological counterparts in requiring binary molecular  counts. However, in some cases we believe that biology may make use of conservation laws (such as having only a single copy of a gene) to allow for chemical networks networks to perform low-cost computations. In the future, we plan to examine these cases in a biological setting as well as generalize our models to higher counts. Developing these CBMs leads us to an important distinction between equilibrium, detailed-balanced CRNs with steady state distributions determined by molecular energies, and CRNs that do not obey detailed balance in the underlying chemistry. The second category includes those that nonetheless appear detailed balanced at the Markov chain level. Physically, this distinction is especially important: a non-detailed balanced CRN will always require some kind of implicit fuel molecule (maintained by a chemostat) to run and the steady state will not be an equilibrium steady state due to the continuous driving from the fuel molecules. A detailed balanced CRN (at the chemical level) requires no fuel molecules: and thus the chemical circuit can act as a sampler without fuel cost. Despite this advantage, working with detailed balanced CRNs presents additional challenges: to ensure that chemical species do not have independent distributions, species counts must be carefully coupled via conservation laws.
Our constructions also highlight important complexity issues underlying CBM design. The number of species, the number of reactions, and the reaction molecularity needed to implement a particular BM are relevant. Trade-offs appear to arise between these different factors, and the thermodynamic requirements of a given design. A breakdown of the main features of each CBM is given in Table 1. Summarizing, the TCBM is by far the simplest construction, using O(N ) species, at most O(N 2 ) reactions, with molecularity ≤ 2. However, this happens at the expense of not being an exact recreation of a BM, and the requirement of a continuous consumption of fuel molecules. The DCBM is the next simplest in complexity terms, using O(N ) species, O(N 2 N ) reactions, and molecularity of at most N . Like the TCBM, the DCBM requires fuel molecules because it is not detailed balanced at the CRN level. The ECBM is considerably more complex than the DCBM, using quadratically more species, O(N 2 ), the same number of reactions, O(N 2 N ) and double the reaction molecularity. The ECBM makes up for this increased complexity by being detailed balanced at the CRN level, meaning that it functions in equilibrium without implicit fuel species.
Finally, we have shown that a broad class of detailed balanced CRNs can be trained using a Hebbian learning rule between a waking phase (clamped) and sleeping phase (free) reminiscent of the gradient descent algorithm for a BM. This exciting finding allows for straightforward optimization of detailed balanced CRNs' distributions.
This work provides a foundation for future investigations of probabilistic molecular computation. In particular, how more general restrictions on reachability classes can generate other "interesting" distributions in detailed balanced CRNs remains an exciting question. We further wonder if the learning rule algorithm can be generalized to certain classes of non-detailed balanced CRNs. Another immediate question is whether these ideas can be generalized to nonbinary molecular counts. From a physical standpoint, plausible implementations of the clamping process and the energetic and thermodynamic constraints require investigation. Indeed, a more realistic understanding of how a CBM might be implemented physically will help us identify when these kinds of inferential computations are being performed in real biological systems and could lead to building a synthetic CBM.

Application of Theorem 2: The Direct CBM Must Use Implicit Fuel Species
Here, we use Theorem 2 to analyze the direct implementation of a CBM and show that it cannot be detailed balanced and thereby requires implicit fuel molecules. First, notice that the the conservation laws used in this construction are of a simple form. The states accessible by (X ON i , X OF F i ) are independent of (X ON j , X OF F j ) for i = j, and therefore the reachability class is a product over the subspaces of each individual node. As a consequence, by Theorem 2, the system must be out of equilibrium and violate detailed balance at the level of the CRN because, by construction, this system is equivalent to a BM and has correlations between nodes i and j whenever w ij = 0. In physical terms, the presence of catalysts cannot influence the equilibrium yield of a species, and therefore a circuit which uses catalysis to bias distributions of species must be powered by a supply of chemical fuel molecules [41,42,43]. It is also worth noting that, as a consequence, this scheme cannot be implemented by tuning of (free) energies; it is fundamentally necessary to carefully tune all of the rate constants individually (via implicit fuel molecules) to ensure that detailed balance is maintained at the level of the Markov chain for the species of interest.

BM Training and TCBM Simulation Details:
We trained a BM using stochastic gradient descent on the MNIST dataset, down sampled to be 10 pixels by 10 pixels [47]. The BM has 100 visible image units (representing a 10 x 10 image), 10 visible class nodes, and 40 hidden nodes as depicted in Figure 4B. Our training data consisted of the concatenation of down sampled MNIST images and their classes projected onto the 10 class nodes. The weights and biases of the trained BM were converted to reaction rates for a CBM using the Taylor series approximation. This CBM consists of 300 species, 300 unimolecular reactions and 22350 bimolecular reactions. The resulting CBM was then compared side-by-side with the trained BM on image classification and generation. The BM was simulated using custom Gibbs sampling written in Python. The CRN was simulated on a custom Stochastic Simulation Algorithm (SSA) [40] algorithm written in Cython. All simulations, including network training, were run locally on a notebook or on a single high performance Amazon Cloud server.
Classification was carried out on all 10000 MNIST validation images using both the BM and the CBM. Each 10 by 10 gray-scale image was converted to a binary sample image by comparing the gray-scale image's pixels (which are represented as real numbers between 0 and 1) to a uniform distribution over the same range. The network's image units were then clamped to the binary sample and the hidden units and class units were allowed to reach steady state. This process was carried out 3 times for each MNIST validation image, resulting in 30000 sample images being classified. Raw classification scores were computed by averaging the class nodes' outputs for 20000 simulation steps after 20000 steps of burn-in (Gibbs sampling for the BM, SSA for the CBM). Max classification was computed by taking the most probable class from the raw classification output. Raw classification and max classification confusion heatmaps, showing the average classification across all sample images as a function of the true label are shown in Figure 4 panels C and D for a BM and in Figure 4 panels F and G for a CBM.
Image generation was carried out by clamping the class nodes with a single class, 0...9, taking the value of 1 and all other classes being 0, and then allowing the network to reach steady state. Generated images were computed by averaging the image nodes over 50000 simulation steps (Gibbs sampling for the BM, SSA for the CBM) after 25000 steps of burn-in. Generation results are shown in Figure 4E for a BM and Figure 4H for a CBM.