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

To better understand the requirement that the response to a CS needs to be distinguishable from the response to non-learned stimuli, consider again Figure 4B. The responses of the GAM to the presentations of the non-learned stimuli A-E at times 0–5 h, respectively, are not identical. These differences are due to the fact that there is stochasticity in the response, resulting from stochasticity in the dynamics of the pseudo-synapses and in the realization of the different CS. Therefore, a memory of a CS is said to be maintained if the distribution of the responses of the GAM to the CS is distinguishable from the distribution of responses to the non-learned stimuli. This notion becomes exact in the next section.

How many CS can be stored in a GAM? Addressing this question using the full dynamical equations (Figure 4) requires extensive simulations that are beyond the scope of this paper. Therefore, we use a binary approximation (see Materials and Methods). The quality of the binary approximation is demonstrated in Figure S1 in the Supporting Information.

As described in the previous section, responses to non-learned stimuli depend on the overlap of the pattern of activation of the stimuli with the pattern of activation of the pseudo-synapses. Because both the stimuli and the dynamics of the pseudo-synapses are stochastic, this response is a stochastic variable. The distribution of the responses to non-learned stimuli (see Eq. (14) in the Materials and Methods) is depicted in Figure 5A (blue). The response of the GAM to learned CS is also a stochastic variable. The distribution of responses to the most recently learned CS is depicted in Figure 5A (black; Eq. (13) and in the Materials and Methods). This distribution is well-separated from the distribution of responses to the non-learned stimuli. Therefore, recently-learned stimuli are distinguished from non-learned stimuli using a simple threshold mechanism (e.g., the dashed line in Figure 4B). The probability of an error depends on the overlap between the two distributions. If the overlap is small, the GAM almost always responds to the most recently learned CS and almost never responds to non-learned stimuli. On the other hand, a large overlap would result in a large number of errors, false positives or misses, depending on the choice of threshold. The difference between the means of the two distributions (black and blue) depends on the transition probabilities. The higher the probabilities, the larger the difference is. Therefore, the higher the transition rates are, the easier it is to distinguish between the most recently learned CS and the non-learned stimuli.

The distribution of responses to the presentation of the second-most recently learned CS (darkest gray) is also to the right of the distribution of responses to non-learned stimuli (blue). Nevertheless, it is shifted to the left relative to the distribution of responses to the most recently learned CS (black). As a result, the overlap of this distribution with the distribution of responses to the non-learned stimuli is larger. The reason for this shift is that as noted in Figure 4B, the newer CS ‘overwrites’ the memory of the older CS, resulting in a decreased overlap between the CS and the pseudo-synapses. The degree of overwriting, manifested as a shift to the left of the distribution of responses to the second-most recently learned CS relative to the most recently learned CS, depends on the US-induced transition rates. The smaller the transition rates, the smaller the overwriting is and therefore the smaller the shift to the left of the distribution.

More generally, the distributions of responses to a CS shift to the left with the ‘age’ of the CS. This is depicted in Figure 5A using grayscale. While the distribution of the several most-recently learned CS is well-separated from the distribution of responses to non-learned stimuli (blue in Figure 5A), the distributions of responses to ‘older’ CS and non-learned stimuli largely overlap, indicating that ‘older’ CS are ‘forgotten’.

More formally, the ability of the GAM to distinguish between a learned CS and a non-learned stimulus depends on the signal to noise ratio (SNR), which is defined as the ratio of the difference in mean responses to the two classes of stimuli, divided by the square root of the sum of variances of the two distributions (Eq. (14) in the Materials and Methods). In general, the larger the SNR, the fewer errors when distinguishing between learned and non-learned stimuli. The SNR, as a function of the ‘age’ of the CS is depicted in Figure 5B: the newer the CS, the larger the SNR. The SNR of the CS (where the numbering of patterns is reversed such that corresponds to the most recent stimulus) is given by Eq. (14) in the Materials and Methods section. The capacity of the GAM can thus be defined as the ‘oldest’ CS such that the corresponding SNR is larger than 1. In other words, the capacity of the GAM is defined as the largest value of such that .

The capacity of the GAM depends on the US-induced transition rates, which determine the transition probabilities. As discussed above, if these rates are high, forgetting is fast. On the other hand, if these rates are too low the GAM cannot reliably retrieve even the most recent CS. The capacity of the GAM is maximal when the US-induced transition rates are intermediate, balancing between these two requirements. The capacity of the GAM as a function of the number of pseudo-synapses (N) is depicted in Figure 5C (blue). The larger N, the larger is the capacity of the GAM. In the Materials and Methods section we show that in the limit of , if the US-induced transition probabilities are optimal, the capacity of the GAM is proportional to the square root of the number of pseudo-synapses, (Eq. (20); red line in Figure 5C). This result is similar to the memory capacity of models of neural networks with binary synapses [47], [48]. However, the learning rule proposed here, even in the binary approximation, differs from the Hebbian synaptic plasticity rule used in neural network models [47], [48].

In the previous sections we studied the ability of a single GRN to learn associations. However in nature, GRNs often do not reside in isolation but in populations comprising of a large number of individual cells of the same type, e.g., as in a colony of bacteria or in a tissue, all exposed to the same external conditions. This raises an interesting question: is the capacity of a population of GAMs to store associations larger than that of a single GAM? The answer is trivially positive if we allow the different GAMs to communicate and form a recurrent network with specialized connections between individual GAMs, similar to neurons in neuronal networks. However, here we ask a different question: is the capacity of a population of non-interacting GAMs to store and retrieve memories different from that of the single GAM?

We consider a population of generalized GAMs as in Figure 4A. All GAMs are identical, exposed to the same sequence of stimuli but differ in their internal stochasticity. In other words, the noise associated with the dynamics of the pseudo-synapses (Eq. (2)) in the different GAMs is assumed to be independent. The population response in our model is assumed to be simply the accumulated response of all individual GAMs.

In order to understand why the capacity of a population of identical GAMs to store memories may be larger than the capacity of a single GAM, we note that a CS of a particular ‘age’ can be retrieved if the overlap between the distributions of responses to the learned and non-learned stimuli is sufficiently small. This overlap is sensitive to the variances of the two distributions (width of the curves in Figure 5A). The larger the variance, the larger is the overlap. Two sources contribute to this variance in the responses. First, there is stochasticity in the realization of CS and non-learned stimuli. Second, there is stochasticity in the encoding process. While the first type of stochasticity is external and thus shared by all GAMs in the population, the second type of stochasticity is independent for each GAM. As a result, when considering the cumulative response of a large population of GAMs, all other parameters being equal, the variance in the distribution of responses is considerably smaller (Eq. (23) in the Materials and Methods). In Figure 6A we plot the distributions of responses to CS of different ‘ages’ (gray, color-coded) and non-learned stimuli (blue).

Similar to the case of a single GAM, the capacity of a population of GAMs depends on the US-induced transition rates. However, because the variance in the responses in the case of the population is considerably smaller than the variance in the case of a single GAM, the US-induced transition probabilities that maximize the capacity of the population are considerably smaller than those that maximize the capacity of a single GAM. In Figure 6B we plot the SNR as a function of the ‘age’ of the CS (solid blue line). Compared to the SNR of a single GAM (dashed blue line, identical to Figure 5B), the SNR of the response of the population of GAMs is larger than 1 for much ‘older’ CS.

The capacity of the population of GAMs as a function of the number of pseudo-synapses (N) is depicted in Figure 6C (solid blue line). The larger N is, the larger is the capacity of the GAMs. More quantitative analysis reveals that for an appropriate choice of parameters, the number of different CS that a large population of GAMs can store is proportional to the number of pseudo-synapses (Eq. (31); solid red line), compared to a capacity that is only proportional to the square root of the number of pseudo-synapses in the case of a single GAM (dashed blue line, identical to Figure 5C).

In this paper, we explored the ability of a general GRN to encode associations. We showed that a GRN that is endowed with bistable elements and stochastic dynamics is capable of storing and retrieving multiple arbitrary and overlapping associations. The capacity of a single GRN in our model, defined as the number of stored associations, is proportional to the square root of the number of bistable elements . This result is reminiscent of Hopfield-like models with bounded synapses, in which the capacity is proportional to the square root of the number of synapses [47], [48], [49]. Remarkably, in a large population of GRNs, as is in a colony of bacteria or in a tissue, this capacity is substantially higher and is proportional to the number of bistable elements.

Despite the similarities between the GAM and the Hopfield model, there are two important differences that are noteworthy. First, the capacity of a single GAM may be limited by the presence of readout noise (e.g., in the dynamics of R). However, this readout noise is not expected to substantially affect the capacity of a population of GAMs because of averaging. Second, the number of neurons available in neuronal networks is much larger than the number of bistable elements in GRNs. Altogether, our model predicts that if the number of bistable elements in the GRN does not exceed several tens, it will be difficult to store more than one or two memories in a single GAM. Therefore, the storage of multiple memories is likely to require a population of GAMs.

The key elements in our model are the bistability and the stochasticity of the dynamics of the GRN. Importantly, bistability and stochasticity are not restricted to the transcriptional machinery. Rather, they are found in various cellular processes, including post-transcriptional regulation (e.g., by non-coding RNA [50], [51]) or post-translational regulation (e.g. phosphorylation and degradation regulation [36], [52], [53], [54]). We modeled associative memory that is based on the interaction of proteins through the transcriptional machinery because these dynamics are better characterized and are more accessible experimentally than other cellular alternatives.

Moreover, the GAM is not restricted to a particular organism. The parameters used in the simulations presented in this paper are biologically plausible for bacteria. However, because the basic elements of the GAM, namely, bistability and stochasticity, are widespread in GRNs of all cells, the potential for associative learning without a nervous system exists for virtually all cell types, including single-celled eukaryotes and plants. Furthermore, this work suggests that even in animals that possess a nervous system, learning that is independent of this nervous system is also possible. In particular, it could be interesting to consider GAM in the immune system, which has evolved to learn to respond to novel pathogens.

Bearing this in mind, we believe that in view of the recent developments of experimental methods that quantitatively measure the expression level of proteins, bacteria, in particular the well characterized E. coli, are the ideal substrate to study the associative learning in GRNs. Each of the components of the GAM module (Figure 3A), namely inducible elements, bistable switches and AND gates, have been established in the E. coli transcription network and therefore a synthetic implementation is achievable [55], [56].

Beyond synthetic implementation, the complexity of the genetic networks suggests that GAM-like modules may exist. A first step in searching for GAMs in known networks should be the identification of plausible candidates for the US, UR and CS. In animals, the US is a stimulus that causes an overt response prior to learning, the UR. Typically the US is a stimulus of biological significance, such as food or a noxious stimulus and the UR is an ecologically-relevant overt response, often in the form of muscle activation. For example, in the eye-blink conditioning experiment (Figure 1) the US is an air puff and the UR is an eye blink that protects the eye from the puff. An important point to consider when searching for associative learning in bacteria is ecological significance. Our model for associative learning, similar to models of associative learning in neuronal networks, does not incorporate any ecological information about the stimuli. However in animals, it is known that the ability to form an association depends on the ecological relevance of the CS to the US. For example, the association of the taste of a certain food (CS) with the symptoms caused by a toxic or spoiled food (US), known as taste aversion, is easily-formed after a small number of repetitions. By contrast, it is substantially more difficult to form an association of a tone with the same US [57]. It is generally believed that this difference results from the fact that typically, taste is more informative about the chemical composition of substances than auditory signals. Therefore, taste-aversion but not tone-aversion has evolved as a specific learning mechanism aimed at preventing the consumption of poisonous substances. Drawing an analogy to associative learning in bacteria, we propose to utilize ecologically-relevant CS rather than arbitrary CS when searching for associative learning in bacteria.

In our model, the strength of association increases with the number of repetitions due to the stochasticity in the encoding process. Such dependence of the strength of association on the number of repetitions is also observed in classical conditioning experiments in animals [58]. Therefore, experiments involving a large number of co-occurrences of the CS and US are more likely to reveal associative learning in GRNs or populations of GRNs. Note that standard experiments studying responses of bacteria are typically short and do not involve repetitions in the presentation of stimuli to the same population of bacteria. Therefore, associative learning in such experiments may have been overlooked. Moreover, we have shown that the learning capacity of the population of bacteria is higher than that of a single GRN. Therefore, the experimental search for associative learning in bacteria should be done at the population level.

More specifically in bacteria, the presence of foreign bacteria is a signal of potential stress. For example, many bacteria produce antibiotics that are harmful to other strains [59]. Other bacteria are sensitive to these damaging antibiotics and respond to their presence by activating a pre-wired stress response, such as the multiple antibiotics response (MAR) [60]. We thus suggest that the R gene in our scheme corresponds to one of the outputs of MAR response, e.g. the micF gene [61]. Note that similar to the blink in the classic eye-blink conditioning that protects the eye from the air puff (Figure 1), the activation of micF prevents the entry of the antibiotics into the cell. Thus, the antibiotics can be considered as a US whereas the stress response can be considered as a UR. However, the production of harmful antibiotics is not present in all bacteria species. Therefore, learning to distinguish between harmful and benign strains of bacteria is of potential great ecological significance because it may allow the bacteria to respond faster. Thus, the presence of foreign bacteria could correspond to the CS in our framework. Indeed, bacteria are able to detect secondary metabolites that are produced by other strains [62]. In that line, we suggest as a candidate for the M protein in the model the MarA gene. MarA is known to positively autoregulate itself, and thus has the potential to be bistable. In addition, the promoter of that gene contains multiple binding sites for transcription factors, allowing for complex regulation of the gene expression including the realization of AND gates.

Experimentally, the UR can be measured using a fluorescent-based reporter that is regulated by a promoter of a stress response gene. The CS in this framework should be stimuli that can be sensed by the bacteria but do not elicit the stress response. These include a change in the concentration of different molecules that does not activate the stress response. Repeated exposure to such conditions can be controlled using a chemostat [63], which can maintain selected growth conditions at a constant level while changing others.

Finally, the benefit of the stress response at the population level can also be found in the induction of the MAR response, as it triggers the activation of genes that inactivate toxic compounds. The benefit of this “pooled response” for the population comes from the decrease in the concentration of the toxic compound [64].

Whether or not associative learning exists in GRNs on a time-scale much shorter than required for evolution is an open question. However, whether considering bacteria that can predict a stress condition or human digestive cells that can predict food intake, associative learning in single and populations of cells seems to have an evolutionary advantage. In view of the computational capabilities of GRNs demonstrated in this paper, we believe that future careful investigations will reveal the existence of associative learning in single and populations of cells.

In this section we provide the rate equations that describe the dynamics of the generalized model (Figure 4A). The kinetic reactions that underlie these dynamics and the derivation of the rate equations from the kinetic reactions are described in Text S1 in the Supporting Information. The single pathway equations correspond to the generalized model with .

The equations that describe the dynamics of the pseudo-synapses are given by:(4) where . The nonlinear positive feedback term, , is described by where the parameters are kinetic parameters, and n is the Hill coefficient, corresponding to the cooperativity of binding. The second term in Eq. (4) denotes the protein degradation, where is a parameter. The third term in Eq. (4) describe the effect of and U, where are kinetic parameters. The last term in Eq. (4) models the stochasticity of the dynamics, and we assume that are independent white noise such that , where is Kronecker's delta function such that if and otherwise and is a parameter.

The reactive equation that describes the dynamics of R is given by:(5)where is the degradation rate of R, and where and are kinetic parameters.

In this section we compute the capacity of the GAM to learn associations. To that goal, we consider a binary approximation of the dynamics of the pseudo-synapses. Because the dynamics of M spend most of the time near the attractors of the deterministic dynamics, Eq. (1), it can be approximated using a two-state Markov chain, where each state corresponds to one attractor of the deterministic dynamics. We further assume that the US and CS are presented in discrete “trials” composed of a fixed period of time. Therefore, the response of M to the presentation of the CS and US can be approximated by:(6)where and such that and denote epochs in which and , respectively and and denote epochs in which and , respectively. The variables m and m′ denote the states of the pseudo-synapse before and after the presentation of the external cues and their values; 0 or 1 denote epochs in which and , respectively. The steady state response to the presentation of a pattern is:(7)The selectivity of the response in Eq. (7) depends on the value of the sum of . In response to the presentation of a CS that was learned CSs ago,(8)where(9) and (see Text S1 in Supplementary Information for a more detailed derivation).

Dissociating a learned pattern from non-learned patterns (which we denote as ) is possible only if is significantly different from . The difficulty in dissociating learned and non-learned patterns lies in the fact that the responses to the two types of patterns are stochastic variables that depend on the stochasticity in the realization of the learned and non-learned stimuli as well as the stochasticity in the learning. Therefore, we consider the distribution of responses to the learned and non-learned stimuli.

To compute the distribution of , note that in response to the presentation of a sequence of CS, changes in the state of the pseudo-synapses follow a Markov chain such that(10)From Eq. (10) it follows that at the stationary distribution,(11)where .

Using Eq. (11), and the fact that and , a straightforward calculation yields that the mean and variance of are given by:(12)where(13)Note that for large N is the sum of a large number of independent and identically distributed random variables and therefore according to the central limit theorem is normally distributed.