Probabilistic Inference with Polymerizing Biochemical Circuits

Probabilistic inference—the process of estimating the values of unobserved variables in probabilistic models—has been used to describe various cognitive phenomena related to learning and memory. While the study of biological realizations of inference has focused on animal nervous systems, single-celled organisms also show complex and potentially “predictive” behaviors in changing environments. Yet, it is unclear how the biochemical machinery found in cells might perform inference. Here, we show how inference in a simple Markov model can be approximately realized, in real-time, using polymerizing biochemical circuits. Our approach relies on assembling linear polymers that record the history of environmental changes, where the polymerization process produces molecular complexes that reflect posterior probabilities. We discuss the implications of realizing inference using biochemistry, and the potential of polymerization as a form of biological information-processing.

Probabilistic inference—a procedure for estimating unobserved variables in probabilistic models—has been used to describe various aspects of cognition [1,2]. In this line of work, organisms are thought to build probabilistic models of their world and use these to guide action and perception via inference. Although the work on biological instantiations of probabilistic inference has focused on animal nervous systems [3,4], single-celled organisms also show complex, history-dependent behaviors in changing environments. This history-dependent character has been described in cognitive terms, such as “learning,” “memory,” and “decision-making” [5,6,7,8,9,10,11,12,13,14,15].

The internal state of microbes, for instance, is shaped by past experiences [16,17,18]. B. subtilis populations exposed to distinct environmental perturbations and then grown in the same environment can be distinguished (based on gene expression) for as long as 24 h after the perturbation, suggesting that cells retain a “memory” of past environments [16]. E. coli cells grown in an environment that switches periodically from glucose to lactose eventually adapt to the switches and maintain a consistent growth rate through the switches [17]. Yeast strains anticipate their environment by producing a metabolic program needed to metabolize an alternative nutrient, even before the preferred nutrient is depleted [19]. All these behaviors occur on an ontogenetic timescale, much faster than that of mutation and natural selection.

These observations raise the question of whether probabilistic inference can be realized in cellular biochemistry. In principle, Chemical Reaction Networks (CRNs), commonly used to formalize biochemical systems, can approximate any computable function [20,21]—and therefore implement inference. Several studies have shown how CRNs can implement probabilistic inference [22,23,24,25,26]. For instance, [22] showed it is possible to construct a CRN that, at steady-state, encodes the joint probability distributions of probabilistic models known as factor graphs, and suggested enzyme-free DNA strand displacement [27] as a physical realization. However, strand displacement is an unlikely cellular mechanism. Moreover, cells act as the environment changes [13,24], and while the abundances of cellular components change in a stochastic manner [28,29], which raises an important question: which biologically plausible, cellular mechanisms can support inference under such stochastic conditions? Some studies have explored phosphorylation and transcriptional regulation as vehicles for forming memories and simple associations. For example, [10] designed biochemical circuits, which use either phosphorylation or transcriptional control, that can be conditioned through a unicellular analog of Hebb’s rule, while [30] evolved chemical circuits in silico for associative learning. Such mechanisms could potentially be used to build inference-performing circuits.

Yet biochemistry also includes generative and combinatorial mechanisms, such as polymerization, which are not part of the typical repertoire of molecular mechanisms (e.g., protein phosphorylation or transcriptional regulation) that are thought to “compute” or “process information”. Through polymerization, the cell produces structures such as microtubules and actin cables, which are highly regulated [31] and crucial to a variety of cellular functions. While these polymers are often studied for their mechanical properties [32], it has also been recognized since the early days of molecular biology that polymerization can be viewed as a computational process, capable of implementing logical automata [33] (more recent studies have also explored the computational power and properties of polymers; see [34] and references therein).

In this paper, we follow a similar line of thought, by exploring how the process of polymer assembly can be used to realize probabilistic inference. We show how the process of linear polymer assembly can be used to perform inference in a Markov model. Our circuit uses the polarity of polymers and their constituents to create a molecular record of the environment’s history (its past sequence of changes), which we show can be used to anticipate the environment’s dynamics and regulate other biochemical programs. This suggests that polymerization—due to its ability to produce macromolecules out of many combinations of parts—can be a useful motif of biological computation. The paper is organized as follows. We first introduce the use of probabilistic inference to anticipate an environment that changes according to a Markov model, and derive the representations and operations needed to approximately implement inference in real-time in such an environment. We then describe a circuit that can perform these operations by assembling linear polymers. We describe several properties of the circuit, and show how the probabilistic information it records can be used to regulate a specific response to environmental change. We close with a discussion of the properties of this circuit and potential future directions.

Probabilistic models have been used to represent a variety of dynamic, uncertain environments. One of the simplest probabilistic models is the discrete-time, finite-state Markov model, in which the state of the environment at time t is assumed to depend only on the prior states going back to the t−k time point (where k is the order of the model; when k=1, we have a first-order Markov model). While this Markov assumption is violated by many natural processes [35], it will serve as a useful idealization for understanding how biochemical circuits that perform inference can be constructed.

We consider a changing environment that can be in one of two states, A or B, and where switches between states are driven by the Markov model. Such models are parameterized by two transition probabilities: the probability of switching from A to B, πAB, and the probability of switching from B to A, πBA—as shown in Figure 1A. The transition probabilities are unobserved. Through probabilistic inference, one could estimate the values of these probabilities, and use this information to anticipate the environment’s switches—potentially in a way that would be useful to an organism adapting to a changing environment.

To guide our search for biological mechanisms that could instantiate inference in this model, we analyze the task of inference in this model through David Marr’s three levels scheme [36]. Marr’s scheme is a general “top-down” approach to understanding biological computation that begins by describing the computational level (level 1): what are the “inputs” and “outputs” to the task, and constitutes a “correct” solution? Given this computational formulation of the task, the algorithmic level (level 2) asks for a suitable representation and procedure for carrying out the task. Finally, the hardware level (level 3) asks for biologically plausible mechanisms that can realize the representation and algorithm from level 2. The idea is that each level guides and constrains the next level. We will apply this approach to derive the computational representations and operations required for inference in our probabilistic model (levels 1 and 2), and then construct a biological circuit guided by these requirements (level 3).

We first constructed a Markov chemical environment whose dynamics are driven by the Markov model shown in Figure 1A. We can think of this chemical environment as a bioreactor in which organisms grow, and where A and B are analogous to nutrients which are flowed into the reactor by an experimenter and consumed by circuit. The Markov chemical environment is created as follows. We used sequences sampled from the discrete-time Markov model to add/remove A and B at fixed time intervals, as follows. Initially, a state X0∈{A,B} is sampled from the Markov model and a fixed amount of X0 is added to a reactor which contains only our circuit (no A or B). After a fixed time interval, we sample Xt+1 from the model given the previous state Xt. If Xt+1=Xt, no perturbation is performed; otherwise, we remove all A and B present in the reactor and add Xt+1. An environment generated by this procedure using a Markov model where πAB=πBA=0.95 is shown in Figure 1B.

Assuming the environment’s perturbations are generated by a Markov model, probabilistic inference can be used to anticipate the environment’s states. In particular, it is useful to compute the probability of encountering A or B next given past observations of the environment’s states, history=Xt,Xt−1,Xt−2,…, where Xt corresponds to the state of the environment at time t. In Bayesian terms, anticipation of the environment means computing the posterior predictive distribution, P(Xt+1∣history). This computation is complicated by the fact that the transition probabilities πAB and πBA are unknown. One option in such cases is to place prior probabilities on these parameters and then integrate them out to calculate, ∫P(Xt+1∣history,πAB,πBA)dπABdπBA. We take an alternative and sometimes simpler strategy of estimating the unobserved transition probabilities, πAB and πBA, through Bayesian inference, and use these estimates to anticipate the next state: (1)P(πAB,πBA∣history)∝P(πAB,πBA∣history)P(πAB,πBA)

If a mathematically convenient prior distribution is chosen for P(πAB,πBA), then Equation (1) can be solved exactly (see Appendix A). Using the posterior distribution, we can then obtain estimates of the transition probabilities, π^AB,π^BA, and use these to anticipate the next state by sampling Xt+1∼P(Xt+1∣Xt,π^AB,π^BA):P(Xt+1=A∣Xt,π^AB,π^BA)=1−π^ABifXt=Aπ^BAifXt=B

An important complication in our setting is that cells act while the environment is changing, so the posterior distribution (Equation (1)) must be estimated in real-time. However, this computation simplifies considerably when we consider the representation and algorithm needed to estimate the posterior distribution (Marr level 2). Crucially, the sequence of observations about the environment, history, does not need to be stored in full; it can be compressed into a matrix T of transition counts (Figure 1C). Assuming n possible states of the environment, T is an n × n matrix whose ith row indicates the number of times the environment switched from the ith state to each of the states. Since the sum of the ith row must equal the total number of transitions from state i that have been observed, the matrix can be summarized by n2 − n entries. In an environment with two states, the sufficient statistics are simply two counts: the number of times the environment switched from A to B, and from B to A (along with the total number of transitions; see Appendix A for details).

Using these counts as the representation, the following elegant algorithm for computing the posterior distribution in real-time emerges. As the environment changes, update the relevant counts in T. Then anticipate the next state given the current state i by: (1) normalizing the ith row of T to convert counts into probabilities, and (2) sampling a state from that row. This algorithm relies on a representation of counts (as stored in T), and two computational operations. First, it is necessary to update the counts by a transition update operation (Figure 1D) that distinguishes A to B transitions from B to A transitions. Second, in order to use these counts in downstream computations, we need to be able to convert the counts into probabilities by normalization and sample from the resulting distribution; this is the normalization-sampling operation shown in Figure 1E. With these three ingredients—the transition counts from T, an operation that updates the transition counts, and an operation that normalizes and samples from a row of T—the posterior distribution can be estimated in real-time.

To realize inference biochemically based on the above analysis, we need a biochemical circuit that can perform the operations shown in Figure 1C–E. Such a circuit would need a representation of the transition matrix T, which encodes the number of switches between relevant states of the environment. Two operations on T would then need to be realized molecularly: (1) the ability to count directionally, i.e., to record when the environment has switched from state A to B (as opposed to from B to A) and store this information in T, and (2) the ability to access a row of counts in T, normalize it, and sample from the resulting probability distribution. This is challenging since directional counting requires recording a potentially unbounded number of switches in the environment. Some digital counters proposed in the synthetic biology literature, such as [37], have a fixed capacity—determined by the number of genetically components—and are thus not suitable.

We instead constructed a biological realization of the counts matrix and its operations using polymerization. We make use of the idea that the process of assembling a macromolecule such as a polymer can be seen as a computational process. The rules that say how to add a new part (monomer), Xn+1, to an existing polymer X1-X2⋯Xn, depending on the state of the new part Xn+1 and the configuration of the existing polymer’s end, Xn, describe a logical flow [33]. A simple example would be the rule: if Xn+1≠Xn, where Xn+1 is a monomer, then add Xn+1 after Xn. If two types of monomers are available, e.g., Xi=A or Xi=B, then this process would generate alternating polymers such as A-B-A or B-A-B-A. The polymerization process, then, is akin to running a kind of logical automata, and the resulting polymer’s sequence is a record of the computation, similar to the output tape of a Turing machine. In our case, however, the assembly process is not sequential and deterministic (as it is in a Turing machine), but rather driven by concurrent, stochastic biochemical reactions.

We use polymerization to represent directional changes in the environment by building a set of linear polymers, which we call transition history polymers, made up of A and B molecules. Each polymer represents the sequence of transitions that have occurred in the environment. The process of assembling the polymer will allow us to count the number of environmental transitions.

The transition history polymers are assembled as follows. We assume that the A and B molecules have a distinguishable “head” site and “tail” site, and that each polymer is built starting from a T monomer (which also has a head and tail site). Polymerization then proceeds according to two sets of rules: (i) an unbound T can bind the head of a free A or B molecule (i.e., those molecules whose head and tail sites are unbound), forming a TA or TB dimer (Figure 2A, top); and (ii) the unbound tail of a TA or TB dimer can bind the head of a free A or B (Figure 2A, bottom). This second set of rules records the direction of change by producing memory molecules corresponding to each of the four possible transitions: A-To-A memory molecule upon A to A binding, A-To-B upon A to B binding, B-To-A upon B to A binding, and B-To-B upon B to B binding. For simplicity, we assume all these polymerization rules are irreversible and proceed according to the polymerization rate constant k. The process of polymer assembly thus implicitly “senses” the environment, and uses the polarity of the polymer’s constituents to record the direction of change. (Note that if A and B are viewed as nutrients, e.g., two different kinds of sugars consumed by microbial cells a bioreactor, we cannot simply assume that these molecules could polymerize into arbitrary linear polymers. However, we could always posit a third “scaffold” protein, which binds A and B, and can differentially polymerize based on its binding state, using the rules in Figure 2A. For simplicity, we assume A and B have those polymerizing capacities directly).

As a consequence of polymerization, the abundance of A-To-B and B-To-A will encode the number of transitions from A to B and B to A, respectively, while the abundance of A-To-A and B-To-B will correspond to the duration of each state (more on this point below). We will refer to the A-To-B and B-To-A memory molecules as transition-encoding memory molecules, because their abundance tracks switches, and the A-To-A and B-To-B memory molecules as duration-encoding memory molecules, because their abundance tracks the environmental state’s duration (a point we will return to below). We further assume that the memory molecules are degraded at a constant rate. Thus, the polymer assembly process will produce memory molecules that, if relatively stable, reflect the environment’s history of transitions, as we will see below. Note that the number of polymers that the system will construct will depend on the abundance of T monomers (which serve as the starting monomers for polymers). The number of polymers being assembled will in turn affect the abundance of memory molecules A-To-B and B-To-A, a point that we will return to below.

The additional operation needed for realizing inference is normalization-sampling (Figure 2B). Normalization-sampling means normalizing the transition counts, and then sampling from the resulting probability distribution. On a traditional digital computer, this is straightforward to implement sequentially. A standard algorithm for sampling from a probability vector θ (in effect, sampling from a multinomial distribution) is: calculate a cumulative sum up to each element θi, then iterate through θ, draw a uniform random number w∈[0,1] for each entry θi, and return the first θi such that ∑j=1iθj≥w. Within a cell, however, there is no readily accessible notion of sequential iteration to support this procedure; a concurrent mechanism is needed instead.

We can make use of the fact that concentrations lend themselves to encoding probabilities (also discussed in [22]) to derive a non-sequential version of the same computation. Consider two proteins, an activator and a repressor, that have expression levels c1 and c2, respectively, (Figure 2B). We can normalize and sample from the vector θ=[c1c2] by having the activator and repressor bind, in mutually exclusive fashion, a third molecule which we call an integrator. If the integrator is not in excess of the activator and repressor, then the fraction of activator-bound integrator will be proportional to the concentration of activator normalized by the sum of concentrations of activator and repressor (Figure 2B). This means that a molecular interaction that depends on the concentration of activator-bound integrator will occur in proportion to the probability c1c1+c2. This mechanism does not depend on sequential order: the interactions between activator, repressor, and integrator can all take place concurrently (similar integrator schemes have also used by other circuits that sense environmental changes in a way that is sensitive to order of exposure, e.g., [38]).

To encode the posterior probability of encountering A (or B) given the environment’s history, we can use the activator-repressor-integrator scheme to “normalize” the transition counts encoded by memory molecules (Figure 2C). When the environment is in state A, the A-To-A and A-To-B transition molecules will bind, in proportion to their abundance, to the A-To-B-integrator. The fraction of integrator molecules bound by A-To-B (relative to those bound by A-To-A) will represent P(A∣history) (likewise for A-To-B and P(B∣history). One complication is that the duration-encoding memory molecules, A-To-A and B-To-B, whose abundance is determined by the environmental state’s duration, will inevitably be more abundant than A-To-B and B-To-A. Thus, the rate constants governing the interactions between memory molecules and integrators have to be adjusted to account for this bias (see Appendix A).

We have implemented this circuit using Kappa, a rule-based language for describing and simulating stochastic biochemical systems (such simulations correspond formally to continuous-time Markov chains as described in [39,40]). Rule-based representations are especially suited for modeling polymerizing circuits such as this one, which may correspond to an intractably large number of ODEs (due to the combinatorial explosion of species [41]). We next explore several features of this circuit when exposed to changing, probabilistic environments.

In order to simulate our model, we had to choose the number of polymers T that the circuit has to work with. The resulting circuit is sensitive to the environment’s history. When simulated (using Kappa) in a periodic environment (πAB=πBA=0.95), the circuit accumulates transition-encoding memory molecules at each switch (Figure 3A,B), as expected.

We next compared the circuit’s “posterior”—as encoded in the fraction of active integrators—to the posteriors based on the idealized, discrete-time Markov model. The posterior estimates of πAB,πBA in the idealized model are obtained by looking only at the sequence of switches, A,B,A,B,…, with a uniform prior over the two states (see Appendix A).

The fraction of active integrators qualitatively matches the posterior estimates of the idealized model (Figure 3C,D). As the circuit accumulates A-To-B memory molecules, for example, the fraction of active A-To-B integrators—i.e., those integrators bound by the A-To-B memory molecule—increases, following the analytic estimates of the transition probabilities (Figure 3C,D). Moreover, the linear polymers at the end of the simulation reflect the history of periodic switches (Figure 3E).

When exposed to an environment with different statistics, the circuit behaves accordingly. We considered an environment where A frequently switches to B, but switches from B to A are rare, meaning B is a “sticky” state (πAB=0.85,πBA=0.25) (Figure 4A). In this environment, the transition-encoding memory molecules track the switches, but the duration-encoding memory molecule B-To-B reflects that the environment tends to dwell in state B (Figure 4B). As a result, the fractions of active A-To-B and B-To-A integrators qualitatively track the analytic posteriors (Figure 4C,D). For instance, after the onset of B (time 50), the fraction of active B-To-A decreases as the environment stays in B (Figure 4C, time 50–750). The linear polymers assembled in this environment reflect this environment’s history (Figure 4E). Thus, the circuit is not limited to one particular Markov environment, but rather uses the process of polymer assembly to produce memory molecules that reflect the environment’s transition history.

We next compared our circuit’s behavior to the idealized probabilistic model more systematically, by exposing the circuit to many different assignments of the transition probabilities πAB,πBA (Figure 5). For each environment, we first estimated the median of the posterior over πAB at each time point and then took π^AB to be the median of those estimates (and similarly for πBA to obtain π^BA). We then computed a similar quantity for our circuit by taking the median of the fraction of bound A-To-B integrators through time (I^AB) and likewise for B-To-A integrators to obtain I^BA.

Since these quantities are not directly, quantitatively comparable, we compared the ratio of posteriors from the idealized model (π^ABπ^BA) with the ratio of the relevant active integrators (I^ABI^BA). The comparison shows broad agreement (rank correlation ρ=0.94) across different assignments of transition probabilities (each assignment was simulated 100 times, Figure 5). It is important to note that this can only be a qualitative comparison, since the biological realization of the posterior and the analytic model begin from different starting points (zero active in the circuit’s case, versus a 50%-50% prior in the analytic model). Another caveat is these ratios being compared only reflect end-state results, not differences in the temporal trajectory of posterior ratios. We also stress again that the analytic model which operates in discrete-time is, by definition, an idealization that could never perfectly match the biological circuit which operates in continuous time and obviously lacks synchronous switches from one time point to another. In the next section, we describe some properties of our polymerizing circuit and its difference from the idealized discrete-time model.

While our circuit was guided by an analysis of inference in an idealized Markov model, these results indicate an important difference between that abstract model and our circuit as realized using polymerizing biochemical parts. Unlike in the abstract model we had started with, the circuit consumes the very “signal” that it is responding to, as As and Bs get incorporated into the circuit’s linear polymers. Thus, the circuit shapes the environment it inhabits, rather than passively observing its dynamics.