Maximum entropy and population heterogeneity in continuous cell cultures
aa r X i v : . [ q - b i o . M N ] D ec Maximum entropy and populationheterogeneity in continuous cell cultures
Jorge Fernandez-de-Cossio-Diaz ∗ and Roberto Mulet † Group of Complex Systems and Statistical Physics, Department ofTheoretical Physics, University of Havana, Physics Faculty, Cuba Systems Biology Department, Center of Molecular Immunology,Havana, Cuba Group of Statistical Inference and Computational Biology, ItalianInstitute for Genomic Medicine, ItalyDecember 19, 2018
Abstract
Continuous cultures of mammalian cells are complex systems displayinghallmark phenomena of nonlinear dynamics, such as multi-stability, hystere-sis, as well as sharp transitions between different metabolic states. In thiscontext mathematical models may suggest control strategies to steer the sys-tem towards desired states. Although even clonal populations are known toexhibit cell-to-cell variability, most of the currently studied models assumethat the population is homogeneous. To overcome this limitation, we use themaximum entropy principle to model the phenotypic distribution of cells ina chemostat as a function of the dilution rate. We consider the coupling be-tween cell metabolism and extracellular variables describing the state of thebioreactor and take into account the impact of toxic byproduct accumulationon cell viability. We present a formal solution for the stationary state of the ∗ [email protected] † mulet@fisica.uh.cu hemostat and show how to apply it in two examples. First, a simplifiedmodel of cell metabolism where the exact solution is tractable, and then agenome-scale metabolic network of the Chinese hamster ovary (CHO) cellline. Along the way we discuss several consequences of heterogeneity, suchas: qualitative changes in the dynamical landscape of the system, increas-ing concentrations of byproducts that vanish in the homogeneous case, andlarger population sizes. uthor Summary While the advantages of continuous culture in the biotechnological industry havebeen widely advocated in the literature, its adoption over batch or fed-batch modesstalls due to the complexities of these systems. In particular, continuous cell cul-tures display hallmark nonlinear phenomena such as multi-stability, hysteresis,and sharp transitions between metabolic phenotypes. Moreover, the impact of theheterogeneity of a cell population on these features is not well understood. We usethe maximum entropy principle to model the phenotypic distribution of an hetero-geneous population of cells in a chemostat. Given the metabolic network and thedilution rate, we obtain a self-consistent solution for the stationary distribution ofmetabolic fluxes in cells. We apply the formalism in two examples: a simplifiedmodel where the exact solution is tractable, and a genome-scale metabolic net-work of the Chinese hamster ovary (CHO) cell line widely used in industry. Wedemonstrate that heterogeneity may be responsible for qualitative changes in thedynamical landscape of the system, like the disappeareance of a bistable regime,the increase of concentrations of byproducts that vanish in the homogeneous sys-tem and larger number of cells. We explain the causes behind these phenomena.3 ntroduction
Recombinant protein production requires suitable cell hosts and culture conditions[1]. For this purpose mammalian cells are often grown in chemostat-like cultureswhere a continuous flow of incoming fresh media replaces culture liquid contain-ing cells and metabolites. Alternative processes such as batch or fed-batch are alsoadopted by many industrial facilities, but the advantages of the continuous modehave been predicted to drive its wide adoption in the near future [2, 3, 4, 5, 6, 7].However, experiments have demonstrated that continuous cultures exhibit hall-mark phenomena of nonlinear dynamics, such as multiple steady states underidentical external conditions [8, 9, 10, 11] and hysteresis loops [12, 8, 13]. So-phisticated control strategies are then required to drive the system towards desiredsteady states.In this context, mathematical modeling has been used with some success [13,14, 15], but most of these works deal with simple cell populations, consisting atmost of a few competing species [15, 16]. Although it is known that no two cellsin culture are alike [17], the effects of individual cell-to-cell variability are seldomconsidered [18]. Attempts to model heterogeneity in cell cultures are often basedin population balance models [19] or similar approaches, which require prior pos-tulation of the mechanism driving the heterogeneity and depend on more quan-titative parameters than homogeneous modeling. These models are affected inpart by the limited availability of quantitative data [20], but also by an incompleteunderstanding of the role played by different mechanisms driving heterogeneity.Indeed, many complex processes contribute to heterogeneity in a cell population,including gene expression noise [21], partition errors at cell division [22, 23], mu-tations [24], size variability [25], as well as environmental gradients in the culture[26].It is then important to understand what features change or are missed in amodel of continuous culture by treating the cell population as an homogeneoussystem. In this work we propose to apply the maximum entropy principle [27,28, 29, 30, 31, 32, 33] to model cell heterogeneity in a continuous bioreactor.Although we build our theory based on the model of a chemostat presented in[13], and discuss specific metabolic networks, the framework can accommodateeasily other models of continuous cell cultures and metabolic networks.4 aterials and Methods
General framework
We study the steady states of a population of cells inside a well-mixed bioreactor,where fresh medium continuously replaces culture fluid at a given dilution rate D . Each cell will be described by vectors r ( h ) and u ( h ) , giving the flux of everyreaction in the metabolic network of the cell line under study and the metaboliteuptake rates, respectively. Here and in what follows, the super-index h denotesan individual cell. If N ik denotes the stoichiometric coefficient of metabolite i inreaction k ( N ik > if metabolite i is produced in the reaction, N ik < if it is con-sumed, and N ik = 0 if it does not participate in the reaction), then cell h producesmetabolite i at an overall rate P k N ik r ( h ) k . This production must balance the cel-lular demands for metabolite i . In particular we consider a constant maintenancedemand e ( h ) i which is independent of growth [34, 35], and the requirements y ( h ) i ofeach metabolite for the synthesis of biomass [36, 37]. If biomass is synthesized ata rate z ( h ) , we obtain the following overall balance equation for each metabolite i : X k N ik r ( h ) k + u ( h ) i = e ( h ) i + y ( h ) i z ( h ) (1)For simplicity we will assume that the overall macromolecular composition of thecell is constant, i.e. that e ( h ) i = e i and y ( h ) i = y i are independent of h . The netgrowth rate of a cell depends both on the rate of biomass synthesis, z , and on theconcentration of metabolites in the media. We assume that some metabolites aretoxic (such as lactate and ammonia) and affect the growth rate of cells. In order toaccommodate the most common dependencies found in the literature, we assume: λ ( h ) = z ( h ) × K ( s ) − σ ( s ) (2)where s is the vector of metabolite concentrations in the culture, λ ( h ) is the netgrowth rate, K ( s ) is a growth inhibition factor, and σ ( s ) a death rate.In addition to Equation 1, the fluxes r ( h ) k , u ( h ) i must respect physico-chemicalconstrains arising from thermodynamics and molecular crowding. These con-strains will also depend on the extracellular conditions of the culture, specificallythe concentrations s of metabolites. They can be summarized as: X k α k | r ( h ) k | ≤ C (3)5 b k ≤ r ( h ) k ≤ ub k (4) − L i ≤ u ( h ) i ≤ U i (5)where α k are the enzymatic flux costs per unit of flux through reaction k , C is themaximum enzymatic cost available to the cell, lb k , ub k are the lower and upperbounds of reaction k , L i = 0 for metabolites that cannot be excreted from the celland L i = ∞ otherwise, U i is the maximum uptake rate of metabolite i , and X the total number of alive cells. The value of U i will depend on the conditions ofthe culture and will be specified below. The reaction bounds lb k , ub k are −∞ and ∞ respectively for reversible reactions, while lb k = 0 , ub k = ∞ for irreversiblereactions.The constrains (1–5) define a convex polytope of feasible metabolic states[38, 39] that we denote by P s . Each point within this space consists of coordinates v = ( r, u ) which fully specify the metabolic state of a cell in the model. Let P s ( v ) be the density of cells with metabolic fluxes v . To determine the form of P s ( v )d v ,we adopt the Principle of Maximum Entropy (MaxEnt), which in this context canbe stated as follows [31, 27]:Given the set of allowed metabolic states ( P s ), the dependency of thecellular growth rate with metabolic fluxes ( λ s ( v ) ), and the averagegrowth rate of the cell population h λ i , the distribution of cells within P s has the form P s ( v ) ∝ e βλ s ( v ) , where β is chosen so that the expec-tation of λ under P s ( v ) coincides with h λ i .For β = 0 every point in the space P s is equally likely. In this case cells exploreuniformly the space of allowed solutions defined by Equations (1-5). For β = ∞ the distribution is a Dirac delta with infinite mass on the points of maximumgrowth rate. In the language of Statistical Mechanics this is the ground state of thesystem, also known in the Systems Biology literature as the Flux Balance Analysis(FBA) solution [38]. Simple metabolic model
To gain some insight, we analyze first a simplified metabolic model that exhibitsa switch from oxidative pathways at low growth rates to fermentative pathwaysat higher growth rates [40]. The model consists of the following stoichiometricconstrains: v g + v l − v o = 0 (6)6 TP ATP
Figure 1:
Simplest metabolic network exhibiting a switch
A primary carbonsource (glucose) is consumed at a rate v g . It is processed into an intermediate P,generating energy (ATP). The intermediate can be excreted generating waste (W)at a rate − v l , or it can be further oxidized via the respiratory pathways (rate v o ).The respiratory pathway is capped, v o ≤ R . v g + 18 v o = v atp (7) v g ≤ min { V g , c g D/X } (8) v l ≤ , ≤ v o ≤ R (9)where v g is the velocity of consumption of glucose, v l the excretion of lactate, v o the velocity of oxidative phosphorylation, and v atp the synthesis of ATP (seeFigure 1). We only consider the ATP requirements of growth and maintenance andignore additional biomass components. Therefore the rate of biomass synthesis( z ) is an affine function of the rate of ATP synthesis, z = ( v atp − e ) /y , where y is the biomass yield and e is the constant maintenance demand for ATP [34]. Therate of oxidative phosphorylation is caped by the limited capacity of mitochondria, R [41]. This simplified model or close variations have been extensively discussedin the literature, e.g. [42, 40, 13, 43].We assume that the secreted lactate is toxic, so that the net growth rate isgiven by λ = z − σ , where the death rate σ is some increasing function of theextracellular concentration of lactate, w . For simplicity we use a linear form, σ = τ w , where τ is the death rate per unit of lactate. One can think of τ as theslope of a more realistic description of the dependence of the death rate on theconcentrations of toxic compounds, thus avoiding the introduction of too manyparameters into this simplified model. 7arameters are based on experimental values obtained for mammalian cells inthe literature. The ATP maintenance demand e = 1mmol / gDW / h has been mea-sured for mouse LS cells [34] (see [44] for similar values in other cell types). Thevalue R = 0 . mmol/gDW/h was calculated in [13], based on measurements ofthe glucose uptake threshold that triggers a switch to fermentation in mammaliancells [40]. The maximum velocity of glucose uptake, V g = 0 . mmol/gDW/h isbased on a value measured for HeLa cells [45]. The value y = 348 mmol/gDWwas then adjusted so that the maximum growth rate matches typical duplicationrates of mammalian cells of / day . From a linear approximation of the deathrate dependence on lactate concentration measured in a mammalian culture [46]we obtain τ = 0 . − mM − . The concentration of glucose was set at a valuetypical of mammalian cell culture media, c g = 15 mM. Metabolic network of CHO-K1
Motivated by the fact that most therapeutic proteins requiring complex post-translationalmodifications are produced in Chinese hamster ovary (CHO) host cell-lines [1],our second example is a metabolic model of this cell line. We employed a CHO-K1 line-specific metabolic model, based on the latest reconstruction of CHOmetabolism available at the time of writing [47]. The network recapitulates exper-imental growth rates, essential enzyme requirements and cell line specific aminoacid auxotrophies. In order to enforce (3), we complemented this network witha set of reaction costs of the form α i = M W i /a i , where M W i is the molecularweight of the enzyme catalyzing reaction i and a i the specific activity. These val-ues were mined by T. Shlomi et. al [48] from public enzyme data repositories.For reactions where the corresponding value could not be found, we set its fluxcost to the median of all the flux costs available. Finally C = 0 . / mgDW isthe mass fraction of metabolic enzymes in the dry weight of a typical mammaliancells, which can be estimated from protein abundance measurements [49, 48].The predictions of combining flux balance analysis with a crowding constrainobtained in this manner have been shown to exhibit significant correlation withenzyme mRNA expression levels in mammalian cells [48].A constant maintenance demand was added as an ATP hydrolysis flux at arate mmol/gDW/h, based on measurements for mouse LS cells [34]. The max-imum glucose uptake was set at V g = 0 . mmol/gDW/h, a value measured inHeLa cells [45]. Unfortunately kinetic parmaeters to estimate V i for many othermetabolites are not available. However, based on multiple reports in the literature[50, 46, 51] we estimated that amino acid uptakes are typically tenfold slower than8ugar uptake, and therefore set V i = V g / for amino acids. In the simulationswe used Iscove’s modified Dulbecco’s medium (IMDM) to set the values of c i ,and set V i = ∞ and infinite concentrations for water, protons and oxygen (seeSupplementary Materials for full specification). Since lactate and ammonia arethe most commonly recognized toxic byproducts in mammalian cell cultures, weset λ = K × z with K = (1 + s nh4 /K nh4 ) − (1 + s lac /K lac ) − (10)with K nh4 = 1 . mM and K lac = 8 mM, based on the values and functional formreported by [52]. Expectation propagation approximation for the computation ofmoments
In order to compute the steady state concentrations s i , it is necessary to computethe expected value of u i under the MaxEnt distribution. Although in the simplemodel this poses no difficulty, in a genome-scale metabolic network such as theCHO-K1 considered in this paper the vector v has hundreds of components and anexact computation [53] becomes intractable. A hit and run Monte Carlo approachhas been used for moderately sized networks [31, 54]. Although these methods areguaranteed to converge to an uniform sample, this is only true in the asymptoticlimit of an infinite number of steps. Unfortunately the geometrical shape of themetabolic flux space is highly elongated in some directions but very compressedin others [54]. In practice it becomes very hard to determine how long the MonteCarlo computation should run to achieve convergence, particularly so for largemetabolic networks.A better approach is to use message-passing algorithms [55]. Recently, Ex-pectation Propagation (EP) [56] has been successfully employed to compute avery good approximation of the marginal flux distributions in absence of selection( β = 0 ) [57]. In [57] the reader can find an exhaustive assessment of the qualityof this approximation in a variety of metabolic networks. In Supplementary Ma-terials we describe how the same method can be used to approximate the marginalflux distributions for non-zero β . Additional details
All model simulations were carried out in Julia [58]. The CHO-K1 metabolicnetwork [47] was loaded and setup using the COBRApy package [59]. The ex-9ectation propagation implementation was taken from https://github.com/anna-pa-m/Metabolic-EP [57].Since taking the absolute value in (3) is not a linear operation, we must replacereversible reactions in the model by two reactions, one in the forward and anotherin the backward direction. This ensures that all reaction fluxes are non-negativevariables and (3) becomes a linear inequality. However, this almost doubles thenumber of reactions in the CHO-K1 model, which makes Expectation Propagationextremely slow to converge or in some cases it even fails to do so. In order toobtain a reduced tractable network, we first solved our model at β = ∞ , whichcan be done with standard linear programming packages [13]. Then we removedfrom the CHO-K1 model all reactions that were identically zero for all values of ξ .This results in a reduced model with 380 metabolites and 401 reactions (providedas supplementary materials), where Expectation Propagation converges robustlyand fast. Results
General solution of the model
In order to determine P s we need to know the concentrations s i , the function λ s ( v ) giving the growth rate for each feasible metabolic state, and the averagegrowth rate h λ i . The later is easy in the chemostat: h λ i = D (the dilution rate)[15]. It is much more difficult to obtain experimental data for all the relevant s i . In the worst case in which no information about the s i is available, we noticethat when the chemostat is in steady state, the input flux of metabolite i mustbalance its consumption by the cells and its output flow. The input and outputfluxes per culture volume in the chemostat are given by Dc i and Ds i , respectively,where c i is the concentration of metabolite i in the external feed. Since the rate ofconsumption of individual cells is u ( h ) i , the total consumption by the populationof cells can be estimated as: X h u ( h ) i ≈ X Z u i ( v ) P s ( v )d v (11)where X is the total number of cells. In the limit of a large number of cells thisequation becomes exact. In steady state, the concentration s i is constant, which10eads to the following mass-balance equation: D ( c i − s i ) − X Z u i ( v ) P s ( v )d v = 0 (12)To obtain the steady state values of s i , equation (12) must be solved self-consistentlytogether with the MaxEnt form of the distribution P s ( v ) . Moreover since thesteady state value of s i must be non-negative, we can extract from this equa-tion an approximate form for the maximum uptake rate that further simplifies thecomputation (see Supplementary Materials for a derivation of this equation fromMichaelis-Menten kinetics), U i = min { V i , c i X/D } = min { V i , c i ξ } (13)where V i is the maximum uptake rate of metabolite i . The solution of the modelin the case β = ∞ (FBA) is described in detail in Ref. [13], where it was arguedthat the ratio ξ = X/D determines the steady state of the culture and links steadystates in the chemostat with those in perfusion systems. It has units (cells × time/ volume), and can be interpreted as the number of cells maintained alive for aperiod of time per unit of volume of fresh media supplied into the culture.In the case of finite β , the value of ξ only determines the shape of the polytopeof phenotypic states that cells can adopt. According to the maximum entropyprinciple, the distribution of cells within this polytope is of the form: P ( v ) = e β ′ z ( v ) R P e β ′ z ( v ) d v (14)where β ′ = βK ( s ) and the parameter β quantifies the level of heterogeneity inthe population of cells. It is a macroscopic representation of the underlying noisyprocesses that sustains cell-to-cell variability in the population. Small values of β lead to an almost uniform distribution P ( v ) over all possible states. This cor-responds to a highly heterogeneous population. In contrast, at larger values of β the population concentrates around the FBA solution maximizing the growth rate.This corresponds to a highly homogenous population.From (14), we compute the expected values of the exchange fluxes h u i i usingthe Expectation Propagation algorithm described in the appendix [57]. Next, thevalues of the metabolite concentrations are obtained from (12): s i = c i − h u i i ξ (15)11iven s , we compute K ( s ) , σ ( s ) and then β = β ′ /K ( s ) . Similarly, h λ i = h z i K ( s ) − σ ( s ) , which then determines the value of the dilution rate consistentwith this solution, D = h λ i . Finally, the total number of cells in the steady stateis given by X = ξD . Simple metabolic model
Within our framework the simple metabolic model admits an exact solution thatprovides important clues about the role of heterogeneity in more realistic models.For a given value of ξ , in this case the MaxEnt distribution takes the form: P ( v atp , v g ; ξ ) = e − βv atp /y /Z ( ξ ) (16)for ( v atp , v g ) ∈ P ξ , where P ξ is the polytope defined by the constrains (6)–(9)(after eliminating v l , v o ), and Z ( ξ ) = Z P ξ d v atp d v g e − βv atp /y (17)Notice that constant terms, including the additive death rate σ and the mainte-nance demands, cancel upon normalization. Due to the low-dimensionality of thismodel, the moments of (16) can be evaluated by numerical integration to any de-sired accuracy. Then the steady state concentrations of glucose and lactate can becalculated using (15).Figure 2 shows typical flux distributions of v atp and v glc in the populationof cells, for different values of ξ and β . As β increases (more homogeneouspopulations), cells cluster around the maximum feasible rate of ATP production,which also coincides with maximum glucose consumption. In the β → ∞ limit(FBA) the distribution becomes a Dirac delta (purple line) localized in this optimalpoint. On the other hand, when β → cells distribute uniformly in the space offeasible metabolic states, favoring states with low growth rate. This observationhas been interpreted as implying that higher growth rates require active regulationin the cell [31] (see Supplementary Material for a study of an evolutionary modelconsistent with this explanation within our framework).Figure 3 shows the solution of the model as a function of ξ , for selected rep-resentative values of β . For comparison the β = ∞ solution is shown in purple.As the crossing curves in Fig. 3a indicate, the effect of decreasing β is not sim-ply to decrease the average growth rate of the population of cells. At first sightthis seems to contradict the fact that the β = ∞ solution is where cells adopt12 ) b) ∞∞ Figure 2:
Distributions of the fluxes v atp and v g . Distribution among cells ofthe glycolytic flux ( v g ) and the ATP synthesis reaction flux ( v atp ). Since β hasunits of inverse time, we use the dimensionless parameter λ m β , where λ m is themaximum growth rate at ξ = 0 , β = ∞ .a phenotype with maximum growth rate (the FBA solution). However, due tothe accumulation of toxic byproducts in the culture, maximizing z may result inhigher toxicities, and the net effect is to decrease the overall growth rate of thepopulation. This translates into the fact that the heterogenous population mayhave higher cell numbers than the homogenous one (where the red curve is higherthan the purple in Fig. 3b). This explanation is confirmed in Fig. 3d, which showsthat the concentration of lactate reaches the highest levels when β = ∞ .A second feature connected to heterogeneity is that even for high values of β and ξ , the concentration of lactate in the culture increases with ξ . However in thestrictly homogeneous limit ( β = ∞ ) it goes to zero. This striking difference hasimportant implications in the interpretations of bulk metabolic measurements inpopulations of cells. For example, the observation of an increasing concentrationof lactate could be interpreted as the result of a selective pressure, pushing cellstowards fermentation. On the contrary, our model shows that in a chemostat thisis a natural consequence of the heterogeneity of the population.Moreover notice that, for high enough β , the curve of D versus ξ (Fig. 3a) dis-plays multistability. To the same value of D there may correspond more than onevalue of ξ . A theorem proved in Ref. [13] establishes that a steady state is stableif D ( ξ ) is decreasing at ξ . The theorem also holds in the present model, where theheterogeneity only redefines the function D ( ξ ) to which the theorem was origi-nally applied. In this case, an important consequence of the heterogeneity is that it13 ) b) c) d) ∞ ∞∞∞ Figure 3:
Steady states of the simple metabolic network as a function of ξ . Thedifferent panels show the steady state values of D , X , s g and s w , as functions of ξ for different values of β (indicated in the legend). Discontinuous lines indicateunstable steady states. 14 Figure 4:
Cell density versus dilution rate in the simple metabolic model.
Plotof the steady state cell density as a multi-function of the dilution rate, for differentvalues of β . Discontinuous lines indicate unstable steady states.reshapes the dynamical landscape of the system, which is tightly connected to thedecrease of the accumulation of lactate in the system as β decreases (cf. Fig. 3d).For example, increasing the heterogeneity (decreasing β ) abolishes the bistableregime. This is shown clearly in Figure 4, where the steady state values of X areplotted against the dilution rate.Heterogeneity also has the undesirable consequence of decreasing the mediumdepth ( ξ m ), the maximum value of ξ with a non-zero steady state concentrationof live cells. It is important to realize that even if the heterogeneity may helpto increase the number of cells in the presence of toxicity (as discussed in theprevious paragraphs), the medium depth never decreases with β . A plot of ξ m asa function of β is shown in Figure 5 (continuous line). For highly heterogeneoussystems (low values of β ), ξ m is almost insensitive to changes in β . Then there isa sharp slope change after which ξ m steadily increases with β .Since the stability of the system depends on the level of heterogeneity, in Fig-ure 5 we show how the values of ξ corresponding to unstable steady states dependon β (discontinuous lines). This way the ξ, β plane is divided into three regions: infeasible , where no steady state is possible, stable and unstable , according to the15 (10 cells/mL) unstablestableinfeasible Figure 5:
Critical values of ξ versus β . Medium depth (continuous line) andcritical values of ξ separating stable steady states from unstable steady states (dis-continuous lines), versus β . The gradient shows the steady state concentration ofcells ( X ). 16ype of steady state. This confirms the dramatic effect of heterogeneity on the dy-namic landscape of the system, as unstable steady states disappear below a certainthreshold value of β , i.e. after a critical level of heterogeneity. For large β , thesystem becomes more robust in the sense that the range of values of ξ definingunstable states is constant, while the stable regime becomes wider.Figure 5 also shows the density of cells in steady state for each pair of values ( β, ξ ) (color gradient). Notice that in this model, the highest cell densities occurnear unstable states. This should be interpreted as a word of caution, since tryingto increase the number of cells in the industrial setting can have the undesiredeffect of washing the culture. Analysis of a genome-scale metabolic network of the CHO-K1cell line
Next, we study a reconstruction of the CHO-K1 cell line metabolic network [47].Figure 6 shows the steady state concentrations of selected metabolites as functionsof ξ and for certain representative values of β . There are several differences be-tween the homogeneous (shown in purple in the plots) and heterogeneous regimes.Formate, a byproduct of mitochondrial oxidative metabolism secreted by normaltissues and especially cancer cells [60], stabilizes at a constant concentration forlarge ξ when β = ∞ . The selective pressure for growth predicted by an FBAcalculation only favors a mild secretion of this metabolite that is not enough tosupport its accumulation at larger ξ . In contrast, even mild levels of heterogeneityresult in an increasing accumulation of formate at larger values of ξ . Ammoniashows a similar behavior. An heterogeneous population shows an accumulationof formate and ammonia that does not originate from the selective pressures forfaster growth acting on individual cells. These differences resemble the behaviorof lactate in the simple model.In the case of histidine and valine, we have the opposite situation, where thehomogeneous model predicts non-vanishing concentrations at large ξ , but evenmild levels of heterogeneity drive the concentrations of these metabolites to zero.Other metabolites show less significant differences, such as glutamate, serine andglycine, where the level of heterogeneity only controls the rate of depletion as ξ increases but does not seem to produce any qualitative differences. Glucose andglutamine are almost insensitive to variations in β . This means that the structureof the metabolic network itself favors maximal consumption of these metabolites,even in absence of active regulation. 17 (10 cell s day/mL)10 -4 λ m β =1000 λ m β =200 λ m β =Inf λ m β =50000 glucose ξ (10 cells day/mL)10 -4 ammonia ξ (10 cells day/mL)10 -4 formate ξ (10 cells day/mL)10 -4 methionine ξ (10 cells day/mL)10 -4 glutamine g l u t a m a t e ξ (10 cells day/mL)10 -4 ξ (10 cells day/mL)10 -4 histidine ξ (10 cells day/mL)10 -4 lactate ξ (10 cells day/mL)10 -4 glycine serine ξ (10 cells day/mL)10 -4 ξ (10 cells day/mL)10 -4 valine ξ (10 cells day/mL)10 -4 p r o li n e Figure 6:
Steady state metabolite concentrations as functions of ξ for theCHO-K1 model. Steady state concentrations of selected metabolites as functionsof ξ , for the simulations of the CHO-K1 metabolic network. Representative valuesof β are plotted (see legend). The purple line corresponds to the β → ∞ limit. Toreduce clutter, this figure does not distinguish between stable and unstable steadystates. 18 ( c e ll s / m L ) ξ (10 cells day/mL)10 -4 ξ (10 cells day/mL)10 -4 X ( c e ll s / m L ) D (1/day) D ( / d a y ) λ m β =1000 λ m β =200 λ m β =Inf λ m β =50000 (c) ((cid:0)(cid:1)(cid:2)(cid:3)(cid:4) Figure 7:
Dilution rate and cell density in steady states of the CHO-K1metabolic model.
Different curves correspond to different levels of heterogene-ity. In (c), discontinuous line indicates unstable steady states.In Figure 7(a,b) we plot the dilution rate and the cell concentration in steadystate as functions of ξ , for representative values of β . The crossing of curvescorresponding to different values of β indicates that under certain conditions, het-erogeneity may enable larger quantities of cells or faster dilution rates. This issurprising because larger values of β mean that the population of cells concen-trate nearer the point of maximum growth rate. The explanation is that in this caselarger β also leads to higher secretion of toxic byproducts. For example, in Figure7 the red curve has a higher β than the black curve, but at ξ = 0 . the black curvehas higher cell counts. Comparing with Figure 6, we see that the red curve at thispoint also has higher concentrations of the toxic byproducts ammonia and lactate,explaining the reduction in cell growth. This also confirms a prediction alreadymade in the simpler model.Figure 7(c) plots the steady state values of the cell density for each dilutionrate, with unstable steady states shown in dashed lines. As in the simple metabolicmodel, we find that the system admits more than one steady state for some dilutionrates. The dynamical landscape of the system depends on the value of β . Inparticular, if β is too low the unstable regime disappears, again recapitulating thebehavior found in the simple metabolic model.Since many of the enzymatic flux costs used in these simulations ( cf. Equation193)) are unknown or poorly annotated, we repeated these simulations after doingrandom perturbations on all these coefficients of up to 25% relative to the originalvalues. Figures S3 and S4 in the Appendix show that all of the qualitative featuresdiscussed here are preserved under these perturbations.
Discussion
In this work we have developed a framework to model cell heterogeneity based onthe Maximum Entropy principle (MaxEnt) [27]. In contrast to previous applica-tions of MaxEnt to model metabolic heterogeneity [31], our work has focused oncontinuous cell cultures in steady state. As we have shown, in this case it is possi-ble to write down a self-consistent set of equations that determines the distributionof cells in phenotypic space as a function of the dilution rate. The dependence hasa non-trivial character due to the nonlinear nature of the chemostat, in some casesincluding multistable regimes.We applied this framework in a simple metabolic network first, where all thecomputations can be carried out exactly. In this toy model we were able to dis-cuss in detail many qualitative features of the model. We found that heterogeneityenables larger populations of cells because of reduced toxicity, that the dynamicallandscape includes includes a multistable regime which shrinks or even disappearsentirely as the level of heterogeneity increases, and that a byproduct predicted tobe zero if heterogeneity is ignored exhibits increasing concentrations when het-erogeneity is included in the model.In this work we have also shown how the Maximum Entropy approach tometabolic modeling can be applied to large metabolic networks, by taking advan-tage of a recent implementation of the Expectation Propagation algorithm [57].By exploiting this algorithm we were able to obtain a numerical solution of ourmodel, applied a CHO-K1 genome-scale metabolic network reconstruction. Al-though this metabolic model is much more complicated than the toy model, itnonetheless recapitulates all the qualitative findings made in the simpler network,but with a richer phenomenology.In order to obtain a tractable model, we have made some simplifications. Wechose to ignore cell-to-cell heterogeneity in mass composition (the parameters e i , y i in the notation of Methods). Indeed, although cells can vary widely in theirmetabolic arrangements, the mass fractions of major constituents such as total pro-tein, lipids, and nucleotides, are approximately invariant for a given cell type. Thisapproximation is not essential for the formal statement of the model, but it greatly20implifies calculations because then the exponent in the Maximum Entropy distri-bution is a linear function of metabolic fluxes. We have assumed that the availabil-ity of nutrients depended only on the total cell concentration. In a more realisticmodel cells compete for certain metabolites and their concentrations determinesthe uptake bounds. In the literature models of this type have been studied [61, 62]and they have an interesting phenomenology of their own. Our approach can beseen as a “mean-field” approximation where each cell interacts with the entirepopulation instead of with specific nutrients. We believe that our main conclu-sions will remain valid if this approximation is relaxed, but further work is neededin this direction. Finally, incorporating the crowding constrain entails replacingreversible reversible reactions by two reactions in the forward and backward di-rection, almost doubling the total number of reactions. We therefore applied areduction where reactions found to carry no flux identically for all ξ in the β = ∞ regime were removed from the model. Obviously this modifies the solution spaceof the model. On the other hand, an extensive model such as the genome-scaleCHO-K1 model [47] contains many reactions that are inactive under most condi-tions and therefore including them would only add noise to our analysis. SinceFBA has been a successful tool in the detection of relevant metabolic pathways inmultiple organisms [48, 47, 63, 64, 65], we believe that this reduction improvesthe quality of our predictions. Indeed, all the byproducts predicted by this modelhave been observed experimentally to be secreted by mammalian cells [13].Some general qualitative features are common to both the homogeneous regime( β = ∞ ) and the heterogeneous regime (finite β ). In both cases, multistability is aconsequence of negative feedback by toxicity accumulation, although extremeleyheterogeneous populations become monostable as shown in Figures 4, 5 and 7.For a given combination of cell-line and feed media, the ratio of cell density tothe dilution rate ( ξ = X/D , inverse of the CSPR [66]) determines the steady stateof a continuous culture and therefore can be used to connect different modes ofcontinuous cultures (chemostat and perfusion). This conclusion was obtained in[13] in the context of an homogeneous model ( β = ∞ ). In the presence of het-erogeneity the steady state does not consist of single flux values for every reactionin the network. Instead, it must be described by a global probability distributionrepresenting the fraction of cells adopting each metabolic phenotype. The param-eter β quantifies the spread of this distribution. Although following the standardpractice of the MaxEnt principle, β should be determined to match experimentaldata on the average growth rate of a population, for the sake of generality, in thiswork β was treated as a free parameter. This is analogous to studying differenttemperatures in statistical physics. Therefore the selected values of β used in the21gures carry no special significance, except that they are representative of the mostsalient aspects of the general phenomenology of the model. In the Appendix weshow how β can be connected to a simple underlying model of the evolutionarydynamics of a cell population.Interestingly, toxicity also leads to the paradoxical result that more heteroge-neous populations can achieve larger sizes. This happens because at finite β lesscells are secreting the toxic byproduct at maximal rate. Thus, a prediction of ourmodel is that inducing heterogoeneity might be beneficial in an industrial settingwhere cell numbers are limited by the accumulation of toxic byproducts. We alsoshowed that heterogeneity might be responsible for features of the bulk popula-tion not derivable from selective pressures on individual cells. The CHO-K1 re-construction exhibits divergence of formate accumulation, while for the simplifiedmetabolic model it is lactate that accumulates. This difference is not surprisingbecause in the toy model lactate is the only allowed byproduct, while the CHO-K1 model has many possible byproducts. In both cases an homogeneous modelwould predict zero or constant concentrations.Compared to population balance models of heterogeneity, our approach hasseveral advantages. First, by using MaxEnt, assumptions about the detailed mech-anisms behind cellular heterogeneity are not required. As a consequence, the de-scription of heterogeneity is simple and only one additional parameter β sufficesfor this purpose. This simplicity enables the study of genome scale metabolicnetworks such as the CHO-K1, which are not accessible to detailed populationbalance models. Acknowledgements
The authors warmly thank Anna Paola Muntoni for her assistance with the Expec-tation Propagation code. 22 eferences [1] Florian M Wurm. Production of recombinant protein therapeutics in culti-vated mammalian cells.
Nature biotechnology , 22(11):1393, 2004.[2] Rolf G Werner, Franz Walz, Wolfgang No´e, and Alois Konrad. Safety andeconomic aspects of continuous mammalian cell culture.
Journal of biotech-nology , 22(1-2):51–68, 1992.[3] JB Griffiths. Animal cell culture processes-batch or continuous?
Journal ofbiotechnology , 22(1-2):21–30, 1992.[4] A Kadouri and RE Spier. Some myths and messages concerning the batchand continuous culture of animal cells.
Cytotechnology , 24(2):89, 1997.[5] Rolf G Werner and Wolfgang No´e. Letter to the editor.
Cytotechnology ,26:81–82, 1998.[6] Matthew S Croughan, Konstantin B Konstantinov, and Charles Cooney. Thefuture of industrial bioprocessing: Batch or continuous?
Biotechnology andbioengineering , 112(4):648–651, 2015.[7] Konstantin B Konstantinov and Charles L Cooney. White paper on continu-ous bioprocessing may 20–21 2014 continuous manufacturing symposium.
Journal of pharmaceutical sciences , 104(3):813–820, 2015.[8] Anna F. Europa, Anshu Gambhir, Peng-Cheng Fu, and Wei-Shou Hu. Mul-tiple steady states with distinct cellular metabolism in continuous culture ofmammalian cells.
Biotechnology and Bioengineering , 67(1):25–34, 2000.[9] Claudia Altamirano, Andres Illanes, A. Casablancas, X. G´amez, Jordi JoanCair´o, and C. G`odia. Analysis of CHO cells metabolic redistribution ina glutamate-based defined medium in continuous culture.
BiotechnologyProgress , 17(6):1032–1041, 2001.[10] Paul M. Hayter, Elisabeth M a Curling, Anthony J. Baines, Nigel Jenk-ins, Ian Salmon, Phillip G Strange, Jeremy M Tong, and Alan T. Bull.Glucose-Limited Chemostat Culture of Chinese Hamster Ovary Cells Pro-ducing Recombinant Human Interferon- γ . Biotechnology and Bioengineer-ing , 39:327–335, 1992. 2311] Anshu Gambhir, Rashmi Korke, Jongchan Lee, Peng-Cheng Fu, Anna F. Eu-ropa, and Wei-Shou Hu. Analysis of cellular metabolism of hybridoma cellsat distinct physiological states.
Journal of bioscience and bioengineering ,95(4):317–27, 2003.[12] Brian D Follstad, R Robert Balcarcel, Gregory Stephanopoulos, andDaniel IC Wang. Metabolic flux analysis of hybridoma continuous culturesteady state multiplicity.
Biotechnology and Bioengineering , 63(6):675–683,1999.[13] Jorge Fernandez-de Cossio-Diaz, Kalet Leon, and Roberto Mulet. Charac-terizing steady states of genome-scale metabolic networks in continuous cellcultures.
PLoS Computational Biology , 13(11):e1005835, 2017.[14] Andrew Yongky, Jongchan Lee, Tung Le, Bhanu Chandra Mulukutla, Pro-dromos Daoutidis, and Wei-Shou Hu. Mechanism for multiplicity of steadystates with distinct cell concentration in continuous culture of mammaliancells.
Biotechnology and bioengineering , 112(7):1437–1445, 2015.[15] Hal L Smith and Paul Waltman.
The theory of the chemostat: dynamics ofmicrobial competition , volume 13. Cambridge university press, 1995.[16] Gyun Min Lee, Amit Varma, and Bernhard O. Palsson. Application of Pop-ulation Balance Model to the Loss of Hybridoma Antibody Productivity.
Biotechnology Progress , 7(1):72–75, January 1991.[17] Daniel J Kiviet, Philippe Nghe, Noreen Walker, Sarah Boulineau, VandaSunderlikova, and Sander J Tans. Stochasticity of metabolism and growth atthe single-cell level.
Nature , 514(7522):376, 2014.[18] Frank Delvigne, Quentin Zune, Alvaro R. Lara, Waleed Al-Soud, andSøren J. Sørensen. Metabolic variability in bioprocessing: Implications ofmicrobial phenotypic heterogeneity.
Trends in Biotechnology , 32(12):608–616, 2014.[19] Michael A Henson. Dynamic modeling of microbial cell populations.
Cur-rent Opinion in Biotechnology , 14(5):460–467, October 2003.[20] Rebeca Gonz´alez-Cabaleiro, Anca M Mitchell, Wendy Smith, Anil Wipat,and Irina D Ofit¸eru. Heterogeneity in pure microbial systems: experimentalmeasurements and modeling.
Frontiers in microbiology , 8:1813, 2017.2421] Michael B Elowitz, Arnold J Levine, Eric D Siggia, and Peter S Swain.Stochastic gene expression in a single cell.
Science , 297(5584):1183–1186,2002.[22] Jorge Fernandez-de Cossio-Diaz, Roberto Mulet, and Alexei Vazquez. Cellpopulation heterogeneity driven by stochastic partition and growth optimal-ity. arXiv preprint arXiv:1805.07768 , 2018.[23] Dann Huh and Johan Paulsson. Random partitioning of molecules at celldivision.
Proceedings of the National Academy of Sciences , 108(36):15004–15009, 2011.[24] Jiguang Wang, Emanuela Cazzato, Erik Ladewig, Veronique Frattini,Daniel IS Rosenbloom, Sakellarios Zairis, Francesco Abate, Zhaoqi Liu,Oliver Elliott, Yong-Jae Shin, et al. Clonal evolution of glioblastoma undertherapy.
Nature genetics , 48(7):768, 2016.[25] Amit Tzur, Ran Kafri, Valerie S LeBleu, Galit Lahav, and Marc W Kirschner.Cell growth and size homeostasis in proliferating animal cells.
Science ,325(5937):167–171, 2009.[26] Alvaro R Lara, Enrique Galindo, Octavio T Ram´ırez, and Laura A Palo-mares. Living with heterogeneities in bioreactors.
Molecular biotechnology ,34(3):355–381, 2006.[27] Edwin T Jaynes. Information theory and statistical mechanics.
Physicalreview , 106(4):620, 1957.[28] E. T Jaynes.
Probability Theory: The Logic of Science . Cambridge Univer-sity Press, Cambridge, UK; New York, NY, 2003.[29] John Harte and Erica A. Newman. Maximum information entropy: A foun-dation for ecological theory.
Trends in Ecology & Evolution , 29(7):384–389,July 2014.[30] Elad Schneidman, Michael J. Berry, Ronen Segev, and William Bialek.Weak pairwise correlations imply strongly correlated network states in aneural population.
Nature , 440(7087):1007–1012, April 2006.[31] Daniele De Martino, Anna Andersson, Tobias Bergmiller, C˘alin C Guet, andGaˇsper Tkaˇcik. Statistical mechanics for metabolic networks during steadystate growth.
Nature Communications , 9:2988, 2018.2532] Daniele De Martino, Fabrizio Capuani, and Andrea De Martino. Growthagainst entropy in bacterial metabolism: the phenotypic trade-off behind em-pirical growth rate distributions in e. coli.
Physical biology , 13(3):036005,2016.[33] Daniele De Martino, Fabrizio Capuani, and Andrea De Martino. Quan-tifying the entropic cost of cellular growth control.
Physical Review E ,96(1):010401, 2017.[34] DG Kilburn, MD Lilly, and FC Webb. The energetics of mammalian cellgrowth.
Journal of cell science , 4(3):645–654, 1969.[35] Kashif Sheikh, Jochen F¨orster, and Lars K Nielsen. Modeling hybridomacell metabolism using a generic genome-scale metabolic model of mus mus-culus.
Biotechnology progress , 21(1):112–121, 2005.[36] Adam M Feist and Bernhard O Palsson. The biomass objective function.
Current opinion in microbiology , 13(3):344–349, 2010.[37] Adam M Feist and Bernhard O Palsson. What do cells actually want?
Genome biology , 17(1):110, 2016.[38] Bernhard O. Palsson.
Systems Biology: Constraint-Based Reconstructionand Analysis . Cambridge University Press, Cambridge, United Kingdom,2015.[39] Matthias P. Gerstl, Stefan M ¨uller, Georg Regensburger, J¨urgen Zanghellini,and Oliver Stegle. Flux tope analysis: Studying the coordination of reactiondirections in metabolic networks.
Bioinformatics , 2018.[40] Alexei Vazquez, Jiangxia Liu, Yi Zhou, and Zolt´an N. Oltvai. Catabolicefficiency of aerobic glycolysis: The Warburg effect revisited.
BMC SystemsBiology , 4:58, May 2010.[41] Jorge Fernandez-de-Cossio-Diaz and Alexei Vazquez. Limits of aerobicmetabolism in cancer cells.
Scientific Reports , 7(1):13488, October 2017.[42] Jorge Fernandez-de-Cossio-Diaz, Andrea De Martino, and Roberto Mulet.Microenvironmental cooperation promotes early spread and bistability of aWarburg-like phenotype.
Scientific Reports , 7(1):3103, June 2017.2643] Douwe Molenaar, Rogier van Berlo, Dick de Ridder, and Bas Teusink. Shiftsin growth strategies reflect tradeoffs in cellular economics.
Molecular Sys-tems Biology , 5(1):323, January 2009.[44] Michael Lynch and Georgi K. Marinov. The bioenergetic costs of a gene.
Proceedings of the National Academy of Sciences , 112(51):15690–15695,December 2015.[45] Sara Rodr´ıguez-Enr´ıquez, Alvaro Mar´ın-Hern´andez, Juan Carlos Gallardo-P´erez, and Rafael Moreno-S´anchez. Kinetics of transport and phosphoryla-tion of glucose in cancer cells.
Journal of Cellular Physiology , 221(3):552–559, December 2009.[46] Sanjeev Dhir, K. John Morrow, R. Russell Rhinehart, and Theodore Wies-ner. Dynamic optimization of hybridoma growth in a fed-batch bioreactor.
Biotechnology and Bioengineering , 67(2):197–205, January 2000.[47] Hooman Hefzi, Kok Siong Ang, Michael Hanscho, Aarash Bordbar, DavidRuckerbauer, Meiyappan Lakshmanan, Camila A. Orellana, Deniz Baycin-Hizal, Yingxiang Huang, Daniel Ley, Veronica S. Martinez, Sarantos Kyr-iakopoulos, Natalia E. Jim´enez, Daniel C. Zielinski, Lake-Ee Quek, TuneWulff, Johnny Arnsdorf, Shangzhong Li, Jae Seong Lee, Giuseppe Paglia,Nicolas Loira, Philipp N. Spahn, Lasse E. Pedersen, Jahir M. Gutier-rez, Zachary A. King, Anne Mathilde Lund, Harish Nagarajan, AlexThomas, Alyaa M. Abdel-Haleem, Juergen Zanghellini, Helene F. Kilde-gaard, Bjørn G. Voldborg, Ziomara P. Gerdtzen, Michael J. Betenbaugh,Bernhard O. Palsson, Mikael R. Andersen, Lars K. Nielsen, Nicole Borth,Dong-Yup Lee, and Nathan E. Lewis. A Consensus Genome-scale Re-construction of Chinese Hamster Ovary Cell Metabolism.
Cell Systems ,3(5):434–443.e8, November 2016.[48] Tomer Shlomi, Tomer Benyamini, Eyal Gottlieb, Roded Sharan, and EytanRuppin. Genome-Scale Metabolic Modeling Elucidates the Role of Prolif-erative Adaptation in Causing the Warburg Effect.
PLOS ComputationalBiology , 7(3):e1002018, March 2011.[49] Christine Vogel, Raquel de Sousa Abreu, Daijin Ko, Shu-Yun Le, Bruce A.Shapiro, Suzanne C. Burns, Devraj Sandhu, Daniel R. Boutz, Edward M.27arcotte, and Luiz O. Penalva. Sequence signatures and mRNA concentra-tion can explain two-thirds of protein abundance variation in a human cellline.
Molecular Systems Biology , 6(1):400, January 2010.[50] Sadettin S. Ozturk, Mark R. Riley, and Bernhard O. Palsson. Effects ofammonia and lactate on hybridoma growth, metabolism, and antibody pro-duction.
Biotechnology and Bioengineering , 39(4):418–431, February 1992.[51] C. Altamirano, C. Paredes, J. J. Cair´o, and F. G`odia. Improvement of CHOCell Culture Medium Formulation: Simultaneous Substitution of Glucoseand Glutamine.
Biotechnology Progress , 16(1):69–75, 2000.[52] M. A. Bree, P. Dhurjati, R. F. Geoghegan, and B. Robnett. Kinetic modellingof hybridoma cell growth and immunoglobulin production in a large-scalesuspension culture.
Biotechnology and Bioengineering , 32(1067–1072):69–75, 1988.[53] D Avis and K Fukuda. A pivoting algorithm for convex hulls and vertex enu-meration of arrangements and polyhedra.
Technical Report , B(237), 1990.[54] Daniele De Martino, Matteo Mori, and Valerio Parisi. Uniform Sampling ofSteady States in Metabolic Networks: Heterogeneous Scales and Rounding.
PLOS ONE , 10(4):e0122670, April 2015.[55] Jorge Fernandez-de Cossio-Diaz and Roberto Mulet. Fast inference of ill-posed problems within a convex space.
Journal of Statistical Mechanics:Theory and Experiment , 2016(7):073207, 2016.[56] Thomas P. Minka. Expectation Propagation for Approximate Bayesian In-ference. In
Proceedings of the Seventeenth Conference on Uncertainty inArtificial Intelligence , UAI’01, pages 362–369, San Francisco, CA, USA,2001. Morgan Kaufmann Publishers Inc.[57] Alfredo Braunstein, Anna Paola Muntoni, and Andrea Pagnani. An analyticapproximation of the feasible space of metabolic networks.
Nature Commu-nications , 8:14915, April 2017.[58] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: Afresh approach to numerical computing.
SIAM review , 59(1):65–98, 2017.2859] Ali Ebrahim, Joshua A Lerman, Bernhard O Palsson, and Daniel R Hyduke.Cobrapy: Constraints-based reconstruction and analysis for python.
BMCsystems biology , 7(1):74, 2013.[60] Johannes Meiser, Anne Schuster, Matthias Pietzke, Johan Vande Voorde,Dimitris Athineos, Kristell Oizel, Guillermo Burgos-Barragan, Niek Wit,Sandeep Dhayade, Jennifer P. Morton, Emmanuel Dornier, David Sump-ton, Gillian M. Mackay, Karen Blyth, Ketan J. Patel, Simone P. Niclou, andAlexei Vazquez. Increased formate overflow is a hallmark of oxidative can-cer.
Nature Communications , 9(1):1368, April 2018.[61] Benjamin Dickens, Charles K. Fisher, and Pankaj Mehta. Analyticallytractable model for community ecology with many species.
Physical Re-view E , 94(2):022423, August 2016.[62] Christos Josephides and Peter S Swain. Predicting metabolic adaptationfrom networks of mutational paths.
Nature Communications , 8(1):685,2017.[63] Jochen F¨orster, Iman Famili, Patrick Fu, Bernhard Ø Palsson, and JensNielsen. Genome-Scale Reconstruction of the Saccharomyces cerevisiaeMetabolic Network.
Genome Research , 13(2):244–253, January 2003.[64] Laurent Heirendt, Sylvain Arreckx, Thomas Pfau, Sebasti´an N. Mendoza,Anne Richelle, Almut Heinken, Hulda S. Haraldsd´ottir, Jacek Wachowiak,Sarah M. Keating, Vanja Vlasov, Stefania Magnusd´ottir, Chiam Yu Ng, Ger-man Preciat, Alise ˇZagare, Siu H. J. Chan, Maike K. Aurich, Catherine M.Clancy, Jennifer Modamio, John T. Sauls, Alberto Noronha, Aarash Bord-bar, Benjamin Cousins, Diana C. El Assal, Luis V. Valcarcel, I˜nigo Apao-laza, Susan Ghaderi, Masoud Ahookhosh, Marouen Ben Guebila, AndrejsKostromins, Nicolas Sompairac, Hoai M. Le, Ding Ma, Yuekai Sun, LinWang, James T. Yurkovich, Miguel A. P. Oliveira, Phan T. Vuong, LemmerP. El Assal, Inna Kuperstein, Andrei Zinovyev, H. Scott Hinton, William A.Bryant, Francisco J. Arag´on Artacho, Francisco J. Planes, Egils Stalidzans,Alejandro Maass, Santosh Vempala, Michael Hucka, Michael A. Saunders,Costas D. Maranas, Nathan E. Lewis, Thomas Sauter, Bernhard Ø Palsson,Ines Thiele, and Ronan M. T. Fleming. Creation and analysis of biochemi-cal constraint-based models: The COBRA Toolbox v3.0. arXiv:1710.04038[q-bio] , October 2017. 2965] Rafael U. Ibarra, Jeremy S. Edwards, and Bernhard O. Palsson. Escherichiacoli K-12 undergoes adaptive evolution to achieve in silico predicted optimalgrowth.
Nature , 420(6912):186, November 2002.[66] Sadettin S. Ozturk. Engineering challenges in high density cell culture sys-tems.
Cytotechnology , 22(1-3):3–16, January 1996.30 upporting Information Legends • S1 Appendix.
Supplementary text • S2 Data.