**William Poole** (1), Tom Ouldridge (2), Manoj Gopalkrishnan (3), and Erik Winfree (1)
(1) Caltech
(2) Imperial College London
(3) IIT Bombay

*To Navigate. Down Arrow: More Details. Right/Left Arrows: Prev/Next Section. Escape: Poster View.*

An version of this notebook can also be found on google colab

Stochastic Chemical Reaction Networks (CRNs) represent a biophysically realistic and expressive modeling framework for biological systems. Detailed Balanced CRNs (dbCRNs) are a subclass of CRNs which admit an equilibrium distribution. This poster builds off our previous work on Chemical Boltzmann Machines (Poole et al. DNA 2017) where we have shown that dbCRNs can exactly implement a stochastic neural network model called a Boltzmann Machine. This relationship illustrates that stochasticity can be advantageous because it allows a CRN to compactly represent complex distributions. We have shown that these CRNs are capable of inference, defined as computing conditional distributions, by “clamping” a subset of their species to a particular value. Finally, we have shown that the equilibrium distributions of dbCRNs can be learned by an external in silico non-CRN algorithm that tunes the energies of individual molecular species.

In new work, we have shown that the above results regarding inference and learning apply to any dbCRN, not just the Chemical Boltzmann Machine constructions. In this interactive poster, we examine physical processes that can tune species’ energies and are potentially implementable in a cell-like environment. The first uses a special kind of detailed balanced chemostat to control carefully designed chemical potentials that maintain detailed balance. The second uses a simple feedback circuit implemented as a non-detailed balanced CRN. These constructions provide equilibrium and non-equilibrium physical implementations that control the inference and learning process in dbCRNs. We then analyze these processes thermodynamically to provide lower bounds on the costs of inference andlearning. These results illustrate possible mechanisms whereby a biochemical system in a smallvolume, such as a cell, can represent and adapt to its environment.

Boltzmann Machines (BMs) are an inspiration for this work. Briefly, they are a class of energy based machine learning models. BMs link machine learning and physics and can be defined equivalently in terms of graphical models, generalized Ising models, or stochastic Hopfield networks. In terms of this work, BMs provide two important sources of inspiration:

P(V, H) | P(V | H) | P(H | V) |

Chemical Reaction Networks (CRNs) are a set of species $S$ and reactions $R$ of the form $\sum_i I_i S_i \to \sum_i O_i S_i$ where $I_i, O_i \in \mathbb{Z}_{\geq 0}$ are integer number of inputs and outputs of species $S_i$.

In this work, we consider stochastic chemical reaction networks which model species' counts on the integer lattice as a Markov jump process as opposed to modeling the concentrations of molecular species. This is appropriate when considering small volumes with relatively few molecules, where stochasticity is known to be important.

We further are primarily interested in a subclass of CRNs called detailed balanced (dbCRNs). There are many important theorems related to dbCRNs. In particular, they are amenable to thermodynamic and statistical mechanical analysis because they admit an equilibrium distribution.

Note: we do not restrict ourselves to bimolecular CRNs because (with some $\epsilon$ of error) high order reactions can be represented as multiple bimolecular reactions.

- The
**Reachability Class**of a CRN are the set of states $\Omega(x_0)$ reachable via any series of reactions from an initial condition $x_0$.

- The
**Stoichiometric Compatability Class**of a CRN is the kernel of the stoichiometric matrix $\mathbb{S} = O - I$, where $O^r_i$ and $I^r_i$ are the outputs and inputs of species $i$ in reaction $r$, respectively.

- A CRN is detailed
**Detailed Balanced**if:- that every species $S_i$ has an associated energy $E_i$.
- all reactions are reversible $\sum_i I^r_i S_i \to \sum_i O^r_i S_i \implies \sum_i O^{r'}_i S_i \to \sum_i I^{r'}_i S_i$.
- reaction rates rates are set by the detailed balanced condition $\frac{k_r^+}{k_r^-} = e^{\Delta E_r}$ where $\Delta E_r = \sum_i I^r_i E_i - O^r_i E_i$ where $r$ is the index of a reaction with forward and reverse rates $k_r^+$ and $k_r^-$.

dbCRNs can produce 3 classes of distribution: (Left) **product poisson distributions** cover the entire integer lattice. (Middle) **unimodal correlated distributions** are created restricting the stoichiometric compatability class of the CRN. (Right) **complex multimodal distributions** are created by carefully controlling the reachability class. These CRNs are examined in more detail below.

This product poisson distribution is generated by the CRN:

It has the cannonical poisson form:

Notice that this distribution is in product form meaning all species are independent at equilibrium.

*Changing the energies of each species will only change the mean of that species. Changing the initial condition will not change the distribution because the entire integer lattice is reachable.*

This restricted product poisson distribution is generated from the CRN:

Although the equilibrium distribution can also be written in cannonical "poisson" form, the partition function is no longer factorizable but instead must be written as a function of the reachability class $\Omega(x(0))$ which in turn depend on the initial condition. In this case, the reachability class is the stoichiometric subspace.

Due to the particular reaction network chosen, species will be anti-correlated with eachother.

*Changing the energies of each species will change mean of that species and also induce changes in the means and variances of the other species. Changing the initial condition will change which subset of the integer lattice is reachable also changing the distribution. However, the general shape of this distribution: 3 anti-correlated variables, is robust to changes of the energies and initial condition.*

The previous example hinted that one species in a dbCRN can strongly influence another species' distribution. This example uses that premise to design a CRN with many hidden species organized into a heirarchical network in order to generate a very specific distribution. This example is a modified from Capalletti et al. 2019 where dbCRNs are constructed to produce arbitrary distributions with finite support. The CRN is divided into three sets of reactions

Happy Reactions: $C_h + L_h^{(i, j)} + i X_1 + j X_2 \leftrightharpoons C_h + L_h^{(l, m)} + l X_1 + m X_2$

Sad Reactions: $C_s + L_s^{(i, j)} + i X_1 + j X_2 \leftrightharpoons C_s + L_s^{(l, m)} + l X_1 + m X_2$

Transition Reactions: $C_h + L_h^{(i, j)} + i X_1 + j X_2 \leftrightharpoons C_s + L_s^{(l, m)} + l X_1 + m X_2$

Here, $C_h$ and $C_s$ are binary "control" species which determine if a happy or sad face is formed. $L_s^{(i, j)}$ and $L_h^{(i, j)}$ denote happy and sad "location" species which represent the $i, j$ pixels seen in each icon. The energies of the locations are counteract the entropies of the $X_1$ and $X_2$ counts in order to exactly produce the distribution.

Conservation laws: $C_h + C_s = 1 \quad \sum_{i, j} L_s^{i, j}+L_h^{i, j} = 1$

The energies of $C_h$ and $C_s$ can be tuned to bias the distribution towards happy or sad. Changing $E_1$ and $E_2$ will further skew the distributions. Finally, this construction is reliant on there only being a single location species and a single control species present at any time. If these conservation laws are broken, the reachability class rapidly expands and breaks the control over the equilibrium distribution.

Note that although this CRN uses high order reactions, it could be accureately implemented using only biomolecular reactions for example by iteratively forming molecular complexes which then transition.

We have demonstrated that changing species' energy pushes that species towards a higher or lower value. This form of energy modulation is also a form of inference which we call **clamping**. For example, by decreasing the energy of a species we are making it more favorable and increasing its mean - effectively "holding" that species at a higher value. We have proven that this is equivalent to computing a conditional distribution conditioned upon the "rare event" that a different mean was sampled from the equilibrium distribution.

What this means, practically, is that varying the energies of a dbCRN's species is a form of inference and can be used to compute conditional distributions!

CRN: $\emptyset \leftrightharpoons X_1 \quad \quad \emptyset \leftrightharpoons X_2 \quad \quad \emptyset \leftrightharpoons X_3$

CRN: $X_1 \leftrightharpoons X_2 \leftrightharpoons X_3$

CRN: See previous slide.

What does it mean to change the energy of a species? Evolutionarily, changes in sequences (nucleic, peptide etc.) could modulate binding energies, providing a connection between mathematical inference and evolution.

There are also more immediate ways to control the species' energies. One implementation comes from coupling each species $X_i$ to an external chemical potential $P_i$ such that whenever $X_i$ is changed in any reaction, $P_i$ is changed with it.

In chemistry and physics, the term $E^P_i + \log [P_i]$ is commonly known as the chemical potential $\mu$. This illustrates that inference via clamping via changing the energies of a detailed balanced CRN can be accomplished by modulating a chemical potential.