Stochasticity, bistability and the wisdom of crowds: a model for associative learning in genetic regulatory networks

It is generally believed that associative memory in the brain depends on multistable synaptic dynamics, which enable the synapses to maintain their value for extended periods of time. However, multistable dynamics are not restricted to synapses. In particular, the dynamics of some genetic regulatory networks are multistable, raising the possibility that even single cells, in the absence of a nervous system, are capable of learning associations. Here we study a standard genetic regulatory network model with bistable elements and stochastic dynamics. We demonstrate that such a genetic regulatory network model is capable of learning multiple, general, overlapping associations. The capacity of the network, defined as the number of associations that can be simultaneously stored and retrieved, is proportional to the square root of the number of bistable elements in the genetic regulatory network. Moreover, we compute the capacity of a clonal population of cells, such as in a colony of bacteria or a tissue, to store associations. We show that even if the cells do not interact, the capacity of the population to store associations substantially exceeds that of a single cell and is proportional to the number of bistable elements. Thus, we show that even single cells are endowed with the computational power to learn associations, a power that is substantially enhanced when these cells form a population.

Almost all animals can associate neutral stimuli and stimuli of ecological significance [1]. An extensively studied example is eye-blink conditioning (Figure 1) [2], [3]. Naïve rabbits respond to an airpuff to the cornea (Unconditioned Stimulus, US) with eyelid closure (Unconditioned Response, UR). By contrast, a weak auditory or visual stimulus (Conditioned Stimulus, CS) does not elicit such an overt response. Repeated pairing of the CS and the US forms a cognitive association between the CS and the US such that the trained animal responds to the CS with eyelid closure, a response known as Conditioned Response (CR). Two important characteristics of associative learning are (1) specificity and (2) generality. The CR does not reflect a general arousal. Rather, the animal learns to respond specifically to the CS. The generality is reflected by the fact that a large family of potential stimuli can serve as a CS if paired with the US.

Neuronal networks are particularly adapted to performing this association and in the last few decades there has been considerable progress in understanding the ways in which experience-based changes in synapses in the nervous system underlie this associative learning process [4], [5]. Neural network models for associative memory, which explain how both specificity and generality are maintained, are typically based on three elements: (1) Synapses are the physical loci of the memory; (2) synaptic plasticity underlies memory encoding; (3) neural network dynamics, in which the activities of neurons depend on the synaptic efficacies, underlie the retrieval of the learned memories in response to the CS.

Genetic regulatory networks (GRN) describe the interaction of genes in the cell through their RNA and protein products [6], [7], [8]. Previous studies have pointed out the similarity between the dynamics of GRNs and the dynamics of neural networks [9]. For example, GRNs, like neural networks, can implement logic-like circuits, where the concentration of a protein (high or low) corresponds to the binary state of the gate [10], [11], [12]. These findings prompted us to evaluate the capacity of GRNs to learn associations.

Considering associative learning in animals, the US is typically a stimulus of biological significance, such as food or a noxious stimulus that elicits a response (UR) in the naïve animal, either in the form of muscle activation or gland secretion. The GRN correlate of a pain-inducing stimulus is stress. Stressful conditions such as heat, extreme pH, or toxic chemicals often result in a substantial change in the expression level of many different proteins in the cell. For example, Escherichia coli (E. coli) bacteria respond to a variety of stress conditions by a general stress response mechanism in which the master regulator controls the expression of many genes [13]. These stressful conditions can be regarded as a US and the resultant change in the expression level of the proteins can be regarded as a UR. By contrast, other stimuli may result in a narrow or absence of a response of the cell and in that sense can be referred to as potential CS. Learning in this framework would correspond to the formation of an association between these potential CS and US such that following the repeated pairing of the CS and US, the presentation of the CS would elicit a UR-like response (CR).

The responsiveness of the GRNs to different stimuli has been shown to change over time in response to evolutionary pressure in a manner that resembles associative learning [14], [15]. These changes take place on time scales that are substantially longer than the lifetime of a single cell and in contrast to associative learning in animals, entail modifications of the genome through mutations. On a shorter timescale, there is some evidence that the single-celled Paramecium can learn to associate a CS with a US within its lifetime [16]. However, these findings have been disputed [17] and the question of whether Paramecia can learn associations and the characteristics of this learning await further experimental validation. The capacity of GRNs to learn associations in shorter, non-evolutionary time-scales has also been studied theoretically using GRN models. Learning in these models is restricted to a small subset of predefined stimuli [18], [19], [20], [21] and thus the computational capabilities of these GRN models are limited compared to neural network models.

Here we show that a GRN based on bistable elements and stochastic transitions can learn associations while retaining both specificity and generality. We further compute the capacity of the network and show that the number of different learned associations that the network can simultaneously retain is proportional to the square root of the number of bistable elements. Moreover, this capacity is substantially enhanced when considering a clonal population of GRNs. These results imply that even bacteria are endowed with the capacity to learn multiple associations.

Our Genetic Associative Memory model (GAM) for associative learning is based on three components: (1) a memory module that provides the long time-scale necessary for the maintenance of memories for long periods of time; (2) a mechanism for encoding the desired memories and (3) a response mechanism for the readout or retrieval of the stored memories in response to the relevant stimuli. We describe the three components separately in the simpler case of a predefined association and then generalize to the case of multiple associations and to the case of a population of GAMs.

A necessary condition for associative learning is the ability to maintain memories. Memories require a long time-scale, which characterizes multistable dynamics [22], [23]. For example, the ability of flip-flop devices in electronic circuits to maintain memories is based on their bistable dynamics. Bistability naturally emerges in dynamical systems if two conditions are fulfilled: positive feedback and saturation [24], [25], [26]. Both these requirements characterize many GRNs, and bistability has been found in both artificially engineered [27], [28], [29], [30] and natural GRNs [31], [32], [33], [34], [35], [36]. For example, the response of the lactose promoter in E. coli to intermediate induction levels was shown to be bistable because of the positive feedback loop on the import of the inducer in the cell [37].

In our GAM, we assume a positive feedback loop between a gene and its protein product. The gene encodes for a protein M which binds cooperatively, as a transcription factor, to the promoter of that gene, resulting in further synthesis of M. The kinetic reactions describing the dynamics of M appear in Text S1 in the Supporting Information, and their deterministic approximation [7], [8] is equivalent to those of positive feedback loops such as the lac system [38], [39]:(1)where reflects the nonlinear positive feedback (see Eq. (4) in the Materials and Methods). The second term in Eq. (1) denotes the protein degradation, where is a parameter. The third term models the contribution of external factors to the dynamics of M (see below).

The function , depicted in Figure 2A (top), is N-shaped and is characterized by three zero-crossings. The two outermost zero-crossings (red arrows in Figure 2A, top) correspond to the two stable states (or fixed-points): a low expression level of M, (left) and a high expression level, (right). It can be readily shown that the intermediate zero-crossing (black arrow) corresponds to an unstable fixed-point of the dynamics. Thus, the dynamics of Eq. (1) converge to one of the equilibrium values, the low or high expression level of M, depending on the initial conditions. This bistability of the dynamics of M endows the GRN with the capacity to store binary memories in the form of the level of expression of M. For this reason we refer to M as a ‘pseudo-synapse’.

It is useful to rewrite Eq. (1) using an ‘energy’ function such that where . The energy function (Figure 2A, bottom) is characterized by two minima (red arrows) and one local maximum (black arrow). The two minima correspond to the two stable fixed points of the dynamics of M and the maximum corresponds to the intermediate unstable fixed point.

The differences between the value of the energy function at the maximum and the values at the minima are known as the energy gaps and are denoted by . In Figure 2A the two energy gaps are approximately equal. However, an increase in the value of raises the function (Figure 2B top), resulting in a smaller energy gap for the fixed point, and a larger energy gap for fixed point (Figure 2B bottom). By contrast, a decrease in the value of lowers the function (Figure 2C top), resulting in the opposite effect: a larger energy gap for , and a smaller energy gap for (Figure 2C bottom).

It should be noted that Eq. (1) is a deterministic approximation of the biochemical dynamics. Biochemical processes such as the bursting activity of the transcriptional machinery are stochastic [40], [41], [42], [43]. One way to account for the stochasticity is by adding white noise to Eq. (1) such that(2)where is a Gaussian white noise, , , is the magnitude of the noise and is the Dirac delta “function”.

One consequence of this stochasticity is that the noise is expected to occasionally induce transitions between the two fixed points. A well-known result from the field of stochastic processes is that if the noise is sufficiently weak, the rate of transitions from one minimum to the other is exponentially dependent on the energy gap, , (e.g., see [44], and see also [45] for a more accurate approximation). Consequently, even a small change in the energy gap is expected to result in a large change in the transition rate. Thus, although the three energy functions in Figure 1A–C are qualitatively similar, they represent very different dynamics. For sufficiently weak noise, the rate of transition from to in Figure 2B is negligible compared to the rate of transition from to . Similarly, the rate of transition from to in Figure 2C is negligible compared to the rate of transition from to . Moreover, the rates of transition between the two stable states in Figure 2A are both negligible compared to the rate of transition from to in Figure 2B or the rate of transition from to in Figure 2C. Thus, the transitions between the states are highly dependent on the value of . We utilize these results when modeling the memory encoding in the next section.

In associative learning, memory is encoded in response to the contiguity of the CS and the US. To implement this idea in the framework of the GAM, we assume that the value of is determined by external cues, the CS and US. Formally, we assume that the CS and the US induce the expression of proteins C and U, respectively. The value of is determined by the concentrations of C and U such that the US is effectively a repressor of M but the co-occurrence of US and CS activates the expression of M (see Eq. (4) in the Materials and Methods). In other words, U in isolation decreases but when bound to C it increases . This mode of regulation has already been observed; e.g., in the osmotic response regulatory system in yeast [46].

For simplicity we assume in our model that the expression levels of C and U are binary, and , reflecting the presence or absence of the CS and US, respectively. Moreover, we assume that independently of the external cues, the value of is such that the dynamics of the pseudo-synapse are bistable (as in Figure 2). Thus, the co-occurrence of the CS and US increases the transition rate to the high expression level of M (as in Figure 2B) whereas an exposure of the GAM to the US alone increases the transition rate to the low expression level of M (as in Figure 2C).

The computational implications of these dynamics are that a repeated exposure of the GAM to the co-occurrence of the CS and US is expected to result in a high state of M, whereas a repeated exposure of the GAM to the US in the absence of the CS is expected to result in a low state of M. In this sense, the state of M is the physical correlate of the memory of the association between the CS and US and a high level of M indicates an association between the CS and US. Assuming that in the absence of the US, the two energy gaps are high (as in Figure 2A), the transition rates between the two states of the pseudo-synapse, in both directions, would be low. Thus, in the absence of the US, information about the existence of an association, as well as its absence, would be maintained for long periods of time. These dynamics are reminiscent of a multiplexer. A multiplexer is a device that selects one of several input signals and forwards the selected input into a single line. In the dynamics of M, the US selects whether M will be maintained (in the absence of the US) or whether the value of M is determined by the CS (in the presence of the US), as in [18]. However, in contrast to a standard multiplexer, transitions in our model are stochastic. Thus, the dynamics depicted in Figure 2 resemble a stochastic multiplexer. This difference implies that multiple repetitions are needed in order to change the state of M with a high probability.

The last component of our GAM is a readout scheme that decodes the state of M in the presence of a CS such that the CS elicits a response if and only if the expression level of M is high. To implement this, we assume that the UR and the CR manifest in the GAM as the production of a response protein R. The expression of R is regulated by two mechanisms: the US regulates the expression of R through the binding of U to a promoter of R, and the CS-pseudo-synapse pair regulates the expression of R by cooperative binding of C and M to another promoter (Figure 3A). The kinetic reactions describing the dynamics of the expression of R appear in Text S1 in the Supporting Information, and their deterministic approximation is given by:(3)where is the degradation rate of R and the functions and describe the dependence of the expression of R on U, C and M (see Eq. (5) in the Materials and Methods). A high level of U in Eq. (3) results in a high value of . This elicits the expression of R, independently of the values of M and C, corresponding to the UR. By contrast, a high level of C results in a high level of only when M is in its high expression level. Thus, in the absence of the US, the stimulus substantially increases R only when M is in its high expression level, corresponding to the CR.

The dynamics of Eqs. (1–3) describe a GAM that can learn the association between a CS and a US. This is demonstrated in Figure 3B. Initially, at time t = 0, M is in the low expression level state, corresponding to the ‘naïve’ state of the network prior to learning. In this state, a US (orange rectangle, t = 1 h) elicits a response (UR), but a CS (open blue rectangle, t = 2 h) does not elicit any response. Following two pairings of the CS and US (t = 3 and 4 h), the state of M does not change but in response to the third pairing (t = 5 h) the state of M changes to the high state level. In this state, the GAM is responsive to both a CS (t = 6 and 7 h) and a US (t = 8 h). Three presentations of the US in the absence of a CS (t = 8–10 h) do not elicit any change in the state of M but in response to another presentation of the US in the absence of a CS (t = 11 h), the state of M reverts to its low value, and the GAM is no longer responsive to CS (at t = 12 and 13 h).

In the previous section, we demonstrated that a GRN can learn to associate a CS with a US (Figure 3). However, this learning is limited, as it is specific to a single, predefined CS. This GAM can be trivially generalized to enable the learning of several different associations by postulating that the GAM is characterized by a number of memory elements, each associated with a single CS. However, this generalized GAM is still limited in its ability to learn associations because only those predefined CS can be learned. This limitation contrasts sharply with neural network models, which are capable of learning general associations. In this section we generalize the model presented in Figure 3A and show that similar to neural network models, GRNs are also endowed with the capacity to learn a large number of arbitrary, overlapping CS.

Consider the network described Figure 4A. In contrast to the single-pathway model (Figure 3A), in which a CS induces the expression of a single protein C, in the generalized model we assume that the CS are complex stimuli that activate N different receptors, , . Each receptor is associated with a single pseudo-synapse . The dynamics of each of the pseudo-synapses follow the same equations as in the single-pathway model (not shown in Figure 4A, see Eq. (4) in the Materials and Methods).

The last component of our generalized GAM is the readout scheme. We assume that similar to the single-pathway model, the UR and the CR manifest in the generalized model as the production of a response protein R. We assume two independent promoters that regulate the expression of R. The response to the presentation of the US is described by (Eq. (3)) and the response to the presentation of the CS is regulated by the cooperative binding of and , where different pairs of and independently regulate R (Eq. (5) in the Materials and Methods and Text S1 in Supplementary Information).

For simplicity, we assume in our analysis that the patterns of expression of the proteins that define the stimuli are random and independent. In this case, the statistics of the stimuli are fully determined by the sparseness of the stimuli, the probability that is in its high expression level, .

To gain insights into the ability of the generalized GAM to learn multiple associations, we consider a naïve GAM, in which the values of the pseudo-synapses are random (Figure 4B, bottom, t = 0). The responses of the GAM to five different stimuli, denoted by A, B, C, D and E, presented to the GAM at times t = 0, 1, 2, 3 and 4 h, respectively, are relatively small. This is due to the random, and hence relatively small overlap between the pattern of activation of pseudo-synapses (color coded) and the pattern of activation of the receptors of the five stimuli ( is denoted by an open blue rectangle in Figure 4B).

In response to the pairing of C, B and A with the US (at times t = 5, 6 and 7 h, respectively), the expression levels of some of the become more similar to that of the in A, B and C, respectively. As a result, the GAM responds more vigorously to the presentation of A, B and C (at times t = 8, 9 and 10 h, respectively) but not to the presentation of D or E (at times t = 11 and 12 h, respectively). However, as a result of a repeated association of pattern E with the US (at times t = 13, 14 and 15 h), the GAM vigorously responds to the presentation of pattern E (at time t = 17 h) but not to pattern D (at time t = 16 h). This example demonstrates that a GAM can selectively learn to associate several arbitrary CS patterns with a US.

A careful analysis of Figure 4B reveals that after learning, the magnitude of the responses to the three learned CS is not equal. The response to stimulus C (t = 10 h) is smaller than the response to stimulus B (t = 9 h) and the response to B is smaller than the response to A (t = 8 h). This difference reflects the fact that the order of association affects the magnitude of the response to a CS. This is because learning a new pattern may change the expression level of a pseudo-synapse that participates in the encoding of an older pattern. For example, consider pseudo-synapse 4 in Figure 4B. In response to the presentation of stimulus C (at time t = 5 h), the state of the pseudo-synapse has changed to the high expression level, in line with the expression level of in CS C. However, the association of the US with A (at time t = 7 h) has reverted the state of the pseudo-synapse to the low expression level, decreasing the overall response to the CS C. In other words, the association with the CS A has overwritten the information stored in pseudo-synapse 4 concerning the CS C. More generally, because of the overwriting of memories by more recent memories, the magnitude of response to a CS is expected to decrease with the number of subsequent CSs. After the encoding of a large number of patterns, the response to an ‘old’ CS is expected to diminish to an extent where it is no longer distinguishable from the response to non-learned stimuli. In this case the CS is said to have been extinguished (a more precise definition of “distinguishable” appears below). By contrast to the diminishing of the response to a pattern following the overwriting by other patterns, the repeated co-occurrence of the same pattern with the US (at times t = 13, 14 and 15 h) augments the strength of association of that pattern with the US, as demonstrated by the response to pattern E at time t = 17 h.

The magnitude of the order effect depends on two probabilities: the probability p that the co-occurrence of U and a high level of expression of would induce a transition from to in the corresponding pseudo-synapse and the probability q that the co-occurrence of U and a low level of expression of the corresponding would induce a transition from to in the corresponding pseudo-synapse . The probabilities p and q are determined by the two rates of the US-induced transitions and the duration of co-occurrence of the US and CS, T (assuming that the rates of all other transitions are negligible, see above) such that and , where and are the low-to-high and high-to-low transition rates, respectively. The larger the transition rates and the longer the duration, the larger the transition probabilities are.

If all pseudo-synapses are determined by the most recent CS and the pattern of expression level of the different pseudo-synapses corresponds to the pattern of activation of the receptors in that CS. As a result, the response to the most recent CS is substantially larger than the response to a non-learned stimulus. However, this comes at a price. The most recent CS overwrites the memory trace of all previously encoded CS and therefore the responses to all these ‘older’ CS are indistinguishable from the responses to the non-learned stimuli. Thus, if , the GAM cannot store more than a single association. The smaller the values of p and q (e.g., due to smaller US-induced transition rates), the fewer pseudo-synapses change in the process of learning a CS, allowing the GAM to maintain information about previously-learned CSs.

However, the transition probabilities should not be too small because the smaller these probabilities are, the weaker is the encoding. If these probabilities are too small, the response of the GAM even to the most recently stored GAM is too small to be distinguishable from non-learned stimuli. Therefore, in order for the GAM to be able store a large number of CS, the values of the US-induced transition rates should be sufficiently large to allow for a sufficiently large response to the learned-CS but sufficiently small to minimize the overwriting of old memories by new memories.