Stochastic modelling reveals mechanisms of metabolic heterogeneity

TL;DR

This paper introduces a stochastic model demonstrating that fluctuations in enzyme expression can lead to diverse metabolic states within genetically identical cell populations, revealing new mechanisms behind metabolic heterogeneity.

Contribution

The study presents a novel stochastic modeling approach showing how enzyme fluctuations cause complex metabolic heterogeneity, including multimodal metabolite distributions.

Findings

Metabolic heterogeneity can arise from enzyme expression fluctuations.

Clonal populations can split into distinct metabolic subpopulations.

Metabolite distributions can be bimodal or multimodal despite unimodal enzyme expression.

Abstract

Phenotypic variation is a hallmark of cellular physiology. Metabolic heterogeneity, in particular, underpins single-cell phenomena such as microbial drug tolerance and growth variability. Much research has focussed on transcriptomic and proteomic heterogeneity, yet it remains unclear if such variation permeates to the metabolic state of a cell. Here we propose a stochastic model to show that complex forms of metabolic heterogeneity emerge from fluctuations in enzyme expression and catalysis. The analysis predicts clonal populations to split into two or more metabolically distinct subpopulations. We reveal mechanisms not seen in deterministic models, in which enzymes with unimodal expression distributions lead to metabolites with a bimodal or multimodal distribution across the population. Based on published data, the results suggest that metabolite heterogeneity may be more pervasive than previously thought. Our work casts light on links between gene expression and metabolism, and provides a theory to probe the sources of metabolite heterogeneity.

Tables

Table 1. Table 1: Nominal parameters for the stochastic model. Parameters correspond to the simulations of Figure 1 A. We use realistic enzyme kinetic parameters Bar-Even  et al. ( 2011 ) and fast promoter switching according to measured ranges So  et al. ( 2011 ) . The dilution rate δ 𝛿 \delta constant corresponds to a doubling rate of approximately 46 minutes, typical in the Escherichia coli bacterium. The distributions in Figure 1 C were obtained with perturbed promoter switching parameters k on = 0.01 ​  s − 1 subscript 𝑘 on 0.01 superscript  s 1 k_{\text{on}}=0.01\text{ s}^{-1} and k off = { 0.03 ,  0.01 ,  0.0001 } ​  s − 1 subscript 𝑘 off 0.03 0.01 0.0001 superscript  s 1 k_{\text{off}}=\left\{0.03,\,0.01,\,0.0001\right\}\text{ s}^{-1} . The parameter values for Figures 2 and 3 are shown in the respective captions.

	value	unit
[]$n_{\text{s}}$	3000	molecule
$k_1$	1	${\text{s}}^{-1}{\text{molecule}}^{-1}$
[]$k_{-1}$	1000	${\text{s}}^{-1}$
$k_{\text{cat}}$	3.6	${\text{s}}^{-1}$
[] $k_{\text{rev}}$	0.01	${\text{s}}^{-1}{\text{molecule}}^{-1}$
$k_{\text{c}}$	0.02	${\text{s}}^{-1}$

Equations

$$\mathbf{P}(n_{\text{p}})=\sum _{n_{\text{etot}}=0}^{\infty }\underset{\text{gene expression}}{\underbrace{\mathbf{P}(n_{\text{etot}})}}\underset{\text{catalysis}}{\underbrace{\mathbf{P}(n_{\text{p}}|n_{\text{etot}})}}.$$

$$\begin{array}{rl}{\text{k}}_{\text{birth}}^{\text{eff}}& ={\mathrm{k}}_{\text{cat}}\mathbb{E}(n_{\text{c}}|n_{\text{etot}},n_{\text{p}})\\ & \approx {\mathrm{k}}_{\text{cat}}\frac{k_{-1}}{(k_{-1}+k_1{n}_{\text{s}})}{n}_{\text{etot}},\\ {\text{k}}_{\text{death}}^{\text{eff}}& =k_{\text{rev}}\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}})+k_{\text{c}}\\ & \approx k_{\text{rev}}\frac{k_1{n}_{\text{s}}}{(k_{-1}+k_1{n}_{\text{s}})}{n}_{\text{etot}}+k_{\text{c}},\end{array}$$

$$\mathbf{P}(n_{\text{p}}|n_{\text{etot}})$$

$$\lambda (n_{\text{etot}})$$

$${\lambda }_{\infty }=n_{\text{s}}\frac{k_{\text{cat}}{k}_1}{k_{\text{rev}}{k}_{-1}},\text{and}K$$

$$\lambda (n_{\text{etot}})=\frac{k_1{n}_{\text{s}}{k}_{\text{cat}}}{(k_1{n}_{\text{s}}+k_{-1})k_{\text{c}}}{n}_{\text{etot}},$$

$$\text{Substrate}+\text{Enzyme}$$

$${\rightleftarrows }_{k_{\text{rev}}}^{k_{\text{cat}}}\text{Metabolite}+\text{Enzyme}$$

$${\to }^{k_{\text{tl}}}\text{mRNA}+\text{Enzyme}$$

$${\text{DNA}}_{\text{on}}$$

$${\text{DNA}}_{\text{off}}$$

$${\to }^{k_{\text{c}}}\varnothing$$

$${\to }^{k_{\text{deg}}}\varnothing$$

$${\to }^{\delta }\varnothing$$

$$\mathbf{P}(n_{\text{etot}})$$

$${\left(\frac{b}{1+b}\right)}^{n_{\text{etot}}}{\left(1-\frac{b}{1+b}\right)}^{{\alpha }_{+}}\times$$

$$\text{textsubscript}2F_1\left({\alpha }_{+}+n_{\text{etot}},\gamma -{\alpha }_{-},\gamma +n_{\text{etot}};\frac{b}{1+b}\right),$$

$$\varnothing {\rightleftarrows }_{{\text{k}}_{\text{rev}}^{\text{eff}}}^{{\text{k}}_{\text{cat}}^{\text{eff}}}\text{Metabolite}{\to }^{{\mathrm{k}}_{\text{c}}}\varnothing ,$$

$$\begin{array}{rl}{\text{k}}_{\text{cat}}^{\text{eff}}& ={\mathrm{k}}_{\text{cat}}\mathbb{E}(n_{\text{c}}|n_{\text{etot}},n_{\text{p}}),\\ {\text{k}}_{\text{rev}}^{\text{eff}}& =k_{\text{rev}}\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}}),\end{array}$$

$$\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}})+\mathbb{E}(n_{\text{c}}|n_{\text{etot}},n_{\text{p}})=n_{\text{etot}}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}})=$$

$$-k_1{n}_{\text{s}}\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}})+k_{-1}\mathbb{E}(n_{\text{c}}|n_{\text{etot}},n_{\text{p}})$$

$$+k_{\text{cat}}\mathbb{E}(n_{\text{c}}|n_{\text{etot}},n_{\text{p}})-k_{\text{rev}}\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}})n_{\text{p}}.$$

$$\begin{array}{rl}\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}})& \approx \frac{k_{-1}}{(k_{-1}+k_1{n}_{\text{s}})}{n}_{\text{etot}},\\ \mathbb{E}(n_{\text{c}}|n_{\text{etot}},n_{\text{p}})& \approx \frac{k_1{n}_{\text{s}}}{k_{-1}+k_1{n}_{\text{s}}}{n}_{\text{etot}},\end{array}$$

$$\frac{{\lambda }_{\infty }}{K}=\frac{k_{\text{cat}}}{k_{\text{c}}}\times \underset{\text{saturation }ϵ}{\underbrace{\frac{k_1{n}_{\text{s}}}{k_1{n}_{\text{s}}+k_{-1}}}},$$

Full text

Introduction

Cellular heterogeneity is ubiquitous across all domains of life. In microbes, clonal populations display phenotypic variability as a result of multiple factors such as fluctuations in the microenvironment, stochasticity in gene expression, or asymmetric partitioning at cell division Elowitz, Levine, and Siggia (2002); Paulsson (2005); Raj and van Oudenaarden (2008). Variability is well recognised at the transcriptional and translational levels. Yet various single-cell phenomena result from the emergence of distinct metabolic states within a clonal population. For example, metabolic heterogeneity plays a key role in antibiotic toleranceBalaban et al. (2004); Lewis (2007); Shan et al. (2017), heterogeneous nutrient uptake Vilhena et al. (2018); Nikolic et al. (2017), and variations in growth rate Kiviet et al. (2014); Thomas et al. (2018). It has also been shown that nutrient shifts can cause populations to split into two van Heerden et al. (2014); Kotte et al. (2014) or more Şimşek and Kim (2018) subpopulations with distinct growth abilities. The emergence of subpopulations has been theorised as a bet-hedging strategy that gives an evolutionary advantage for survival in adverse environments Balaban et al. (2004); Acar, Becskei, and van Oudenaarden (2005).

A central challenge to quantify metabolic variability is the lack of techniques for measuring metabolites with single-cell resolutionTakhaveev and Heinemann (2018). In contrast to single-cell measurements of protein expression, for which powerful reporter systems have been developed Golding et al. (2005); Taniguchi et al. (2010), quantification of metabolites in single-cells remains a major challenge. Some of the techniques employed so far include Förster resonance energy transfer (FRET) sensors Lemke and Schultz (2011), metabolite-responsive transcription factors Xiao et al. (2016); Mannan et al. (2017), RNA sensors Paige et al. (2012), and mass-spectrometry Ibanez et al. (2013), yet most of these technologies are in early stages of developmentTakhaveev and Heinemann (2018). As a result, metabolic heterogeneity is typically quantified indirectly via measurements of metabolic enzymes or growth rate in single-cells Ozbudak et al. (2004); Kotte et al. (2014); Kiviet et al. (2014).

Our objective in this paper is to characterise heterogeneity in metabolites as a result of stochastic enzyme expression and catalysis. Metabolic models traditionally assume that enzymatic reactions behave deterministically on the basis that both enzymes and metabolites appear in high molecule numbers Cornish-Bowden (2004). However, single-cell proteomics in Escherichia coli show that metabolic enzymes are as variable as any other member of the proteome Taniguchi et al. (2010), while metabolomics data suggest that average metabolite abundances span several orders of magnitudeBennett et al. (2009). The few datasets on single-cell metabolite abundance already suggest substantial variability in some metabolites in E. coli Yaginuma et al. (2014); Xiao et al. (2016). Such evidence casts doubt on the traditional assumption of metabolism being a purely deterministic process, suggesting a link between fluctuations in enzyme expression and metabolites.

The role of stochastic gene expression in protein variability has been well studied Paulsson (2005); Raj and van Oudenaarden (2008); Shahrezaei and Swain (2008); Labhsetwar et al. (2013); Thomas, Popović, and Grima (2014), but the impact of such randomness on metabolic reactions remains much less understood. Various theoretical studies have analysed the impact of fluctuations in the supply and consumption of metabolites Levine and Hwa (2007); Thomas, Straube, and Grima (2011); Gupta, Milias-argeitis, and Khammash (2017), or the propagation of enzyme noise to a metabolite Oyarzún, Lugagne, and Stan (2015). However, despite the advancing experimental evidence of stochasticity in metabolism, mathematical models still lack the sufficient detail to integrate the processes that are known to shape protein heterogeneity, such as stochastic promoter switching and transcriptional bursting.

In this paper we propose a model for metabolite heterogeneity in single-cells. The model integrates stochasticity in enzyme catalysis Cornish-Bowden (2004) and expression Shahrezaei and Swain (2008), two well-established processes that so far have been studied in isolation. Our approach includes a stochastic formulation of various relevant mechanisms in enzymatic reactions, including reversible catalysis, stochastic switching of promoter activity, fluctuations in mRNA transcripts, and consumption of the enzymatic product by downstream processes.

We probe the model for various sources of stochasticity using simulations and analytical solutions for the stationary distribution of the metabolite. The analysis reveals intricate patterns of heterogeneity that translate into bimodal and multimodal distributions for the number of metabolite molecules. These phenomena arise from the interplay between a lowly abundant enzyme and its catalytic parameters. Under the separation of timescales typical of metabolic reactions, we show that metabolite distributions can be accurately approximated by a Poisson Mixture Model across large regions of the parameter space. The mixture model can be readily adapted to a wide class of gene expression models and provides a quantitative tool to predict metabolite variability from enzyme measurements in single-cells.

Results

Stochastic model of an enzymatic reaction

We consider a model that combines enzyme kinetics and enzyme expression into a single stochastic description (Figure 1A). The model includes an enzymatic reaction with standard Michaelis-Menten kinetics, in which substrate and enzyme bind reversibly to form a complex that undergoes reversible catalysis into a metabolite. We assume that enzyme expression follows the well established three-stage model for gene expression Raj and van Oudenaarden (2008); Shahrezaei and Swain (2008), where a single copy gene switches stochastically between an inactive state ($D_{\text{off}}$) and active state ($D_{\text{on}}$). In the active state, mRNAs are transcribed and translated into protein. The model also includes consumption of the metabolite by downstream pathways, degradation of mRNA transcripts, and dilution by growth of all species. Since metabolic reactions operate far from thermodynamic equilibrium, we assume that the substrate pool remains constant so that the system reaches a non-zero flux, e.g. when the substrate is a highly abundant extracellular carbon source or a slowly varying intracellular metabolite. The model reactions are shown in equations (R1)–(R9) in the Methods section.

To investigate the emergence of metabolic heterogeneity, we need to compute the stationary probability distribution of metabolite molecules ($n_{\text{p}}$) for relevant combinations of model parameters. Figure 1B shows a typical simulation of the model obtained with Gillespie’s algorithm Gillespie (1976). A key challenge for such simulations, however, is the multiscale nature of enzymatic reactions: not only do metabolic reactions operate in a much faster timescale (milliseconds) than enzyme expression (tens of minutes)Levine and Hwa (2007); Cao, Gillespie, and Petzold (2005); Lugagne, Oyarzún, and Stan (2013), but also the average number of enzymes is much lower than the number of metabolites. These multiple scales result in reaction propensities that differ by several orders of magnitude, thus leading to extremely slow simulations which make the exploration of the parameter space infeasible. An alternative is to use simulation algorithms that exploit the separation of scales to increase computational speed, such as tau-leaping or slow-scale approximationsGillespie (2007). Yet in our case it is unclear how such numerical approximations impact the predictions drawn from the simulations.

To determine the impact of genetic and catalytic parameters on metabolic heterogeneity, we obtained an analytic approximation for the distribution of metabolite molecules that can be evaluated efficiently without expensive stochastic simulations. Our solution allows the exploration of parameter space to characterise the different regimes promoting metabolic heterogeneity. The approximation follows from exploiting time scale separation in the Chemical Master Equation of the stochastic process van Kampen (1992). In physiological regimes the model has three timescales: a fast metabolic time scale, in which substrate and enzyme bind and unbind; an intermediate time scale associated with the catalysis of the metabolite ($n_{\text{p}}$); and a slow timescale associated with the expression of the enzyme and dilution by cell growth.

The total amount of enzyme (free and substrate-bound, denoted as $n_{\text{e}}$ and $n_{\text{c}}$ respectively) varies in the slowest timescale, and therefore the binding/unbinding of substrate and enzyme equilibrates quickly. As a result, in the timescale of gene expression, the metabolite can be assumed to depend directly on the the total enzyme $n_{\text{etot}}=n_{\text{e}}+n_{\text{c}}$ rather than on $n_{\text{e}}$ and $n_{\text{c}}$ individually. Under this approximation, it is convenient to use the law of total probability:

[TABLE]

The formula in (1) decomposes the distribution of metabolite $\mathbf{P}(n_{\text{p}})$ into stochasticity originating from enzyme expression, $\mathbf{P}(n_{\text{etot}})$, and from fluctuations in the catalytic reaction itself, described by the conditional distribution of metabolite given the amount of total enzyme, $\mathbf{P}(n_{\text{p}}|n_{\text{etot}})$. In the timescale of metabolite fluctuations, the total enzyme can be assumed to be in a quasi-stationary state. Further, exploiting the fast binding/unbinding between substrate and enzyme, we showed that the metabolite follows a birth-death process with effective propensities (details in Methods section):

[TABLE]

where $\mathbb{E}(n_{\text{e}}|n_{\text{etot}},n_{\text{p}})$ and $\mathbb{E}(n_{\text{c}}|n_{\text{etot}},n_{\text{p}})$ are the conditional expectations of the free enzyme ($n_{\text{e}}$) and complex ($n_{\text{c}}$) given the total enzyme and metabolite. In equation (2), $n_{s}$ is the constant number of substrate molecules, the parameters $k_{1}$, $k_{-1}$, $k_{\text{cat}}$, and $k_{\text{rev}}$ are the rate constants of the Michaelis-Menten mechanism (defined in Figure 1A), and $k_{\text{c}}$ is an effective first-order rate constant of metabolite consumption by downstream pathways. The conditional distribution needed in the model (1) can then be computed explicitly:

[TABLE]

with Poisson parameter

[TABLE]

and $(\lambda_{\infty},\,K)$ are two effective kinetic parameters

[TABLE]

The parameters $\lambda_{\infty}$ and $K$ are in units of molecules$/$cell and depend on the interplay between substrate abundance, enzyme kinetics, and downstream processes.

As illustrated in Figure 1B, the distribution in (1) is a Poisson Mixture ModelChaturvedi et al. (1977); Iyer-Biswas, Hayot, and Jayaprakash (2009); Dattani and Barahona (2017) (PMM) that convolves the enzyme distribution $\mathbf{P}(n_{\text{etot}})$ with various Poisson modes $\mathbf{P}(n_{\text{p}}|n_{\text{etot}})$ arising from the catalytic activity. In our model, the analytical distribution of the total enzyme abundance follows the standard solution of the three-stage model for gene expression Shahrezaei and Swain (2008), which can be computed explicitly in terms of model parameters. In certain limits, the three-stage model produces approximately Gamma or Normal distributions depending on the the mean expression level and the half lives of mRNAs and proteinsRaj and van Oudenaarden (2008); Shahrezaei and Swain (2008).

The decomposition in equation (1) shows that the PMM is not limited to the model for gene expression we have considered here. Other models may be used, either by using closed-form expressions for $\mathbf{P}(n_{\text{etot}})$, or by inferring the enzyme distribution directly from single-cell protein expression data such as flow cytometry or single-cell microscopy Taniguchi et al. (2010); Elowitz, Levine, and Siggia (2002). The PMM thus provides a versatile tool to predict metabolite heterogeneity from modelled or measured enzyme heterogeneity viewed as an upstream source of variationDattani and Barahona (2017). Figure 1C shows that the PMM distribution provides a good approximation to Gillespie simulations computed with typical parameter values.

Qualitative features of the Poisson Mixture Model

At the heart of the metabolite PMM is the interplay between variability from gene expression and that originating from enzyme kinetics. Specifically, the Poisson parameter $\lambda(n_{\text{etot}})$ in equation (4) controls the density and dispersion of the Poisson modes, which in turn shape the overall pattern of variability. As shown in Figure 1B, there are several cases of interest. For example, for irreversible reactions ($k_{\text{rev}}=0$), the Poisson parameter simplifies to

[TABLE]

which scales linearly with the enzyme abundance and thus the Poisson modes have equidistant means. In reversible reactions, on the other hand, the Poisson parameter saturates and causes the Poisson modes to concentrate around $\lambda_{\infty}$. This effect is stronger for strong reversibility (high $k_{\text{rev}}$), in which case the kinetic parameter $K$ is small. Note also that in either case, as the enzyme number $n_{\text{etot}}$ grows, the Poisson modes spread out since $\lambda(n_{\text{etot}})$ controls both their mean and variance.

From the construction of the PMM in (1), we observe that the enzyme distribution weighs the various Poisson modes, potentially producing metabolite distributions that are unimodal, bimodal or even multimodal. For example, for highly expressed reversible enzymes, the distribution $\mathbf{P}(n_{\text{etot}})$ is non-negligible when $n_{\text{etot}}$ is large. Hence most Poisson modes do not contribute to the final metabolite distribution, except the mode centred at $\lambda_{\infty}$, which leads to a unimodal metabolite distribution with a mean close to the deterministic average.

Conversely, for lowly expressed enzymes, there is a non-negligible probability of enzymes not being expressed, and thus the first term of the PMM, i.e. $\mathbf{P}(0)\text{Poisson}(n_{\text{p}},0)$, causes the metabolite distribution to peak at zero. However, the metabolite distribution may also display a second peak at $\lambda_{\infty}$ if, for example, the $\lambda(n_{\text{etot}})$ parameter causes many Poisson modes to concentrate around $\lambda_{\infty}$. This results in a bimodal metabolite distribution, whereby an isogenic population splits into metabolite producers and non-producers. Similar reasoning can be used to understand the emergence of multimodal metabolite distributions, which correspond to three or more subpopulations with varying metabolic activities. This qualitative analysis suggests that metabolic subpopulations can emerge even in cases where enzymes display unimodal distributions across the population. Crucially, this also indicates that metabolic subpopulations emerge through mechanisms that do not follow trivially from transcriptional heterogeneity alone, as we explore in more detail in next section.

Mechanisms for metabolic bimodality

First, we explored the impact of stochastic promoter switching on the emergence of metabolite bimodality. Figure 2A shows the summary of calculations when evaluating the PMM for variations in two key characteristics of the promoter across several orders of magnitude: the switching time scale ($(k_{\text{on}}+k_{\text{off}})/\delta$) and the promoter activity ($k_{\text{on}}/k_{\text{off}}$), for various values of the kinetic parameter $\lambda_{\infty}$. We found three qualitatively distinct parameter regimes for the metabolite distribution that emerge from the combination of stochastic switching and catalysis: (1) a regime where both enzyme and metabolite have unimodal distributions, akin to the results shown earlier in Figure 1C; (2) a regime where both enzyme and metabolite have bimodal distributions; and (3) a regime in which the enzyme is unimodal but the metabolite is bimodal.

It can be shown that the deterministic version of our model in equations (R1)–(R9) has a single steady state. Hence regime (1) can be thought of as a stochastic correction consisting of unimodal distributions around a deterministic steady state. This is the expected behaviour under the traditional assumptions of high abundance of enzyme and metabolite molecules.

The other two regimes, however, correspond to alternative routes of noise-induced bimodality that cannot be explained using deterministic modelsOyarzún and Chaves (2015); Lipshtat et al. (2006); To and Maheshri (2010). Regime (2) is a highly stochastic regime dominated by the slow stochastic switching of the promoter, which drives and entrains the metabolic response. Hence we term it switching-induced bimodality. Slowly switching promoters are known to produce bimodal gene expression Thomas, Popović, and Grima (2014); Dattani and Barahona (2017), and thus this regime corresponds to a case in which bimodality propagates from enzymes to metabolites. Figure 2A shows that this behaviour appears robustly for slow switching and high promoter activity across values of $\lambda_{\infty}$ kinetic parameter.

Regime (3), the second route for metabolite bimodality, originates from a unimodal but weakly expressed enzyme (low $k_{\text{on}}/k_{\text{off}}$) expressed from fast switching promoters. In this case, the birth of a small number of enzyme molecules is sufficient to kick-start catalysis and make it rapidly settle in a quasi stationary regime. This distinct phenomenon is a result of the separation of time scales between enzyme expression and catalysis, and we refer to it as catalytically-induced bimodality. From Figure 2A, we observe that this form of bimodality appears for a narrow range of promoter switching parameters corresponding to fast switching genes with medium to low promoter activity. This behaviour disappears altogether for a low $\lambda_{\infty}$ parameter, for example in case of strong reversibility.

To validate the predictions of the PMM approximation, we ran full Gillespie simulations over a long time horizon for different parameter sets. Figure 2B shows the simulation time courses and resulting histograms. For switching-induced bimodality, we observe how slowly switching promoters causes a single cells to lack the enzyme over several cell cycles, a period during which the metabolite is not produced. In the case of catalytic-induced bimodality, however, fast switching combined with a low average expression level cause the metabolite abundance to drop for shorter but more frequent intervals. In both cases the PMM provides an excellent approximation to the bimodal histograms obtained from the stochastic simulations. Furthermore, we observe that the bimodal metabolite distributions both regimes are almost indistinguishable from each other, yet they are produced by enzymes with substantially different time courses and distributions. These regimes therefore correspond to distinct forms of bimodality, arising from fundamentally different mechanisms.

Emergence of metabolic multimodality

To explore the emergence of multimodality, we examined the analytical formula of the PMM in (1) to identify kinetic regimes associated with distinct enzyme distributions. A necessary condition for the emergence of multiple modes is that the Poisson components do not overlap and are sufficiently spaced from each other. From the definition of the $\lambda(n_{\text{etot}})$ parameter in (5), this happens when the kinetic parameter $K$ is large. As discussed earlier, depending on the distribution of the enzyme, the Poisson modes may appear or cancel in the final metabolite distribution. We thus swept the parameter $K$ and evaluated the PMM across various enzyme expression levels, including low expression with a skewed distribution and high expression with a normally distributed enzyme.

As shown in Figure 3, we found intricate patterns of multimodal distributions, depending on the interplay between the heterogeneity of the enzyme, $\mathbf{P}(n_{\text{etot}})$, and the enzyme kinetics encapsulated by the $K$ parameter. Multimodality appears when the enzyme expression levels are low as compared to the parameter $K$. For instance, the values of $K$ in Figure 3 are approximately 5-, 20-, and 100-fold those used in the bimodal examples in Figure 2. For enzymes expressed at intermediate levels, in the order of tens of molecules/cell on average, we found metabolite distributions that are unimodal but highly skewed. In the case of highly expressed enzymes, metabolites followed approximately normal distributions for a wide range of kinetic parameters.

The predictions are confirmed by Gillespie simulations of the full stochastic model, which display a striking match with the PMM approximation, even for complex multimodal distributions. The simulation time courses (shown in the insets of Figure 3) show that the multiple modes for weakly expressed enzymes correspond to cells remaining in a fixed metabolic state over the scale of the cell cycle but fluctuate across other states over longer time scales. For intermediate enzyme expression and large values for $K$, the metabolite does not settle in the quasi-stationary states and displays a long-tailed distribution. A decrease in $K$ suppresses the tail of the distribution driving the PMM towards an approximately normal distribution. Altogether, these results indicate that the relation between enzyme expression and the kinetic parameters $\lambda_{\infty}$ and, in particular, $K$ are key determinants for the emergence of multimodality. This underscores the utility of the PMM to guide the prediction of qualitative and quantitative features of metabolite distributions for a wide range of parameter combinations.

Discussion

Metabolic reactions are the powerhouse of living systems, fuelling the activity and dynamics of most cellular functions. Yet metabolism has been traditionally considered as a static process isolated from the rest of the cellular machinery. Currently, the accepted notion is that due to the large number of molecules involved, metabolism is a deterministic process at the cellular level, modulated by potentially random extrinsic factorsWehrens et al. (2018); Kotte et al. (2014). Here we integrated enzyme kinetics and enzyme expression to propose a theoretical model for variability of metabolites in single cells. The model suggests that cell-to-cell metabolite variation can also arise as a result of intrinsic sources such as stochastic fluctuations in enzyme expression.

The majority of work on non-genetic heterogeneity has focused on stochastic gene expression and the resulting variability in protein levels Elowitz, Levine, and Siggia (2002); Paulsson (2005). This has produced a wealth of single-cell data and models to understand the variability in transcription and translation observed in clonal populations. Metabolite heterogeneity, however, remains poorly understood theoretically and has been observed only indirectly (e.g., through measurements of metabolic enzymes Ozbudak et al. (2004); Kotte et al. (2014) or growth rate Kiviet et al. (2014)) due to the lack of techniques to measure metabolite abundance in single cells.

Using the separation of time scales characteristic of metabolic reactions, we found that the stationary distribution of a metabolite follows a Poisson Mixture Model (PMM). The PMM can be efficiently evaluated across large domains of the parameter space and provides excellent approximations to the distributions computed from full stochastic simulations. Importantly, the model can be readily adapted to include different stochastic models for enzyme expression, beyond the three-stage model considered here Raj and van Oudenaarden (2008); Corrigan et al. (2016), or even stochastic and time-dependent enzyme expression modelled as upstream drivesDattani and Barahona (2017). The model can also be parameterised from experimentally measured distributions for enzyme levels in single-cellsTaniguchi et al. (2010). In combination with the enzyme kinetic parameters, the PMM could provide a powerful tool to predict metabolite variability from single-cell protein data obtained with established techniques such as flow cytometry or time-lapse microscopy.

We found complex patterns of metabolite heterogeneity depending on the interplay between the timescale of promoter activation/deactivation, the enzyme expression level, and the enzyme kinetics. The model predicts that bimodal and multimodal metabolite distributions can emerge in various parameter regimes. In such regimes single-cells spend several cell cycles in a constant metabolic state, but in timescales as long as tens of cell cycles, they switch stochastically across different states. Such long-term fluctuations in single cells result in highly heterogeneous populations containing several subgroups of metabolically distinct cells.

Bimodal metabolic phenotypes have been observed as a result of transcriptional regulation Ozbudak et al. (2004); Acar, Becskei, and van Oudenaarden (2005), post-translational control van Heerden et al. (2014) and stochastic effects triggered by environmental shifts Kotte et al. (2014). Our model reveals two distinct regimes in which metabolites display bimodality. One regime, which we call switching-induced bimodality, corresponds to the intuitive case in which a bimodal enzyme produces a bimodal metabolite. In agreement with previous studies on stochastic gene expression, this type of bimodality appears as a result of slow switching between promoter statesThomas, Popović, and Grima (2014); Ge et al. (2018). In addition, we identified a fundamentally different mechanism of catalytically-induced bimodality, in which a unimodal enzyme produces a bimodal distribution of metabolite. This phenomenon results from a combination of slow fluctuations of a weakly expressed enzyme and the comparatively faster timescale of enzyme catalysis. Catalytic timescales are typically in the order of seconds or faster, so that slow fluctuations in enzyme expression levels produce two quasi-stationary metabolic states in single cells. At a population level, this leads to two distinct subpopulations of metabolite producers and non-producers.

As shown in Figure 4A, single-cell measurements in E. coli suggest that metabolic enzymes appear in low copy numbers across most cellular pathwaysTaniguchi et al. (2010). In the specific growth conditions of that experiment, the data did not reveal bimodal expression of enzymes, which precludes the emergence of switching-induced bimodality in the metabolites they catalyse. However, as illustrated by the three representative distributions in Figure 4A, a number of enzymes have a low mean and a long-tailed distribution, akin to those required for catalytically-induced bimodality and multimodality. This suggests that enzyme distributions found in nature have the characteristics needed for the emergence of subpopulations with two or more distinct metabolite abundances.

Further requirements for metabolite bimodality and multimodality involve conditions on the parameters $\lambda_{\infty}$ and $K$ in equation (5). However, their computation requires rate constants ($k_{1}$, $k_{-1}$ and $k_{\text{rev}}$) that are rarely measured separately, and instead enzymology data typically provides values for $k_{\text{cat}}$ and $K_{\mathrm{M}\,=}(k_{\text{cat}}+k_{-1})/k_{1}$ only Schomburg et al. (2013). In the Methods we show that the ratio $\lambda_{\infty}/K$ can be expressed as $\lambda_{\infty}/K=\epsilon\times k_{\text{cat}}/k_{\text{c}}$, where $\epsilon$ is the saturation level of the enzyme and $k_{\text{c}}$ is the first-order rate constant of metabolite consumption. As illustrated in Figure 4B, the ratio $\lambda_{\infty}/K$ corresponds to a straight line in a $(\lambda_{\infty},\,K)$-space, and a specific enzyme (i.e. with specific values for $k_{1}$, $k_{-1}$ and $k_{\text{rev}}$) corresponds to a single point on the line. In Figure 4B we compare model predictions for a lowly abundant enzyme with different $\lambda_{\infty}/K$ ratios computed for $k_{\text{cat}}$ constants measured in E. coliBar-Even et al. (2011). Considering the large spread in measured $k_{\text{cat}}$ values, of up to seven orders of magnitude, plus the multiple combinations of metabolite consumption rates and enzyme saturation, the analysis suggests that catalytically-induced bimodality and multimodality are plausible within physiological regimes. Further validations of our predictions require measuring metabolite distributions directly, but this is still subject to a number of challenges in single-cell measurement technologies Takhaveev and Heinemann (2018).

Our analysis shows that metabolite heterogeneity depends on a delicate interplay between enzyme expression and enzyme kinetics. It is reasonable to expect that energy-critical enzymes, such as those in central carbon metabolism, filter away fluctuations through post-translational regulatory mechanisms commonly found in metabolism. However, this may not be the case in pathways that are dynamically regulated in response to changes in the environment or cellular context. For example, transcriptional regulation in response to nutrient shifts may steer enzyme levels into regimes of low copy numbers where heterogeneity may dominate the resulting phenotypes. Such mechanism has been already shown to produce growth bimodality in the gluconeogenic switch of E. coliKotte et al. (2014), while a similar mechanism could underpin the large variability observed in single-cell measurements of S-Adenosyl methioninePaige et al. (2012). Noise-induced phenomena also have implications for the design of dynamic control systems for heterologous pathways, which are focus of much research in synthetic biology and metabolic engineering Liu et al. (2018).

In our efforts to build a theory that includes components shared by most enzymatic reactions, we have purposely overlooked a number of processes that can shape metabolic activity. For example, we have not addressed the impact of feedback mechanisms that control enzyme activity, including e.g. product inhibition and allostery, or transcriptional mechanisms that control enzyme expression in response to metabolites. Since post-translational regulation operates on timescales much shorter than enzyme expression, with similar timescale separation arguments it should be possible to express the metabolite distribution as a mixture model akin to ours. In such case the mixture components are not necessarily Poisson and their distribution will depend on the particular mechanism under study. We expect that bimodal and multimodal responses are likely to emerge in this setting, but the precise parameter conditions would have to be studied on a case-by-case basis. Transcriptional feedback control can also display various mechanisms depending on the particular pathway under study. One common motif relies on transcription factors (TF) that up- or down-regulate enzyme expression upon binding to a specific metaboliteMannan et al. (2017). These mechanisms have been shown to play important roles on metabolic activityChaves and Oyarzún (2019), but they also bring to the fore subtle questions that require detailed examination, for example, on the role of fluctuations coming from TF expression itself, or the impact of negative TF autoregulationFang et al. (2017). Our study paves the way for these and other questions to be addressed and raises exciting prospects for the future research in metabolic heterogeneity.

In this paper we laid theoretical foundations to study metabolism in conjunction with stochastic enzyme expression. We brought together classic models for gene expression and enzyme kinetics, and discovered a rich array of distinct stochastic phenomena that underpin the emergence of metabolic subpopulations. Our theory provides a quantitative basis to draw testable hypotheses on the sources of metabolite heterogeneity, which together with the ongoing efforts in single-cell metabolite measurements, will help to re-think metabolism as an active source of phenotypic variation.

Methods

Stochastic modelling and simulation

We built a fully stochastic model for the reaction scheme describing a metabolic reaction coupled with gene expression (Fig. 1A):

[TABLE]

All reactions are assumed to follow mass action kinetics. Model simulations were computed with Gillespie’s algorithm Gillespie (1976) over long time horizons, in the order of hundreds of cell cycles for all simulations. Because of the complex multimodality observed, long simulations are needed to obtain accurate approximations of the stationary molecular distributions. The time courses shown in figures correspond to a small time window of the overall simulation. Unless mentioned in the figure captions, all parameter values were fixed to their nominal values shown in Table 1. The parameters are selected in a physiologically realistic range respecting the scale separation of molecule numbers between mRNA ($\sim 1-5$ molecules), total enzyme abundance ($\sim 100$ molecules) and metabolites ($\sim 1000$ molecules).

To identify bimodality in Figure 2, we detected the existence of one or two peaks in a distribution and defined it as bimodal if the height of the smaller peak is a least 10% of the larger peak and the trough between peaks is at most 10% of the height of the smaller peak.

Analytical expressions for the metabolite distribution

To derive an analytic approximation for the probability to observe $n_{\text{p}}$ metabolites in a cell, we first use the law of total probability as shown in (1).

Distribution of the total enzyme —

Because free enzymes and complexes degrade at the same rate and $n_{\text{etot}}=n_{\text{e}}+n_{\text{c}}$ is conserved by the metabolic reaction, in the slow timescale the enzyme distribution $\mathbf{P}(n_{\text{etot}})$ follows the standard solutionShahrezaei and Swain (2008) of the three-stage model for gene expression:

[TABLE]

where $\Gamma$ is the Gamma function and $\textsubscript{2}F_{1}$ is the ordinary hypergeometric function. The parameters are $\gamma=(k_{\text{on}}+k_{\text{off}})/\delta$ and $\alpha_{\pm}=\left(a+\gamma\pm\sqrt{(a+\gamma)^{2}-4ak_{\text{off}}/\delta}\right)/2$, with $a=k_{\text{tx}}/\delta$ and $b=k_{\text{tl}}/k_{\text{deg}}$.

Conditional distribution for the metabolite —

To compute the second term in (1), we observe that enzyme expression occurs on a much longer timescale than enzyme kinetics, and thus metabolites can be considered to be in a quasi-equilibrium state of the catalytic reactions (R1)–(R2) and metabolite consumption (R6).

To explicitly compute the mixture components $\mathbf{P}(n_{\text{p}}|n_{\text{etot}})$, we assume that reversible binding between substrate and enzyme in reaction (R1) is much faster than the catalytic step and metabolite consumptionCao, Gillespie, and Petzold (2005). In this limit, the metabolite number evolves according to the effective reactions:

[TABLE]

where $\text{k}_{\text{cat}}^{\text{eff}}$ and $\text{k}_{\text{rev}}^{\text{eff}}$ are effective propensities averaged over the fast fluctuating variables $n_{\text{c}}$ and $n_{\text{e}}$:

[TABLE]

where $\mathbb{E}$ denotes the expectation operator. The derivation of the effective propensities in (9) corresponds to a particular case of a more general methodology for timescale separation in stochastic chemical systemsGoutsias (2005); Thomas, Popović, and Grima (2014); Melykuti, Hespanha, and Khammash (2014). Note that since the total enzyme levels are conserved in the catalytic timescale, it follows that

[TABLE]

To derive the conditional expectations in (9), we write the first-order moment equation for the free enzyme $n_{\text{e}}$, which according from equations (R1)–(R2) reads

[TABLE]

Under the assumption that the reversible binding of substrate and enzyme is much faster than the other processes, the first two terms dominate the right hand side of equation (11) and determine the enzyme-complex quasi-equilibrium. Equating these two terms and using the conservation relation in (10), we obtain

[TABLE]

and thus both conditional expectations depend on $n_{\text{etot}}$ and are independent of the metabolite abundance. Therefore, the reactions in (8) correspond to a birth-death process with a zero-th order birth propensity $\text{k}_{\text{cat}}^{\text{eff}}$ and two linear death propensities. The mixture components $\mathbf{P}(n_{\text{p}}|n_{\text{etot}})$ are thus Poissonian with parameter $\lambda(n_{\text{etot}})$ as shown in (3)–(4).

Comparison of PMM predictions and measured kinetic parameters

The PMM depends on the effective parameters $\lambda_{\infty}$ and $K$, which are functions of five rate constants ($k_{\text{cat}}$, $k_{1}$, $k_{-1}$, $k_{\text{c}}$ and $k_{\text{rev}}$). Most of these parameters are not available, except $k_{\text{cat}}$ and $K_{\text{M}}=(k_{\text{cat}}+k_{-1})/k_{1}$. From equation (5) it follows

[TABLE]

which allows the computation of $\lambda_{\infty}/K$ for measured $k_{\text{cat}}$ values in different saturation conditions and consumption rate constants. The red line in Figure 4B represents the $\lambda_{\infty}/K$ ratio for a saturated enzyme ($\epsilon=1$), fast consumption ($k_{\text{c}}=100\times\delta$) and the median $k_{\text{cat}}\approx 16.5\text{ s}^{-1}$ in E. coli Bar-Even et al. (2011). The top grey line is the case without consumption, i.e. metabolites are diluted by cell growth ($k_{\text{c}}=\delta$). Lower saturation ($\epsilon=0.2$) moves the red line down the vertical axis (bottom grey line). The enzyme distribution in Figure 4B was produced with promoter switching parameters $\left\{k_{\text{on}},\,k_{\text{off}}\right\}=\left\{1.56,\,3\right\}\times 10^{-4}\text{ s}^{-1}$, $k_{\text{tx}}=0.025\text{ s}^{-1}$, and $k_{\text{tl}}=0.2\text{ s}^{-1}$. The boundaries between unimodal, bimodal and multimodal distributions were computed as follows. Unimodal distributions are those with a single maximum. Bimodal distributions were detected as in Figure 2. Multimodal distributions are those with at least one additional peak higher than a threshold of $1\times 10^{-4}$ and the trough between neighbouring peaks at most 90 % of its height.

Acknowledgements

This work was funded by the Human Frontier Science Program through a Young Investigator Grant awarded to D.O. (RGY0076-2015), The Royal Commission for the Exhibition of 1851 through a Fellowship to P.T., and the EPSRC Centre for Mathematics of Precision Healthcare (EP/N014529/1).

Author contributions

M.T., P.T., M.B and D.O. conceived the study. M.T. carried out the theoretical derivations, simulations and analysed the data. M.B. and D.O. supervised the work. M.T., P.T., M.B. and D.O. wrote the paper.

Competing interests

The authors declare no competing interests.

Code availability

The code for producing model simulations is available from the authors upon request.

Data availability

All data is available from the authors upon request.