Statistical mechanics for metabolic networks during steady-state growth

TL;DR

This paper extends flux balance analysis to a maximum entropy framework at the single-cell level, accurately modeling E. coli metabolism and predicting flux variability, correlations, and growth rate distributions validated by microfluidic data.

Contribution

It introduces a maximum entropy model for metabolic fluxes that generalizes FBA, providing more accurate predictions and insights into single-cell metabolic variability and regulation.

Findings

E. coli metabolism is near, but not at, FBA optimality.

Model accurately predicts flux variability and correlations.

Scaling relations for growth rate distributions are validated by experiments.

Abstract

Which properties of metabolic networks can be derived solely from stoichiometric information about the network's constituent reactions? Predictive results have been obtained by Flux Balance Analysis (FBA), by postulating that cells set metabolic fluxes within the allowed stoichiometry so as to maximize their growth. Here, we generalize this framework to single cell level using maximum entropy models from statistical physics. We define and compute, for the core metabolism of Escherichia coli, a joint distribution over all fluxes that yields the experimentally observed growth rate. This solution, containing FBA as a limiting case, provides a better match to the measured fluxes in the wild type and several mutants. We find that E. coli metabolism is close to, but not at, the optimality assumed by FBA. Moreover, our model makes a wide range of predictions: (i) on flux variability, its regulation, and flux correlations across individual cells; (ii) on the relative importance of stoichiometric constraints vs. growth rate optimization; (iii) on quantitative scaling relations for singe-cell growth rate distributions. We validate these scaling predictions using data from individual bacterial cells grown in a microfluidic device at different sub-inhibitory antibiotic concentrations. Under mild dynamical assumptions, fluctuation-response relations further predict the autocorrelation timescale in growth data and growth rate adaptation times following an environmental perturbation.

Figures

G) has been obtained by inverting the previous equation, $\beta=\beta(\overline{x})$, and performing a Legendre transform De Martino *et al.* ([2016](#bib.bib24))

G) has been obtained by inverting the previous equation, $\beta=\beta(\overline{x})$, and performing a Legendre transform De Martino *et al.* ([2016](#bib.bib24))

G) has been obtained by inverting the previous equation, $\beta=\beta(\overline{x})$, and performing a Legendre transform De Martino *et al.* ([2016](#bib.bib24))

Tables

Table 1. Table 1: Average growth rate, empirical maximum, inferred maximum and level of optimization, β ​ λ max 𝛽 subscript 𝜆 max \beta\lambda_{\rm max} , from growth rate data.

Step	sample size	$\overline{\lambda }$ [min-1]	${\lambda }_{e,\max }$ [min-1]	${\lambda }_{\max }$ [min-1]	$\beta {\lambda }_{\max }$ [adim.]
$1$	$5204$	$0.022$	$0.029$	$0.034$	$66$
$2$	$3300$	$0.017$	$0.024$	$0.029$	$54$
$3$	$2252$	$0.012$	$0.018$	$0.023$	$45$
$4$	$1464$	$0.0085$	$0.013$	$0.017$	$45$
$5$	$1536$	$0.007$	$0.011$	$0.015$	$42$

Equations

$${\dot{c}}_{\mu }=\sum _i{S}_{i\mu }{v}_i.$$

$$\sum _i{S}_{i\mu }{v}_i=0,$$

$$v_i\in [v_i^{\mathrm{LB}},v_i^{\mathrm{UB}}].$$

$$P_{\beta }(\mathbf{v})=\{\begin{array}{ll}\frac{1}{Z}{e}^{\beta \lambda (\mathbf{v})}& \mathbf{v}\in \mathcal{P},\\ 0& \mathbf{v}\notin \mathcal{P}.\end{array}$$

$$\bar{\lambda }(\beta )={\int }_{\mathbf{v}\in \mathcal{P}}d\mathbf{v}\lambda (\mathbf{v})P_{\beta }(\mathbf{v})={\lambda }_{\mathrm{data}}.$$

$$S(\beta )=-{\int }_{\mathbf{v}\in \mathcal{P}}d\mathbf{v}{P}_{\beta }(\mathbf{v})\log P_{\beta }(\mathbf{v}).$$

$${\bar{\lambda }}_{\mathrm{eff}}(\beta )=\bar{\lambda }(\beta )-\alpha K\left[2^{2I(\beta )/K}-1\right],$$

$$\frac{\sigma }{{\lambda }_{\max }}\propto (\beta {\lambda }_{\max }{)}^{-1}\propto ({\lambda }_{\max }-\bar{\lambda });$$

$$\dot{p}(\lambda )=(\lambda -\overline{\lambda })p(\lambda )+D\left[\frac{{\partial }^{2}p}{\partial {\lambda }^2}-\frac{\partial }{\partial \lambda }\left[p(\lambda )\frac{\partial }{\partial \lambda }(\log q(\lambda ))\right]\right],$$

$$\sigma \sim {\tau }^{-1}\sim {\lambda }_{\max }-\overline{\lambda }\sim D^{1/3}.$$

$$\tau \sim O\left(D^2\frac{R^2}{r^2}\right),$$

$$q(x)=(a+1)(1-x{)}^a,$$

$$Z(\beta )={\int }_0^{1}q(x)e^{\beta x}dx=$$

$$=\frac{\gamma (a+1,\beta )e^{\beta }(a+1)!}{{\beta }^{a+1}},$$

$$\overline{x}=\frac{d\log Z}{d\beta }.$$

$$I(\overline{x})\log 2=\Delta S(\overline{x})=\overline{x}\beta (\overline{x})-\log Z(\beta (\overline{x})).$$

$${\lambda }_{\mathrm{eff}}(\beta )=\bar{\lambda }(\beta )-\alpha K\left(2^{2I(\beta )/K}-1\right),$$

Full text

††preprint: APS/123-QED

**Significance: **While metabolic reactions have been studied in detail in model organisms and can today in principle be inferred automatically from genome-scale data, our systems-level understanding of the complete metabolic networks rests primarily on flux balance analysis (FBA). In FBA, known reactions impose physical constraints within which free parameters are chosen to maximize the growth rate. Using ideas from statistical physics, we extend this approach to capture cell-to-cell variability. Our model is easy to infer, provides superior fit to measurements in Escherichia coli, and makes several categorically new predictions connecting metabolism to single cell physiology. We find that E. coli grows close to, but not at, optimality, and further predict and experimentally test the scaling of growth rate fluctuations in clonal bacterial populations.

I Introduction

After the significant developments of molecular biology and biochemistry of the last century, many aspects of cellular physiology could be understood as a result of interactions between identified molecular components. Perhaps the best-characterized example is intermediate metabolism, the set of reactions that enable cell growth by converting organic compounds and transducing free energy. Today, one can infer nearly complete metabolic networks from genome-scale data, but the dynamics and parameter dependence of such networks remain difficult to analyze. Alternatively, one can assume that known reactions only provide physico-chemical constraints within which some adaptive dynamics has maximized the growth rate, e.g., by adjusting enzyme levels and controlling reaction rates Kacser et al. (1995). An influential implementation of this idea for batch cultures under steady state conditions has been the flux balance analysis (FBA) Orth et al. (2010), which has been tested experimentally Ibarra et al. (2002); Edwards et al. (2001), also in mutant strains, strains used for industrial production Majewski and Domach (1990); Varma and Palsson (1994); Edwards and Palsson (2000), as well as phenotypes implicated in disease (e.g., Warburg effect Vazquez et al. (2010)). Using maximum entropy ideas from statistical physics, we extend the application of FBA from batch to single-cell level and show that our extension makes a wide range of predictions, some of which we test experimentally.

Recent measurements at the single-cell level demonstrated the existence of substantial cell-to-cell growth rate fluctuations even in well-controlled steady-state conditions Wang et al. (2010). These fluctuations exhibit universal scaling properties Iyer-Biswas et al. (2014a, b); Kennard et al. (2016), relate to cell size control mechanisms Taheri-Araghi et al. (2015), act as a global collective mode for heterogeneity in gene expression Kiviet et al. (2014); Shahrezaei and Marguerat (2015); Keren et al. (2015), and are ultimately believed to affect fitness Cerulus et al. (2016). To link these observations to metabolism, however, we need to set up a mathematical description not only of the optimal metabolic fluxes and maximal growth rate in batch culture (as in FBA, which permits no heterogeneity across cells), but for the complete joint distribution over metabolic fluxes. Metabolic phenotypes of individual cells growing in steady-state conditions can then be understood as samples from this joint distribution, which would automatically contain information about flux correlations, and, in particular, could directly predict cell-to-cell growth rate fluctuations.

The simplest construction of a joint distribution over metabolic fluxes can be derived in the maximum entropy framework Jaynes (1957). The key intuition is to look for the most unbiased (or random) distribution over fluxes through individual metabolic reactions that is consistent with the given stoichiometric constraints, while matching the experimentally measured average growth rate. The maximum entropy model that we specify below will turn out to be a one-parameter family of distributions, where the single parameter can be fit to match experimental data; all subsequent predictions follow directly, without further fitting. A similar approach has recently been used in diverse biological settings, ranging from neural networks Schneidman et al. (2006); Tkačik et al. (2014), genetic regulatory networks Lezon et al. (2006), antibody diversity Mora et al. (2010), and collective motion of starling flocks Bialek et al. (2012).

In addition to accounting for fluctuations, the maximum entropy construction provides a principled interpolation between two extremal regimes of metabolic network function. In the “uniform” (no-optimization) limit, no control is exerted over metabolic fluxes: they are selected at random as long as they are permitted by stoichiometry, resulting in broad yet non-trivial flux distributions that support a small, non-zero growth rate. In the FBA limit, fluxes are controlled precisely to maximize the growth rate, with zero fluctuations. The existence of these two limits defines a fundamental, and still unanswered, question about metabolic networks: Is there empirical evidence that real metabolic networks are located in an intermediate regime between the two limits where fluctuations are non-negligible De Martino et al. (2016), and if so, what are the properties of this intermediate regime? Here, we address this question using metabolic flux and single-cell physiology data for Escherichia coli.

II Results

We start by considering a set of metabolic reactions in the well-mixed, continuum limit. Let $S_{i\mu}$ be the stoichiometric coefficient of the metabolite $\mu$ (whose concentration is $c_{\mu}$) in reaction $i$, whose flux is $v_{i}$. The metabolic network dynamics is then given by mass balance equations:

[TABLE]

Assuming homeostasis (i.e., steady state, $\dot{c}_{\mu}=0$) and including further constraints from thermodynamics, nutrient availability, and kinetic limits in the form of lower (LB) and upper (UB) bounds on fluxes, we obtain a convex polytope $\mathcal{P}$ of feasible steady states (“metabolic phenotypes”) in the space of fluxes:

[TABLE]

In addition to bona fide chemical reactions, constraint-based models often include a phenomenological “biomass reaction” in the form of a linear combination of metabolite fluxes, $\lambda(\mathbf{v})=\sum_{i}\xi_{i}v_{i}$, where the proportions $\xi_{i}$ are set to mimic cell growth, i.e., the metabolite fluxes necessary to reconstitute the biomass of a new cell in a typical division time.

Flux balance analysis looks for the flux configuration $\mathbf{v}_{\rm max}$ that maximizes growth $\lambda_{\rm max}=\lambda(\mathbf{v}_{\rm max})$ subject to constraints given by Eqs (2), which can be easily found by linear programming. In contrast, our maximum entropy approach starts with a distribution over fluxes with a Boltzmann form, which assumes that the fluxes are as random as possible while achieving a desired average growth rate De Martino et al. (2016):

[TABLE]

The parameter, $\beta$, of the distribution $P$ can then be set to match the predicted average growth rate to the measured growth rate, $\lambda_{\rm data}$:

[TABLE]

Alternatively, $\beta$ could also be fit to match (e.g., in a $\chi^{2}$-sense) some subset of fluxes that can be measured experimentally. Once $\beta$ is fixed, the joint distribution of Eq. (3) can be queried for average fluxes, flux correlations, or other quantities of interest that we discuss later.

The maximum entropy distribution with a constrained average growth rate has two interesting limits, as illustrated in Fig 1. The growth rate, $\bar{\lambda}$, increases with $\beta$ (which we will refer to as an “optimization parameter”) until, in the limit $\beta\rightarrow\infty$, the distribution $P_{\infty}(\mathbf{v})$ collapses into a delta function at $\mathbf{v}_{\rm max}$, lying at the boundary of the polytope $\mathcal{P}$: this is the FBA solution that supports the maximal growth rate $\lambda_{\rm max}$. Conversely, as $\beta\rightarrow 0$, Eq (3) yields a uniform sampling of fluxes over the permitted polytope $\mathcal{P}$: this “uniform” solution is an interesting baseline case for comparison because it incorporates all stoichiometric constraints but postulates no regulation-mediated growth rate optimization. In statistical physics, high-$\beta$ regime (limiting toward the FBA solution) corresponds to the “energy-dominated” regime, while the low-$\beta$ regime (limiting toward the uniform sampling) corresponds to the “entropy-dominated” regime; the optimization parameter $\beta$ corresponds to the inverse temperature.

II.1 Flux predictions for Escherichia coli

We constructed a maximum entropy model for the catabolic core of the E. coli metabolic network from the genome scale reconstruction iAF1260 (see SI Appendix), in a glucose-limited minimal medium in aerobic conditions Orth et al. (2011). The network comprises $N=86$ reactions among $M=68$ metabolites and includes glycolysis, pentose phosphate pathway, TCA cycle, oxidative phosphorylation, and nitrogen catabolism. The dimension of the resulting polytope $\mathcal{P}$ of allowed steady states is $D=23$, from which we can efficiently draw flux configurations according to Eq. (3) using Hit-and-Run Monte Carlo Markov Chain after suitable preprocessing De Martino et al. (2015) (see SI Appendix).

To evaluate the predictions of the maximum entropy model, we compared $N_{f}=23$ measured metabolic fluxes in Escherichia coli from previously published data to our predictions, as shown in Fig 2A,C. As a goodness-of-fit measure we defined $\chi^{2}=N_{f}^{-1}\sum_{i=1}^{N_{f}}\left(\langle v_{i}\rangle-V_{i}\right)^{2}/\left(\sigma^{2}_{i}+E^{2}_{i}\right)$, where $V_{i}$ is the measured flux (relative to glucose uptake) and $E_{i}^{2}$ the associated measurement error variance, while $\langle v_{i}\rangle$ and $\sigma_{i}^{2}$ are the mean (and variance, respectively) of the corresponding flux computed in the maximum entropy model of Eq. (3). We chose $\beta$ that minimized $\chi^{2}$, as shown in Fig 2B, to make our flux predictions. $\chi^{2}$ at the optimal parameter $\beta$ was lower than in either the FBA or uniform limits, and this remained robustly true even with alternative goodness-of-fit measures (e.g., mean-squared-error, MSE=$N_{f}^{-1}\sum_{i=1}^{N_{f}}\left(\langle v_{i}\rangle-V_{i}\right)^{2}$ that does not normalize by a variance, $\sigma^{2}_{i}$, which grows as $\beta$ is decreased).

The maximum entropy model thus provided a better account of the measured metabolic fluxes for wild type E. coli at intermediate values of $\beta$, rather than in the limit $\beta\rightarrow\infty$, which corresponds to the flux balance analysis (FBA) solution. Specifically, for both datasets analyzed in Fig 2, the best-fitting dimensionless combination was $\beta^{*}\lambda_{\rm max}\sim 10^{2}$. In addition to a better quantitative match overall, the maximum entropy model correctly predicted non-zero flux through the glyoxylate shunt, i.e., for the ICL and MALS reactions, which FBA misses qualitatively by setting them to zero. As a consequence, this also leads to a better match of our model with data for reactions ICDHyr and AKGDH that channel pyruvate through the Krebs cycle. Lastly, the maximum entropy model also improved on the FBA predictions for several mutant strains of E. coli that were deleted for various metabolic genes, as shown in Fig 2D.

II.2 Fluxes, variances, and correlations along the optimization trajectory

It is instructive to examine the evolution of the joint distribution over fluxes, $P_{\beta}(\mathbf{v})$, as a function of the optimization parameter, $\beta$. Fig 3A shows how the growth rate approaches the maximal rate achievable, $\lambda_{\rm max}$, with the inferred values of $\beta^{*}$ from Fig 2A,C suggesting an “optimization level” in the range of $\sim 60-80\%$ of the maximum. These levels are reached by adjusting flux values away from what they would have been under uniform sampling from the polytope of the allowed metabolic phenotypes, $\mathcal{P}$. Fig 3B traces the relative changes in all fluxes as a function of $\beta$. Interestingly, in the FBA limit almost half of the fluxes (38 out of 86 fluxes, the upper half of the plot) are forced to zero, whereas at values that fit data best ($\beta^{*}\lambda_{\rm max}\sim 10^{2}$) these fluxes only decrease by about $1/3$ relative to their average value in the uniform sampling limit. Furthermore, the glyoxylate shunt remains active, in agreement with experimental observations. Surprisingly, only for a few reactions the fluxes are predicted to increase with growth rate optimization relative to the uniform sampling (lowest $\sim 5$ fluxes in Fig 3B). These are mainly nitrogen and phosphate transport reactions, and to a lesser extent, MDH and PGI reactions; the latter two reactions are classified as reversible, so the predicted increase may have thermodynamic rather than regulatory reasons.

We separately illustrate three flux behaviors in Fig 3C, for isocytrate lyase (ICL), dehydrogenase (ICDH), and glutamate dehydrogenase (GLUDy). ICL and ICDH track the relative channeling of carbon sources in the Krebs cycle vs glyoxylate shunt; ICL flux is switched off in the $\beta\rightarrow\infty$ limit, whereas ICDH flux remains nearly constant with $\beta$. In contrast, GLUDy reaction is reversible, switching sign at intermediate values of $\beta$, while at high $\beta$ the reaction ultimately gets frozen in the backward direction, implying high levels of ammonia in the cell.

We also evaluated the predicted fluctuations in metabolic fluxes from the maximum entropy model, Eq. (3), at $\beta^{*}\lambda_{\rm max}\sim 10^{2}$, and found a clear division between reactions with high and low coefficients of variation ($CV$). Among tightly controlled fluxes were all glycolytic reactions ($CV<0.3$) with the exception of PGI, as well as all transport reactions related to biomass formation (i.e., for glucose, oxygen, ammonia, carbon dioxide, phosphate ions; $CV<0.11$), the first part of the Krebs cycle, and the irreversible reactions of oxidative phosphorylation.

We next wondered how flux variances scale with the optimization level. In the uniform sampling limit ($\beta=0$) the variances should be large, characteristic of the shape and extent of the permitted polytope, $\mathcal{P}$. While in the FBA limit ($\beta\rightarrow\infty$) the flux fluctuations should vanish, we expect a well-defined scaling regime at high $\beta$ where the variances shrink towards the FBA solution in a manner that is independent of the global polytope properties. This regime is indeed reached for all fluxes at $\bar{\lambda}/\lambda_{\rm max}\gtrsim 0.90$ and for some fluxes much earlier, as shown in Fig 3D: flux fluctuations subsequently decrease with $\beta$ as $\sigma_{i}(\beta)/\sigma_{i}(\beta=0)\propto(1-\bar{\lambda}(\beta)/\lambda_{\rm max})$.

What kind of correlation structure between fluxes does the maximum entropy model predict? While the growth rate $\lambda$ is linear in constituent fluxes in Eq. (3), suggesting that the joint distribution could factorize, correlations between fluxes develop because of the stoichiometric constraints that define the polytope $\mathcal{P}$. A subset of fluxes that we focus on in Fig 3 exhibits a clear structure of strong (anti-)correlation both under uniform sampling (Fig 3E) and in the FBA limit (Fig 3F). The FBA pattern of correlations, in particular, can easily be partitioned into four groups using a clustering algorithm Slonim et al. (2005) so that the groups are strongly enriched for reactions characteristic of glycolysis, glyoxylate shunt, pentose phosphate pathway, and citric acid cycle, respectively. Fluxes in the glycolysis cluster tend to correlate strongly with fluxes in the citric acid cycle cluster, but anti-correlate with glyoxylate shunt and pentose phosphate pathway cluster. Comparison of the FBA correlations (F) with the uniform sampling (E) reveals that stoichiometric constraints alone shape much of the correlation structure, with the exception of anti-correlation between glycolysis and glyoxylate shunt clusters, which is a distinct consequence of the growth rate optimization. More generally, it is intriguing to apply maximum entropy to recover the correlation structure of metabolic fluxes in the FBA limit and use that to identify, automatically via clustering, separate metabolic pathways (see SI for correlation between all fluxes).

II.3 Lower limit to regulatory information required for high growth rates

As growth rate optimization parameter $\beta$ is increased, flux variances shrink (Fig 3D), correlations strengthen (Fig 3F), and the distribution over fluxes within the polytope $\mathcal{P}$ localizes closer to the FBA solution, $\mathbf{v}_{\rm max}$. The degree of localization can be quantified by the entropy of the joint distribution over fluxes:

[TABLE]

Because $P_{\beta}(\mathbf{v})$ is, by construction, a maximum entropy distribution with average growth rate $\bar{\lambda}(\beta)$, the decrease in entropy, $I(\beta)=S(\beta=0)-S(\beta)$, is a measure for the minimal amount of information necessary to control the fluxes and achieve a given average growth rate. This is shown graphically in Fig 3G, where we plot the average growth rate, $\bar{\lambda}$, as a function of information, $I$ (expressed in bits), parametrically in $\beta$. The resulting curve divides the $(I,\bar{\lambda})$ plane into two halves: while it is possible to achieve metabolic phenotypes below the $I(\bar{\lambda})$ curve, the dashed region above the curve is forbidden. This is because no distribution exists that achieves high growth rates $\bar{\lambda}$ without also deviating from the uniform distribution by at least the required number of bits.

Fig 3G suggests that at least $\sim 40$ bits of information are required to control the fluxes and reach growth rates amounting to $\sim 80\%$ of the maximal rate, $\lambda_{\rm max}$, as reported in data for E. coli; higher growth rates call for increasing amounts of information, which formally diverges in the FBA-limit as $\beta\rightarrow\infty$.

Cells control metabolic fluxes through regulatory networks, either indirectly, by regulating the expression of metabolic enzymes, or directly, by modulating the enzymatic activity through various feedback loops; either way, metabolic resources are required to exert this control. This leads to a trade-off: flux control is necessary to support a high growth rate, but itself carries a growth rate penalty. We created a simple toy model to capture this intuition (see SI Apprndix). Here, $K$ regulatory pathways control the fluxes and each pathway is modelled as a Gaussian information channel, so that together, these channels provide $I(\beta)$ bits of necessary information as shown in Fig 3G. The signal-to-noise of each regulatory channel is determined by the number of regulatory molecules: higher molecular counts enable precise control and thus higher information, but impose higher cost. In this model, the cost-free growth rate at given $\beta$ is reduced by the cost to support $K$ channels which control the fluxes, so that the resulting effective growth rate is:

[TABLE]

where $\alpha$ determines the metabolic cost of regulatory molecules, and we estimated the number of regulatory pathways, $K$, to be approximately the number of degrees of freedom of the flux polytope, $K\approx D$. The cost of regulation clearly limits the achievable growth rate, as shown in Fig 3H, where the $\bar{\lambda}_{\rm eff}(\beta)$ curves now develop a maximum rather than increasing monotonically as in the cost-free case of Fig 3A. While our toy model is very simplistic, it does capture properly the scaling of information with the growth rate, as well as the exponential metabolic cost of achieving high information transmission in molecular networks, reported previously Tkačik and Walczak (2011); Tkačik and Bialek (2016). Thus, among many possible constraints acting on a cell, the cost of regulating metabolism itself Kacser et al. (1995) can impose non-negligible limits to growth.

II.4 Experimental test of growth rate fluctuation scaling

Can we test the novel predictions of our theory that extend beyond the domain of validity of the flux balance analysis? While it is experimentally unfeasible to measure metabolic fluxes and their fluctuations at the single cell level, one can tractably measure division times and growth rates for single Escherichia coli cells growing in stable conditions for long periods of time. In our model, such growth measurements directly connect to the biomass producing reaction with its associated “growth flux,” $\lambda(\mathbf{v})$. Fig 3D suggests that flux fluctuations should scale $\propto(\lambda_{\rm max}-\bar{\lambda})$, and since the growth flux is a linear combination of bona fide metabolic fluxes, its fluctuations, too, should follow the same scaling. To verify this explicitly, we computed the fluctuations in growth rate, $\sigma/\lambda_{\rm max}$, as a function of the optimization parameter $\beta$, in Fig 4A. In the range of $\beta\lambda_{\rm max}\gtrsim 40$, characteristic of both wild type E. coli experiments and many mutants, the predicted growth fluctuations indeed obey

[TABLE]

we refer to this range as the “scaling regime.” Beyond variance, the complete distribution of growth rates, $Q(\lambda)$, can be sampled by marginalizing the maximum entropy model, Eq. (3).

Measurements of single-cell growth rates allow us to estimate growth rate distributions and compare them to the predicted $Q(\lambda)$, as well as to empirically extracted $\bar{\lambda}$, $\lambda_{\rm max}$, and the fluctuations $\sigma$, to verify the predicted relation of Eq. (7). We used previously published data Bergmiller et al. (2017) where E. coli cells were stably grown in the “mother machine” microfluidic device while multiple sub-inhibitory steps of concentration of the antibiotic tetracycline were delivered as shown in Fig 4B. Low concentrations of antibiotic allowed us to probe different average growth rates in the same setup, and to construct empirical distributions of growth rates for every antibiotic concentration by pooling data from technical replicates of the multi-step experiments (SI Appendix). We find an excellent match between measured and predicted growth rate distributions in Fig 4C for all five concentrations of the antibiotic used. Looking at many individual lineages in separate microfluidic channels, we can also extract $\lambda_{\rm max}$, $\bar{\lambda}$, and $\sigma$ per lineage empirically and confirm the predicted scaling of growth rate fluctuations, as shown in Fig 4D.

II.5 Fluctuation-response relationship and associated timescales

Our approach so far has been limited to modeling steady state growth. Under mild conditions, we can also study the dynamical response of the network in linear regime under small perturbations. A simple biologically-motivated dynamics is provided by the minimal model of diffusion-replication inside the metabolic space De Martino et al. (2016), described by the one-parameter ($D$, the diffusion coefficient) equation for the growth rate distribution $p(\lambda)$:

[TABLE]

where $q(\lambda)$ is the marginal distribution of the growth rate at uniform sampling ($\beta=0$). An analytical solution of this equation can be obtained in the limit of small $D$, leading to the following scaling laws for the typical response times as well as the growth rate fluctuation autocorrelation time $\tau$ De Martino and Masoero (2016):

[TABLE]

These experimentally testable relations predict a divergent slowing down of the response time with growth rate maximization. As a consequence, growth rate fluctuations could take on a functional role in speeding up the response to environmental perturbations, e.g., to nutritional up-shifts or externally applied stresses. Even if the experimental test of such dynamic predictions is beyond the scope of this paper, our model connects to a wide range of currently ongoing metabolism- and growth-related investigations.

Discussion

In this work we considered maximum entropy distributions at fixed average growth rate in the space of metabolic phenotypes, a straightforward and statistically rigorous extension of the flux balance analysis, which is recovered in the asymptotic limit. Experimental estimates of enzymatic fluxes of the central carbon core metabolism in bulk cultures of E. coli, as well as empirical growth rate distributions of E. coli collected from single cell measurements, are consistent with intermediate level of growth optimization ($\beta\lambda_{\rm max}\sim 10^{2}$ and $\bar{\lambda}/\lambda_{\rm max}\sim 0.6-0.8$). We find that fluctuations can be captured by a simple maximum entropy model, and that the zero-fluctuation FBA limit qualitatively misses important experimental facts, e.g., the observed non-zero fluxes through the glyoxylate shunt.

The improved ability of our model to match the flux measurements is a consequence of a single extra parameter, $\beta$, which can easily be determined from existing experimental data. Beyond a better fit, however, our model also makes a wide range of predictions, extending the domain of metabolic network analysis to the single-cell level. While it is difficult to measure the single-cell metabolic fluxes and their fluctuations in isogenic populations in steady state, such measurements for the growth rate are increasingly available. This connection enables the new predictions of our theory to be tested, and opens up the theory for verifiable extensions. Validating the predicted scaling of growth rate fluctuations in Fig 4 is only the first step, with two broad lines of investigation within reach.

First, our approach is not limited to the core catabolism analyzed here or to bacterial metabolism, but can in principle be extended to other genome-scale networks. In practice, however, we often lack suitable large-scale experimental flux measurements. It is also likely that physico-chemical constraints alone are insufficient to yield quantitatively accurate predictions. Similar issues arise also in the core catabolism for high growth rates that exceed the threshold of the acetate switch, for which additional constraints have to be added in FBA-based approaches Mori et al. (2016). Our method can be extended to accommodate such cases, or systems where strict growth maximization is likely not a suitable objective. Extra objectives or constraints in the maximum entropy would appear as additional terms in the exponent of Eq. (3), where their corresponding parameters would control various trade-offs between the objectives. This flexibility may be required to model metabolic dependencies, cell type heterogeneity, or interactions between cells.

Second, we sketched how the maximum entropy model could be extended to dynamics through fluctuation-response relations. This requires further assumptions that need to be tested separately, but makes a very strong prediction about the link between the autocorrelation time of growth fluctuations and the typical response time to, e.g., nutrient shifts. This link appears fundamental, since the response time is a central biological quantity measurable in bulk, while the fluctuation autocorrelations are microscopic, single-cell properties, which can be measured with recent experimental setups. Interestingly, the predicted response times lengthen with the degree of growth rate optimization, suggesting a trade-off between responsiveness to changes and efficiency in steady state; as a consequence, it is unclear whether the evolutionarily optimal outcome should be equated to complete growth rate optimization with no fluctuations, e.g., the FBA limit. Quantitatively, in stable environments where E. coli grows well and possibly achieves a high degree of growth rate optimization, one could experimentally look for signatures of long-timescale fluctuations, either directly in the growth signal, or by proxy through constitutive gene expression. Curiously, we report that the parameter $\beta$ of our model has the dimension of time, whose best-fit value inferred from E. coli data is of the order of 1 day. A proper investigation into the role of $\beta$ deserves further theoretical work.

Beyond extensions to dynamics, our analysis made two further theoretical contributions. First, it clarified the relative roles of stoichiometric constraints and the growth optimization assumption in flux balance analysis. The maximum entropy model is an explicit construction of a smooth interpolation between the uniform regime (where only stoichiometric constraints are active) and the FBA (where growth is maximized in addition). The uniform limit is a natural baseline—where no control is exerted by the cell—against which to compare the observed fluxes, their fluctuations, and correlations, as we have done in Fig 3. Without this baseline comparison, it is hard to assess how surprising the observations of metabolic optimality should really be Ibarra et al. (2002). Our second theoretical contribution is the observation that a certain minimal information is needed to achieve a desired growth rate (Fig 3G, H). This information is expressed in the same currency (bits) in which we measure the performance of regulatory networks, enabling us to suggest a tradeoff that sets the optimal degree of metabolic control. Contrary to other cellular networks where estimation of information only has been done for single network components or simple pathways Tkačik and Bialek (2016), the metabolic network is the sole case where we could estimate the lower bound on the required number of regulatory bits. Our statistical mechanics approach thus opens a connection between metabolic networks and their regulatory counterparts, which is both of theoretical interest and could also be probed in comparative genomic studies.

Acknowledgements.

We acknowledge the support of the Austrian Science Fund grant FWF P28844 (G.T.) and of the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement no. $[291734]$ (D.D.M)

Appendix

The sampling Monte Carlo Algorithm

We sampled the space of steady states of the model using a hit-and-run Markov Chain Monte Carlo algorithm, subject to an ellipsoidal preprocessing step in order to tackle ill-conditioning. The uniform sampling of convex bodies by means of Monte Carlo methods has been considered a breakthrough in computational convex analysis, e.g., it makes feasible the calculation of the volume which is otherwise a computationally difficult problem Dyer and Frieze (1988). A fast and popular algorithm to sample points inside convex bodies is the hit-and-run (HR) Markov Chain Monte Carlo Turcin (1971); Smith (1996), which works along the following lines. Given a $D$-dimensional convex polytope $\mathcal{P}$, from which one wants to sample, and a point inside the polytope, $x_{k}\in\mathcal{P}$:

• •

Choose a uniformly distributed direction $\theta_{k}$, that is, a point generated from the uniform distribution on the $D$-dimensional unit sphere. This can be done with the Marsaglia method, i.e., by generating $D$ independent gaussian random variables with zero mean and unit variance, and then normalizing the vector to unit length.

• •

Generate $t$ uniformly on the interval $[t_{\rm min},t_{\rm max}]$, where $t_{\rm min}$ ($t_{\rm max}$) is the minimum (maximum) value of $t$ such that $x_{k}+t\theta_{k}\in\mathcal{P}$;

• •

Update $x_{k+1}=x_{k}+t\theta_{k}$, start again.

The starting point can be found, for instance, by interpolating between two vertices obtained by linear programming. The second step requires finding the intersections among a line and $\mathcal{P}$. In order to perform the HR dynamics we should always use a full-dimensional representation of the convex set. The mixing time of the HR, i.e., the time to converge to the desired distribution, scales as a polynomial of the dimensions of the body but the method can suffer from ill-conditioning; more precisely, the mixing time $\tau$ scales like Lovász (1999)

[TABLE]

where $R,r$ are the radii of respectively the minimum inscribing and the maximum inscribed balls. The prefactor $R/r$ can be reduced to a polynomial of $D$ by extracting $\theta_{k}$ from the surface of a matching ellipsoid instead of the unit sphere. The ellipsoid of maximum volume inscribed inside $\mathcal{P}$, i.e., the Loewner-John ellipsoid Ball (1997), would reduce the prefactor to $O(D)$, but the problem of finding it is NP-hard. Below we describe a method due to L. Lovasz Lovász (1987) that finds in polynomial time an approximate matching ellipsoid that reduces the prefactor to $O(D^{3/2})$.

Preprocessing: the Lovasz algorithm

We want to construct a couple of concentric ellipsoids $E,E^{\prime}$ matching the polytope $\mathcal{P}$, i.e., such that $E^{\prime}\subseteq\mathcal{P}\subseteq E$, where $E^{\prime}$ is obtained from E shrinking by a factor $O(D^{3/2})$. This is called weak Loewner–John pair. We define a series of enclosing ellipsoids $E_{k}$, starting with $E$ as the sphere with center in the origin and radius $R$ large enough in order to inscribe the body, according to the following lines:

• •

INPUT: An ellipsoid $E_{k}$ with its center $x_{k}$.

• •

Check if $x_{k}\in\mathcal{P}$, if yes go to 2, if no go to 1.

• •

1. Consider a hyperplane separating $x_{k}$ and $\mathcal{P}$, and the halfspace H enclosing $\mathcal{P}$, calculate the ellipsoid of minimal volume enclosing $H\bigcap E_{k}$; go to OUTPUT 1.

• •

2. Determine the endpoints of the axis of $E_{k}$, shrink the ellipsoid and check if the shrunk ellipsoid $E^{\prime}_{k}$ is inside the body. If yes, go to OUTPUT 2, if no go to 3.

• •

3. Consider the endpoint of an axis of the shrunk ellipsoid outside $\mathcal{P}$, e.g. $x^{\prime}_{k}$; consider a hyperplane separating $x^{\prime}_{k}$ and $\mathcal{P}$, and the halfspace enclosing $\mathcal{P}$; calculate the ellipsoid of minimal volume enclosing $H\bigcap E^{\prime}_{k}$; go to OUTPUT 1.

• •

OUTPUT 1: A new ellipsoid $E_{k+1}$ of lower volume with center $x_{k+1}$, update $k$, repeat from INPUT.

• •

OUTPUT 2: A weak Loewner-John ellipsoid.

Upon calculating the reduction in volume of the enclosing ellipsoid after one step, it can be demonstrated that this series converges in polynomial time to a weak Loewner–John pair. We refer to Bland et al. (1981) for demonstrations and formulae.

Source Code

We have provided in doi:10.15479/AT:ISTA:62 a C++ code implementing the Lovasz preprocessing as well as the Hit-and-Run algorithm and the polytope representation of the metabolic network used in this study. Please refer to the README.txt file for further information.

Metabolic network and flux data

We considered the model of the E. coli catabolic core from the metabolic genome scale reconstruction iAF1260 Orth et al. (2011). Apart from mass-balance (purely stoichiometric) constraints, we have considered the bounds on reaction fluxes from the reversibility assignments provided within the model and default bounds on uptakes reflecting glucose limited aerobic conditions. We considered the dataset from Ishii et al. (2007) that includes extensive flux, enzyme, mRNA and metabolite level measurements of E. coli populations growing in steady-state in a glucose limited minimal medium in several conditions (dilution rates and/or gene knock outs). We considered further experimental estimates of the core carbon metabolism from the database Zhang et al. (2014b), where we collected experiments Farmer and Liao (1997); Siddiquee et al. (2004); Schaub et al. (2008); Peng et al. (2004); Emmerling et al. (2002); Hua et al. (2003); Nanchen et al. (2008); Nikel et al. (2009); Li et al. (2006); Toya et al. (2007); Zhu et al. (2003); Zamboni et al. (2009) on wild type E. coli performed in glucose limited medium and dilution rates not exceeding $0.4$h*-1*, i.e. below the acetate switch Basan et al. (2015)

Data analysis details

Figure 2.

Flux intensities have been normalized with respect to the glucose uptake (i.e. expressed in relative values) in order to compare datasets with different growth rate/dilution (in all cases case below $0.4$h*-1*). While the Lagrange parameter $\beta$ can be used to fix the average growth rate, in this case we consider the best fit for $\beta$ upon comparing relative flux values, i.e., we look for the minimum of $\chi^{2}=N_{f}^{-1}\sum_{i=1}^{N_{f}}\left(\langle v_{i}\rangle-V_{i}\right)^{2}/\left(\sigma^{2}_{i}+E^{2}_{i}\right)$, where $V_{i}$ is the measured flux and $E_{i}^{2}$ the associated measurement error variance, while $\langle v_{i}\rangle$ and $\sigma_{i}^{2}$ are the mean (and variance, respectively) of the corresponding flux computed in the maximum entropy model.

• •

*Figure 2B

*We considered $10^{2}$ values of $\beta$ spanning the interval $\left[10^{-2},10^{4}\right]$ uniformly in logarithmic scale. For each value of $\beta$ we sampled $10^{5}$ flux configurations according to the max entropy distribution with the previously described hit-and-run Monte Carlo algorithm from which we calculate the estimates of the averages and variances. These are then compared to experimental values through the $\chi^{2}$, as well as the mean square error, $MSE=N_{f}^{-1}\sum_{i=1}^{N_{f}}\left(\langle v_{i}\rangle-V_{i}\right)^{2}$.

• •

*Figure 2A,C

*The flux values reported for the max entropy model refer to the value of $\beta$ that minimizes the $\chi^{2}$.

• •

Figure 2D

We simulated an enzyme knock out by removing the corresponding reaction from the model (one equality constraint). This leads to a new metabolic space that has been analyzed along the same lines as the unperturbed one. Flux estimates for knockouts from Ishii et al. (2007) consist of one single experiment with no repeats. The MSE improvement is the fractional difference between the minimum value and the value retrieved by FBA.

Figure 3.

In these plots, flux intensities have been normalized with respect to the glucose uptake, while the growth rate has been rescaled by its maximum.

• •

*Figure 3A,G

*The curves were calculated analytically along the following lines. We fit the marginal probability density of the growth rate, otained numerically by Monte Carlo sampling, with a Beta distribution, i.e. in relative units $x=\lambda/\lambda_{\rm max}$:

[TABLE]

finding simply $a=D-1=22$, where $D$ is the dimension of the polytope, i.e. the growth rate is maximized in a subspace of dimension [math] (a vertex) De Martino and Masoero (2016). The average growth rate as a function of $\beta$ can be calculated from the normalizing factor of the max entropy distribution,

[TABLE]

where $\gamma$ is the lower incomplete gamma function and finally (Figure 3A)

[TABLE]

The minimal entropy reduction of the metabolic space upon fixing the average growth rate (Figure 3G) has been obtained by inverting the previous equation, $\beta=\beta(\overline{x})$, and performing a Legendre transform De Martino et al. (2016)

[TABLE]

Both curves can be calculated numerically, the first directly, upon sampling the average growth rate at fixed $\beta$, while the second can be obtained in parametric form upon integration; in both cases we have found an excellent agreement with the analytical approximation.

• •

*Figure 3B-F

*Averages, variances and correlations among fluxes have been calculated numerically by sampling max entropy distribution of the metabolic space with the previously described Monte Carlo method. We considered $10^{2}$ values of $\beta$ spanning the interval $\left[10^{-2},10^{4}\right]$ uniformly in logarithmic scale. For each value of $\beta$ we sampled $10^{5}$ flux configurations. Figure 3F has been obtained by numerical asymptotic analysis of flux correlations in the limit $\beta\to\infty$.

• •

Figure 3H

We considered a simple model where $K$ regulatory channels together provide the information, $I(\bar{\lambda})$, required to support a given average growth rate, $\bar{\lambda}$. Each channel is modeled as an additive Gaussian information channel whose capacity is given by $I_{1}=\frac{1}{2}\log_{2}(1+SNR)$, where $SNR$ is the ratio of the signal to noise variance. To provide sufficient regulatory information we must have $I(\bar{\lambda})\approx KI_{1}$. We assume that the variance of the signal is proportional to its mean, $\bar{S}$, so that $SNR=\alpha_{1}\bar{S}$, for some constant $\alpha_{1}$; this is approximately true in the simplest models of biochemical reaction pathways where the dominant source of noise is the shot (Poisson) noise due to finite number of signaling molecules. We further assume that the metabolic cost of regulation per channel is also proportional to the number of signaling molecules used, i.e., to the mean signal in the channel, and thus the total cost (written as the decrease in the growth rate) is $\Delta\lambda=-\alpha_{2}K\bar{S}$. We can eliminate $\bar{S}$ by writing it out in terms of $SNR$, which can in turn be expressed in terms of $I(\bar{\lambda})$, to yield an expression for the effective growth rate that takes into account the cost of regulation:

[TABLE]

where $\alpha=\alpha_{2}/\alpha_{1}$ is a constant proportional to the metabolic cost of signaling molecules.

Figure 4.

The growth rate data have been obtained by analyzing images from a microfluidics device (mother machine) setting with $30$ channels each containing $3$ E. coli cells growing in steady-state in glucose-limited conditions. Their growth is monitored under increasing exposure to antibiotics (five concentration steps of tetracycline) specific for ribosome translation inhibition for $5$h. Experiments have been performed on $3$ strains; $3$ technical repeats were carried out. Strains differed only by the fluorescent tags whose effects are not expected to have major contributions to the growth rates.

• •

Figure 4A

The curve has been calculated analytically as described for Figure 3.

• •

*Figure 4C

*The distributions have been obtained as follows. We disregard measurements during the first $15$ time frames from each antibiotic step in order to obtain approximately stationary distributions; we also rejected outliers exceeding $4\sigma$. For each antibiotic step, we lump together growth rate measurements pertaining to different channels, time, cell order in the channel (mother-daugthers), strains and repeats; this results in a sample size of $\simeq 1-5\cdot 10^{3}$ per antibiotic step. We performed a best fit ($\chi^{2}$ minimization) of max entropy marginal growth rate distributions in the analytical approximation outlined in discussion for Figure 3 (above), and we retrieved the values for the two parameters $\beta$ and $\lambda_{\rm max}$ that we report in Table 1, alongside of the average and empirical maximum of the growth rate.

• •

Figure 4D:

Each data point stands for a distribution obtained upon lumping growth rate data measured in different channels ($30$) and time points ($\simeq 3$), i.e. it constitutes a unique sample whose size is $O(10^{2})$. Each data point refers thus to a different antibiotic step ($5$), to a different cell ordering in the channel ($3$), strain ($3$) and repeat ($3$) for a total of $135$ points. We considered the usual estimators for the average and variance, while we have reported in all cases the empirical maximum for the growth rate.