Probabilistic adaptation in changing microbial environments

which means \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$P({S}_{1}=i\mid {\pi }_{{s}_{1}})={\pi }_{{s}_{1}}^{(i)}$\end{document}PS1=i∣πs1=πs1i, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${\pi }_{{s}_{1}}^{(i)}$\end{document}πs1i denotes the ith entry in πs1. For clarity, we will sometimes omit the explicit value assignment for a random variable, and write P(St = i∣πs1) as simply P(S1∣πs1). The switch states at later time points are generated as follows: each switch state i ∈ {1, …, K} has an associated probability vector πi, whose jth entry \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${\pi }_{i}^{(j)}$\end{document}πij determines the probability of transitioning from switch state i to j. The probability of the switch state at time step t > 1 therefore depends on the previous switch state’s value St−1 and the switch state transition probabilities:

Similarly, each nutrient state can take one value i ∈ {1, …, J}, and the initial nutrient state C0 is drawn from a probability distribution, πc0, which is in turn drawn from a Dirichlet prior distribution:

The probability of a nutrient output at time t > 0 depends on the previous nutrient Ct−1 and on the switch state at time t. The switch state value indexes which transition probability distribution to use for the outputs, and as before, the transition probabilities are drawn from a Dirichlet prior:

The posterior predictive distribution, P(Ct+k∣C0:t), is the main quantity of interest (we assume k = 1 as in all of our simulations). This distribution can be calculated recursively in dynamic probabilistic models, a property that we will use in a later section to derive a real-time estimation procedure for this distribution using particle filtering. (For an accessible introduction to recursive estimation of Bayesian posteriors and to particle filtering, see Ch. 1 in Murphy (2002), Ch. 1-3 in Särkkä (2013) or Cappé, Godsill & Moulines (2007)). We derive the posterior predictive distribution in steps. First, consider the posterior distribution over the switch state at time t given the environment history up until and including time t, P(St∣C0:t), called the filtering posterior distribution, which is obtained by marginalizing out the switch state St−1:

Note that the third term is the filtering posterior over the switch state at time t − 1, which can be rewritten as we did above in terms of the filtering posterior at t − 2, and so on—this shows that the posterior can be computed recursively. The base case of the recursion is given by our prior distributions over the initial nutrient state C0 and the initial switch state S1 (as shown in Fig. S2).

The posterior predictive distribution P(Ct+1∣C0:t) can then be written as a product that uses the filtering posterior, marginalizing out the hidden switch states:

The distributions that the filtering posterior decomposes to depend on parameters that are unobserved, such as the transition probabilities. Since Dirichlet-Multinomial distributions are conjugate (Gelman et al., 1995), we can analytically integrate out these parameters. For example, the transition probabilities π for the switch states can be integrated out:

Similarly, the probability of observing a nutrient given the previous nutrient and the switch state, P(Ct+1∣Ct, St+1), can be computed while integrating out the nutrient transition probabilities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${\pi }_{i,j}^{{^{\prime}}}$\end{document}πi,j′:

We discuss in detail how to estimate these distributions from observations in the section on particle filtering.

Values of hyperparameters:In all analyses, we set hyperparameters αs1, αc0, αc to be vectors consisting of all ones. The hyperparameter on switch state transitions αs was set such that self-transitions get the hyperparameter value 2, i.e., for the ith switch state, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${\alpha }_{s}^{(i)}=2$\end{document}αsi=2, and all other entries in αs are set to 1. This encodes a weakly “sticky” prior that slightly favors self-transitions for hidden switch states.

We estimate the posterior predictive distribution above in real-time using particle filtering (Cappé, Godsill & Moulines, 2007) as outlined in Algorithm S1. Particle filtering for our model is implemented in the particlefever library (available on Github.) We used 200 particles for all simulations with particle states initialized from the prior.

In particle filtering, a distribution of hidden state values (in our model, the hidden switch states) is represented using a set of particles. Each particle corresponds to a configuration of the hidden states (configurations are typically assigned from the prior distribution at time t = 0). The particles are assigned weights that are initialized to be uniform. Starting with time t, the particle filtering algorithm works by first predicting hidden state values particles t + 1. When a data point at time t + 1 is observed, the weight of each particle is updated to be proportional to the likelihood of the new data point given the particle’s hidden state configuration. This process is repeated as data points are observed. To prevent particle degeneracy (a case where particles get weights that are too low), a resampling step is used where particles are sampled in proportion to their weights and the weights are reset to be uniform.

We now turn to the representation of our state space that is encoded in each particle. As discussed above, because of the conjugacy of Dirichlet-Multinomial distributions (Gelman et al., 1995), we can analytically integrate out the transition probabilities in our model. This means that these transition probabilities don’t have to be represented in our particles. Instead, the sufficient statistics’ for our model are simply: (1) a matrix of counts S where S(i,j) is the number of times hidden switch state i transitioned to hidden state j in the particle’s trajectory, and (2) a three-dimensional array of counts C where C(i,j,k) is the number of times nutrient i transitioned to nutrient j under switch state k. The predict/update cycles for a particle p = {st, S, C}, where s is the particle’s hidden state at time t, are:

Prediction step:For each particle p we draw a new switch state for t + 1, st+1 ∼ P(St+1∣St = st). By conjugacy, the switch state transition probabilities πst can be integrated out, yielding the posterior predictive distribution for a Dirichlet-Multinomial:

Updating step:When a nutrient ct+1 is observed, we update the weight w of our particle p = {st+1, S, C} in proportion to P(Ct+1∣Ct, St+1): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}\begin{eqnarray*}w\propto P({C}_{t+1}={c}_{t+1}\mid {C}_{t}={c}_{t},{S}_{t+1}={s}_{t+1}). \end{eqnarray*}\end{document}w∝PCt+1=ct+1∣Ct=ct,St+1=st+1.Integrating out the nutrient transition probabilities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${\pi }_{{c}_{t},{c}_{t+1}}^{{^{\prime}}}$\end{document}πct,ct+1′, this also gives the Dirichlet-Multinomial posterior predictive distribution: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}\begin{eqnarray*}w\propto \frac{{\alpha }_{c}+{\mathbf{C}}^{({c}_{t},{c}_{t+1},{s}_{t+1})}}{\sum _{k}({\alpha }_{c}+{\mathbf{C}}^{({c}_{t},k,{s}_{t+1})})} . \end{eqnarray*}\end{document}w∝αc+Cct,ct+1,st+1 ∑kαc+Cct,k,st+1.See Algorithm S1 for remaining details.

We simulated growth with different policies using a simple model of exponential growth. Cells were assumed to grow exponentially with a growth rate determined by the environment’s nutrient state. The initial population size and the time duration of each environment simulated are as described in figure legends. To determine growth rates empirically, we fitted splines to the log of the population sizes across the time and computed the first derivative. All code for fitness simulations is available in the paper’s Github repository.

The nutrient transition counter model was drawn in CellDesigner (version 4.4) (Funahashi et al., 2008), serialized as an SBML file, and simulated in Python using libRoadRunner. The SBML file for the model is available at 10.6084/m9.figshare.3493994. A detailed report of the chemical reactions and rate parameters was generated by SBML2 (Drager et al., 2009) and is available at 10.6084/m9.figshare.3492185.

We first asked whether an adaptive strategy, which exploits the probabilistic structure of the environment, would pay off in terms of growth compared with a non-adaptive strategy. We considered a class of Markov nutrient environments that fluctuate between two nutrients, glucose and galactose, where nutrient changes follow a Markov model. The probabilistic structure of these environments is determined by two parameters: the probability of transitioning from glucose to glucose (θGlu→Glu) and the probability of transitioning from galactose to glucose (θGal→Glu), as shown in Fig. 1A. (Equivalently, the model can be parameterized by θGlu→Gal and θGal→Glu, since θGlu→Gal = 1 − θGlu→Glu.) Intuitively, the higher θGlu→Glu and the lower θGal→Glu, the more likely we are to encounter glucose in the environment. Different settings of the transition probabilities can produce qualitatively different environments (Fig. 2A). We also assume that one nutrient (in this case, glucose) confers a higher growth rate than alternative nutrients, which is typically the case.

In our framework, the behavior of a cell population is determined by a policy: a mapping from the environment’s past state to a future action. We compared the performance of two policies in Markov environments: a posterior predictive policy, in which cells tune to the most probable nutrient (based on θGlu→Glu and θGal→Glu), and a non-adaptive strategy in which cells are constitutively tuned to the preferred nutrient, glucose. The quantity of interest in the posterior predictive policy is the posterior predictive distribution, which is the probability of a nutrient at the next time step given the previously observed nutrients: P(Ct+1∣C0:t), where Ct+1 denotes the nutrient at time t + 1 and C0:t denotes the environment’s nutrient history up until t. The posterior predictive policy chooses the nutrient that maximizes this probability (see ‘Materials & Methods’ for details).

To compare the fitness difference between these policies, we simulated population growth with each policy using a highly simplified growth model, similar to one used in Jablonka et al. (1995). In this model, we assume that cells respond to changes in the environment with a delay, or “lag.” Cells tune to a nutrient at time t and incur a change in growth rate as a consequence of this decision at a later time t + k, where k is the lag parameter (we assume here that k = 1). For growth kinetics, we assume that: (1) cells grow exponentially when their nutrient state matches the environment’s state, (2) there is no switching cost for cells between nutrient states, and (3) when cells are “mismatched” to their environment—i.e., tuned to a nutrient that is not present—their growth rate is zero (this assumption is supported by the observation that cells lacking Gal4, a transcription factor required to activate the GAL pathway, cannot grow in galactose alone (Escalante-Chong et al., 2015)). Altogether, our growth assumptions represent extreme conditions, but they serve as a useful starting point for seeing when an adaptive strategy is beneficial to population growth.

Given these growth assumptions, we plotted for each Markov environment the ratio of expected growth rate using the posterior predictive policy to the expected growth rate under the glucose-only policy in Fig. 2B (the detailed calculation of these ratios is described in ‘Materials & Methods’). To constrain our choice of growth rates, we analyzed growth measurements of 61 yeast strains cultured with different sugars as primary carbon sources (see ‘Materials & Methods’ for details). As expected, median growth rate in glucose was higher than in other sugars (Fig. S1A). Across strains, growth rate in glucose was on average ∼1.5 fold higher than in galactose and in some strains over 3-fold higher (Fig. S1B).

We find that in environments where glucose yields a significantly larger (∼2–4 fold) growth rate than galactose, the posterior predictive policy outperforms the glucose-only policy only in specific regimes of the space of possible Markov environments (red regions in Fig. 2B). For a wide range of environments, the non-adaptive policy is equal to or better than the adaptive policy. As expected, when the difference in growth rate between glucose and galactose gets smaller (Fig. 2B, left to right), the payoff from using the posterior predictive policy is greater. However, as our analysis of growth rates in yeast strains showed, glucose typically confers a substantially higher growth rate than galactose.

This idealized calculation shows the importance of the probabilistic structure of the environment in assessing where an adaptive strategy would pay off. This suggests that in order to highlight the advantage of adaptive strategies, specific types of environmental fluctuations will have to be used.

Our analysis of Markov nutrient environments above doesn’t capture several key aspects of adaptive growth in changing environments. First, our environment’s changes had a simple “flat” structure describable by only two parameters, whereas natural environment may be generated by far more complex underlying mechanisms. Second, our comparison of the adaptive and glucose-only strategies assumed that the transition probabilities governing the environment, θGlu→Glu and θGal→Glu, are known and can be used by cells. In reality, if this information can be used by cells, it has to be learned from the environment and cannot be assumed as given. Third, such information has to be learned in real-time as cells must respond to the environment while it is changing. We now address each of these aspects in turn.

Changes in complex environments may be governed by dynamic processes that are unobservable to cells. Some environments may oscillate between noisy regimes, where the nutrient switches are less predictable, and stable regimes where the nutrient switches are either rare or more predictable. Such environments can be thought of as “meta-changing” in the sense that there’s a change in the way they fluctuate: the probability of being a specific state of the environment (e.g., where a nutrient is available) changes through time, depending on an unobserved condition, such as whether we’re in a noisy or stable regime. The transition from noisy to stable regimes might itself be governed by another time-varying mechanism. As an intuitive example of meta-changing environments, consider the eating routine of animals like us. During feeding, bursts of nutrients that are otherwise scarce may become available in the gut. The fluctuations in nutrient levels within a feeding period will depend on what and how much is being consumed. The separation between meals is also subject to randomness, but can still be predictable, depending on how consistent we are in our eating schedule. This high-level structure can be exploited by adaptive systems to anticipate future changes and to separate noisy fluctuations in nutrients from signals of feeding periods.

To understand the adaptive strategies that may be used for effective growth in such environments, we developed a dynamic Bayesian model of meta-changing environments. In our model, we assume a fixed number of hidden “switch states” that correspond to regimes in the environment, and these states are used to generate the fluctuations in nutrients (Fig. 3A—see ‘Materials & Methods’ for full model description.) The switch states and their dynamics are not observable, and therefore have to be learned from the environment. To do this, our model has a prior distribution over the dynamics of the unobservable switch states and the nutrient transitions in the environment (‘Materials & Methods’). Through experience with the environment, these priors are updated using Bayesian inference to learn the dynamics that drive nutrient fluctuations. Formally, the goal is to compute the posterior predictive distribution over nutrients, P(Ct+1∣C0:t), which depends on the nutrient history C0:t and on the hidden switch state St+1: (1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}\begin{eqnarray*}P({C}_{t+1}\mid {C}_{0:t})=\sum _{{S}_{t+1}}P({C}_{t+1}\mid {C}_{t},{S}_{t+1}=i)P({S}_{t+1}\mid {C}_{0:t}).\end{eqnarray*}\end{document}PCt+1∣C0:t= ∑St+1PCt+1∣Ct,St+1=iPSt+1∣C0:t.

For an inference-based strategy to be biologically plausible, it has to be carried out in real-time since cells respond to the environment while it is changing. To compute the posterior predictive distribution (Eq. (1)) in real-time, we use a particle filtering algorithm (Cappé, Godsill & Moulines, 2007) (described in detail in Algorithm S1 and ‘Materials & Methods’). Rather than naively “memorizing” the environment’s history, in particle filtering a representation of the posterior distribution is continuously updated as the environment changes (‘Materials & Methods’). The uncertainty of the distribution is represented by a set of “particles” (which can be thought of as samples from a distribution). Each particle corresponds to a configuration of the hidden states of our system. For our model, each particle p is a set {s, S, C}, where s corresponds to the value of the hidden switch state, S is a transition matrix that tracks the frequencies of transitions between switch states and C is a transition matrix that tracks the frequency of transitions between nutrients (‘Materials & Methods’).

The particle filtering algorithm can be understood by analogy to evolution through mutation and selection. Initially, all particles are weighted equally. Before the environment changes, particles are probabilistically assigned to new configuration based on our model of the environment (“mutation” step). When a new state of the environment is observed, the particles are re-weighted by their fit to this observation and probabilistically resampled using the updated weights (“selection” step). This process repeats as the environment continues to change. Particles that represent more probable states of the environment will get “selected” for through time, while the “mutation” and resampling steps ensure that diversity is brought into the particle population. As each particle p goes through this process, the hidden state and nutrient state configurations in its trajectory are counted by updating the transition matrices S and C. Since the particles are propagated through this process in parallel, particle filtering takes what naively would be a serial computation (requiring the complete nutrient history) and converts it to a parallel one that can be performed in real-time. This property may make particle filtering amenable to implementation by molecular circuits, as discussed later.

Our model makes a number of predictions about the dynamics of adaptation by systems that represent hidden environmental states. In Fig. 3B, a meta-changing environment is shown that switches between two regimes: one with periodic switches between glucose and galactose, and another where glucose is constant. The posterior predictive probability of glucose in the next time step, P(Ct+1∣Ct), as it gets updated by real-time inference, is plotted along the environment (Fig. 3, top).

The change in the posterior predictive distribution has a number of characteristic features. Starting with a uniform probability over the nutrients, the posterior predictive distribution slowly changes to “learn” the first periodic regime of the environment (Fig. 3B, red line). When the environment shifts to the constant regime (at t = 20, Fig. 3B), the posterior distribution also updates slowly to reflect this. However, when the environment shifts again to the periodic regime (at t = 40), the posterior predictive distribution updates faster to reflect the periodicity. Similarly, when the environment shifts to the constant region for the second time (at t = 60), the posterior predictive distribution changes even more quickly, since the hidden switch state where the environment is periodic has been seen before. More generally, our model predicts that the more familiar stretches of the environment will be adapted to more quickly. By contrast, a Markov model with a “flat” structure that only tracks transition probabilities between glucose and galactose does not show this behavior (Fig. 3B, dotted line). This is one of several subtle predictions about the dynamics of adaptation one would expect from an inference-based strategy that uses a representation of hidden environmental states.