A Van-Der-Waals picture for metabolic networks from MaxEnt modeling: inherent bistability and elusive coexistence
AA Van-Der-Waals picture for metabolic networksfrom MaxEnt modeling: inherent bistability &elusive coexistence.
D.De Martino Institute of Science and Technology Austria (IST Austria), Am Campus 1,Klosterneuburg A-3400, Austria
Abstract.
In this work maximum entropy distributions in the space of steady statesof metabolic networks are defined upon constraining the first and second momentof the growth rate. Inherent bistability of fast and slow phenotypes, akin to aVan-Der Waals picture, emerges upon considering control on the average growth(optimization/repression) and its fluctuations (heterogeneity). This is applied to thecarbon catabolic core of
E coli where it agrees with some stylized facts on the persistersphenotype and it provides a quantitative map with metabolic fluxes, opening for thepossibility to detect coexistence from flux data. Preliminary analysis on data for
E Coli cultures in standard conditions shows, on the other hand, degeneracy for the inferredparameters that extend in the coexistence region.
1. Introduction
Quantitative modeling cell metabolism would give essential insights on reconstructingcell physiology from molecular mechanisms with wide applications, a quest that shall beaddressed within a systemic perspective [1]. Upon modeling cell growth, an approachbased on stochiometric constraints under a steady-state hypothesis for the underlyingchemical network, complemented by optimality principles, it works in semi-quantitativeway at population level in lab conditions [2], it is straightforward to formulate [3] andcomputationally fast to implement [4]. However, several phenomena in biology wouldrequire to assess single-cell heterogeneity [5] and thus an extension of modeling beyondthe population level. Nowadays, recent advancements in microfluidic device can giveaccess to growth rate physiology at single cell level [6]. In particular measures onisogenic populations in identical conditions give rise to growth rate distributions withunavoidable fluctuations [7], that shall reflect metabolic variability [8]. Both the issueof a systemic perspective and the one of addressing heterogeneity could get interestinginsights from statistical mechanics methods. Recent applications of Monte Carlo Markovchain algorithms to integral calculus [9] give rise to the possibility of mapping growthrate distributions into an underlying distribution for the metabolic phenotypes, in themost simple way upon recurring to maximum entropy distributions at fixed average a r X i v : . [ q - b i o . M N ] J u l Van-Der-Waals picture for metabolic networks
E Coli endows inherentbistability of the growth rate [16] and it has been observed that the central carboncore metabolism shall thus endow a bistable phenotypic landscape where slow and fastgrowing cells sit on wells separated by a watersheed barrier whose height and positionsare ruled by two phenomenological parameters related to the growth rate and themetabolic flux[17]. In this work a simple phenomenological approach is proposed tomodel quantitatively bistability and coexistence of fast and slow growing phenotypesfor the central catabolic core of
E. Coli in terms of maximum entropy distributions. Itwill be shown that bimodal distributions are retrievec upon constraining the first twomoments of the growth rate in the metabolic space, with the two lagrange multipliersin the Boltzmann-Gibbs distribution that play the role of the control parameters thatshape the bistable phenotypic landscape akin to a Van-Der-Waals picture. It will beshown that, despite its simplicity, this framework captures stylized facts of the metabolicpersistent phenotype, like an active metabolism and increased rate of respiration andcarbon dioxide production (per unit of carbon intake) with respect to the fast growingphenotype. This framework would ideally provide a mean to infer elusive coexistence ofpersisters from flux data and a preliminary analysis from
E. Coli cultures in standardconditions will be presented, alongside with conclusions.
2. Results
A standard approach to model metabolism during cell growth is to consider theunderlying chemical network in a well-mixed steady state in the continuum limit, thatis approximately valid for timescales slower than diffusion and metabolic turnover. Thereaction network is encoded in a matrix of stoichiometric coefficients S iµ (metabolite µ in reaction i ), that linearly relates the time derivative of concentrations c µ to enzymatic Van-Der-Waals picture for metabolic networks ν i (mass balance)˙ c µ = (cid:88) i S iµ ν i (1)Upon considering steady states and taking into account reversibility, kinetic limits andmedium capacity in terms of flux bounds, a space of feasible states is defined that isgeometrically a convex polytope P (cid:88) i S iµ ν i = 0 (2) ν i ∈ [ ν mini , ν maxi ] . (3)In order to model cell growth, a phenomenological biomass production reaction isadded in models, whose flux will be denoted λ and whose maximum in given conditionscan be evaluated with linear programming.Here we will consider Maximum Entropy distributions in the space of steady stateof a metabolic network upon fixing the first two moments of the growth rate, (cid:104) λ (cid:105) and (cid:104) λ (cid:105) . By a standard variational approach [18], constraints on the moments are enforcedby Lagrange multipliers ( β, γ ) as control parameters in Boltzmann-Gibbs distributions P ( ν ) ∝ e βλ ( ν )+ γλ ( ν ) if ν ∈ P (4)whose analysis is developed in the next section. The marginal distribution q ( λ ) of the growth rate from an uniform sampling of thenetwork under exam, the catabolic core of E.coli, has in a good approximation thesimplex-like form q ( λ ) ∝ ( λ max − λ ) a (5)where λ max is the maximum achievable growth rate, given the conditions (obtainablewith linear programming) and a = D − d −
1, where D is the dimension of the spaceand d the dimension of the subspace where λ (eg d = 0 is a vertex). Upon consideringthe aforementioned constraints on the moments, we have for the marginal distributionof the growth rate p ( λ ) ∝ e F β,γ ( λ ) , where the rate function (from now on we assume λ max = 1 for simplicity of notations) has the form F β,γ ( λ ) = βλ + γλ + a log(1 − λ ) λ ∈ [0 ,
1] (6)The extremal points of this function can be analyzed as a function of the controlparameters in simple terms. There is a maximum in λ = 0 (where F (0) = 0) iff F (cid:48) (0) = β − a <
0. There are other two extremal points where F (cid:48) ( λ ± ) = 0, i.e. λ ± = 2 γ − β ± (cid:112) (2 γ + β ) − γa γ (7) Van-Der-Waals picture for metabolic networks β > √ γa − γ , and positive if γ > γ c = a/
2, with λ + maximumand λ − minimum. We have qualitatively a picture similar to a Van-Der-Waals fluid,with the two “spinodal” lines ( γ > γ c = a/ β = a (8) β = (cid:112) γa − γ (9)that mark the onset of probability maxima that refer to respectively, slow (below (5)) andfast (above (6)) growing phenotypes. The curves (5,6) thus define a coexistence regionwith an equilibrium line given implicitly by the equation F (0) = 0 = F ( λ + ), wherethe maxima have equal heights, that marks discontinuous transitions and hysteresis,shrinking at the point ( γ c = a/ , β c = a ). This picture is summarized in the diagram inFig. 1. Growth rate marginal distributions in particular can be obtained analytically γ -100102030405060 β equilibrium fast growthslow growth coexistence Figure 1.
Phase diagram of metabolic network states in the plane ( γ, β ) of the twolagrange multipliers constraining the first two moments of the growth rate ( a = 20).Red dashed lines: “spinodal” lines that mark the onset of probability peaks for slowand fast growing phenotypes. Blak full line: “equilibrium” curve where the frequenciesof fast and slow phenotypes are equal. by Legendre transform as we show in fig 2, where a comparison with Monte Carlocomputations on the network under exam (see materials and methods section ) is shown.More in general, to each point of the phase diagram it corresponds a given distributionof metabolic fluxes, and in next section a simple analysis will be performed for a pointinside the coexistence region for the model under exam, that is the central carbonmetabolism of E Coli in glucose limited aerobic conditions and low dilution. Van-Der-Waals picture for metabolic networks β = γ =0 pd f λ / λ max montecarloanalytics β =50 γ =0 pd f λ / λ max montecarloanalytics 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 0.2 0.4 0.6 0.8 1 β =0 γ =50 pd f λ / λ max montecarloanalytics Figure 2.
Growth rate distributions for slow growing ( β = γ = 0), fast growing( β = 50 , γ = 0) and coexisting states ( β = 0 , γ = 50 ) from analytical calculations(Legendre transform) and Monte Carlo computations on the network under exam (seematerials and methods) We consider the point ( β = 0 , γ = 50 ) in the coexistence region of the phase diagram.In Fig 3 we show scatter plots of a) the rate of oxygen consumption, b) carbondioxide production and c) ATP synthase flux, all relative to the glucose uptake andscattered with the respect to the growth rate (relative to the maximum achievableone). In qualitative agreement with stylizied facts, the slow growing phenotype resultsmetabolically active as indicated by comparable level of ATP synthase flux with respectto the fast growing phenotype (and increased variance). Further, the level of oxygenconsumption and carbon dioxide production of the slow growing phenotype is higherwith respect to the fast growing phenotype, clearly indicating this as a diversion of thecarbon source intake from the biomass. However the absolute rate of glucose uptake offast and slow growing phenotypes are comparable within this simple framework, at oddswith experiments with induced persistence where slower cells intake less carbon units. Van-Der-Waals picture for metabolic networks
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 λ / λ max -6-5-4-3-2-1 0 C O fl u x ( r e l a t i v e )
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 λ / λ max O fl u x ( r e l a t i v e )
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 λ / λ max A T P s y n t h a s e fl u x ( r e l a t i v e ) d e n s i t y ( a . u . ) Figure 3.
Scatter plots of a) the rate of oxygen consumption, b) carbon dioxideproduction and c) ATP synthase flux, relative to the glucose uptake and scatteredwith the respect to the relative growth rate, for a point ( γ = 50 , β = 0) in the phasediagram in fig. 1. From numerical Monte Carlo calculations on the catabolic core of E.Coli . The flux distributions obtained at varying ( β, γ ) can be then compared toexperimental data in order to infer coexistence of fast and slow growing phenotypesas we stress in the following section.
Here it is examined the problem of inferring the parameters ( β ∗ , γ ∗ ) that maximize thelikelihood of data. For simplicty we consider a Gaussian approximation , i.e. if we havea subset of K experimental fluxes, averages and variances, ν i,exp , σ i,exp we look at the χ as a function of the parameters ( β, γ ) that is χ ( β, γ ) = K (cid:88) i ( ν i,exp − (cid:104) ν i (cid:105) β,γ ) σ i,β,γ + σ i,exp (10)A preliminary analysis has been performed to data for E.Coli cultures in standardconditions (see materials and methods section). In Fig. 4 it is shown the χ ( β, γ ) heat-map comparing data and predicted flux distributions in the plane ( β, γ ). The inferencecorrectly predict growing phenotypes, but with a certain level of degeneracy for the Van-Der-Waals picture for metabolic networks γ β χ coexistenceslow growthfast growth Figure 4. χ ( γ, β ) heat-map comparing predicted flux distributions in the plane( γ, β ) and data referring to E Coli core metabolism for cultures grown in glucoselimited aerobic conditions with growth/dilution rates below 0 . − , averages over 35technical replicates, control experiments retrieved from the database [19]. Dashed line:80 = β + 1 . γ . inferred parameters, with low values of the χ that extend in the coexistence region.More precisely a line of χ minima develops approximately given by 80 = β + 1 . γ .
3. Materials and methods
The network analyzed is the catabolic core of a network reconstruction of
E. Coli metabolism [20], that includes glycolysis, pentose phosphate pathway, Krebs cycle,oxidative phosphorylation and several nitrogen catabolic reactions. The underlyingpolytope of steady states in glucose limited aerobic conditions turns out to be a D = 23 − subdimensional variety in a space with N = 86 (number of reactions)dimensions. States can be efficiently sampled with an hit-and-run Monte Carlo Markovchain upon handling of ill conditioning due to heterogeneous scales [21], with roundingellipsoids, for which a method due to L.Lovazs has been implemented [22]. The bias dueto the Boltzmann-Gibbs factor enforcing the moments constraint has been implementedwith a Metropolis-Hastings rejection rule for the proposed step during the hit-and-run dynamics [23]. The flux data refer to E. Coli cultures grown in glucose limitedminimal medium in aerobic conditions at low dilution/growth rates below 0 . − , i.e.the threshold for the acetate switch. The sample consists of 35 technical replicatescollected from control experiments retrieved in the database [19] (same dataset analyzedin [11]). Van-Der-Waals picture for metabolic networks
4. Conclusions
In this work it has been shown that a simple phenomenological approach in which thefirst two moments of the growth rate in the space of steady states of a metabolic networkare constrained can be used to model quantitatively bistability and coexistence of fastand growing cells in heterogeneous cell populations. This shows that bistability andcoexistence can be an inherent feature of the central carbon core upon assuming controlof the level of the average growth rate (optimization/repression) and its fluctuations(heterogeneity), providing in turn a quantitative map with metabolic flux distributions.It has been shown that this simple framework can account for stylized facts on persistersphenotypic heterogeneity, in particular their active metabolism at the level of ATPproduction, oxygen consumption rate and carbon dioxide production. This simpleframework, however fails to predict lower level of carbon uptakes by the slow growingpehnotype as observed experimentally. This last experimental fact refers however tostressed conditions, something that could be added phenomenologically by constrainingthe average covariance of biomass and uptakes, (cid:104) uλ (cid:105) (with a corresponding term inthe Boltzmann-Gibbs weight ∝ e δλu ). Such Maximum Etrnopy approach can befurther used to infer from flux data elusive coexistence and a preliminary analysis ondata about E Coli cultures grown in glucose limited aerobic minimal medium showdegeneracy of parameters that extends in the coexistence region. This would requirefurther investigations upon complementing flux data with measures of single cell growthphysiology, in particular rare events and in microfluidic devices where slower persistersare not taken over. Finally it is worth to mention that in general a Maximum Entropyapproach shows that constraining the first two moments of any flux in the metabolicspace can lead to a Van-Der-Waals picture with coexistence and bistability for theunderlying flux, in turn providing an effective way to model single cell heterogeneity forany secondary metabolites production. Acknowledgments
The research leading to these results has received funding from the People Programme(Marie Curie Actions) of the European Union’s Seventh Framework Programme(
F P / − n [291734]. The Author thanks Prof.M. Heinemann for interesting discussions and kind hospitality in the University ofGroningen. References [1] H Kacser, , and J Burns. The control of flux. In
Symp. Soc. Exp. Biol. , volume 27, pages 65–104,1973.[2] RA Majewski and MM Domach. Simple constrained-optimization view of acetate overflow in e.coli.
Biotechnology and bioengineering , 35(7):732–738, 1990.[3] Leonid Vitalievich Kantorovich. A new method of solving of some classes of extremal problems.In
Dokl. Akad. Nauk SSSR , volume 28, pages 211–214, 1940.
Van-Der-Waals picture for metabolic networks [4] George Dantzig. Linear programming and extensions . Princeton university press, 2016.[5] Steven J Altschuler and Lani F Wu. Cellular heterogeneity: do differences make a difference?
Cell , 141(4):559–563, 2010.[6] Ping Wang, Lydia Robert, James Pelletier, Wei Lien Dang, Francois Taddei, Andrew Wright, andSuckjoon Jun. Robust growth of escherichia coli.
Current biology , 20(12):1099–1103, 2010.[7] Andrew S Kennard, Matteo Osella, Avelino Javer, Jacopo Grilli, Philippe Nghe, Sander J Tans,Pietro Cicuta, and Marco Cosentino Lagomarsino. Individuality and universality in the growth-division laws of single e. coli cells.
Physical Review E , 93(1):012408, 2016.[8] Daniel J Kiviet, Philippe Nghe, Noreen Walker, Sarah Boulineau, Vanda Sunderlikova, andSander J Tans. Stochasticity of metabolism and growth at the single-cell level.
Nature , 2014.[9] V Turcin. On the computation of multidimensional integrals by the Monte Carlo method.
ThProbab Appl , 16:720–724, 1971.[10] Daniele De Martino, Fabrizio Capuani, and Andrea De Martino. Growth against entropy inbacterial metabolism: the phenotypic trade-off behind empirical growth rate distributions in e.coli.
Physical biology , 13:036005, 2016.[11] Daniele De Martino, Anna Andersson, Tobias Bergmiller, C˘alin C Guet, and Gaˇsper Tkaˇcik.Statistical mechanics for metabolic networks during steady-state growth. arXiv preprintarXiv:1703.01818 , 2017.[12] Vanessa M DCosta, Christine E King, Lindsay Kalan, Mariya Morar, Wilson WL Sung, CarstenSchwarz, Duane Froese, Grant Zazula, Fabrice Calmels, Regis Debruyne, et al. Antibioticresistance is ancient.
Nature , 477(7365):457–461, 2011.[13] Nathalie Q Balaban, Jack Merrin, Remy Chait, Lukasz Kowalik, and Stanislas Leibler. Bacterialpersistence as a phenotypic switch.
Science , 305(5690):1622–1625, 2004.[14] Rosalind Allen and Bartlomiej Waclaw. Antibiotic resistance: a physicists view.
Physical biology ,13(4):045001, 2016.[15] Oliver Kotte, Benjamin Volkmer, Jakub L Radzikowski, and Matthias Heinemann. Phenotypicbistability in escherichia coli’s central carbon metabolism.
Molecular systems biology , 10(7):736,2014.[16] J Barrett Deris, Minsu Kim, Zhongge Zhang, Hiroyuki Okano, Rutger Hermsen, AlexanderGroisman, and Terence Hwa. The innate growth bistability and fitness landscapes of antibiotic-resistant bacteria.
Science , 342(6162):1237435, 2013.[17] Jakub Leszek Radzikowski, Silke Vedelaar, David Siegel, ´Alvaro Dario Ortega, Alexander Schmidt,and Matthias Heinemann. Bacterial persistence is an active σ s stress response to metabolic fluxlimitation. Molecular Systems Biology , 12(9):882, 2016.[18] Edwin T Jaynes. Information theory and statistical mechanics.
Physical review , 106(4):620, 1957.[19] Zhengdong Zhang, Tie Shen, Bin Rui, Wenwei Zhou, Xiangfei Zhou, Chuanyu Shang, Chenwei Xin,Xiaoguang Liu, Gang Li, Jiansi Jiang, et al. Cecafdb: a curated database for the documentation,visualization and comparative analysis of central carbon metabolic flux distributions exploredby 13c-fluxomics.
Nucleic acids research , page gku1137, 2014.[20] Jeffrey D Orth, Ronan MT Fleming, and Bernhard Ø Palsson. Reconstruction and use of microbialmetabolic networks: the core escherichia coli metabolic model as an educational guide.
EcoSalplus , 4(1), 2010.[21] Daniele De Martino, Matteo Mori, and Valerio Parisi. Uniform sampling of steady states inmetabolic networks: Heterogeneous scales and rounding.
PLoS ONE , 10(4):e0122670, 2015.[22] L´aszl´o Lov´asz.
An algorithmic theory of numbers, graphs and convexity , chapter 2, page 59. SIAM,1984.[23] W Keith Hastings. Monte carlo sampling methods using markov chains and their applications.