Continuous cultures of mammalian cells are complex systems displaying hallmark phenomena of nonlinear dynamics, such as multi-stability, hysteresis, as well as sharp transitions between different metabolic states. In this context mathematical models may suggest control strategies to steer the system towards desired states. Although even clonal populations are known to exhibit cell-to-cell variability, most of the currently studied models assume that the population is homogeneous. To overcome this limitation, we use the maximum entropy principle to model the phenotypic distribution of cells in a chemostat as a function of the dilution rate. We consider the coupling between cell metabolism and extracellular variables describing the state of the bioreactor and take into account the impact of toxic byproduct accumulation on cell viability. We present a formal solution for the stationary state of the chemostat and show how to apply it in two examples. First, a simplified model of cell metabolism where the exact solution is tractable, and then a genome-scale metabolic network of the Chinese hamster ovary (CHO) cell line. Along the way we discuss several consequences of heterogeneity, such as: qualitative changes in the dynamical landscape of the system, increasing concentrations of byproducts that vanish in the homogeneous case, and larger population sizes.
Recombinant protein production requires suitable cell hosts and culture conditions [1]. For this purpose mammalian cells are often grown in chemostat-like cultures where a continuous flow of incoming fresh media replaces culture liquid containing cells and metabolites. Alternative processes such as batch or fed-batch are also adopted by many industrial facilities, but the advantages of the continuous mode have been predicted to drive its wide adoption in the near future [2–7]. However, experiments have demonstrated that continuous cultures exhibit hallmark phenomena of nonlinear dynamics, such as multiple steady states under identical external conditions [8–11] and hysteresis loops [8,12,13]. Sophisticated control strategies are then required to drive the system towards desired steady states. In this context, mathematical modeling has been used with some success [13–15]. Already in Ref. [13], we have shown how a model of a homogeneous continuous cell culture can explain these phenomena in the context of a detailed metabolic model, while predicting numerous metabolic transitions as a function of the ratio between cell density and dilution rate (also known as the inverse cell-specific perfusion rate [16]). However, most of these works deal with simple cell populations, consisting of identical cells (as in Ref. [13]), or at most of a few competing species [15,17]. Although it is known that no two cells in culture are alike [18], the effects of individual cell-to-cell variability are seldom considered [19]. Attempts to model heterogeneity in cell cultures are often based in population balance models [20] or similar approaches, which require prior postulation of the mechanism driving the heterogeneity and depend on more quantitative parameters than homogeneous modeling. These models are affected in part by the limited availability of quantitative data [21], but also by an incomplete understanding of the role played by different mechanisms driving heterogeneity. Indeed, many complex processes contribute to heterogeneity in a cell population, including gene expression noise [22], partition errors at cell division [23,24], mutations [25], size variability [26], as well as environmental gradients in the culture [27]. It is then important to understand what features change or are missed in a model of continuous culture by treating the cell population as an homogeneous system. In this work we propose to apply the maximum entropy principle [28–34] to model cell heterogeneity in a continuous bioreactor. This approach necessitates much less parameters than a microscopically detailed model based on population balance. We show how the model can be resolved computationally in large metabolic networks by exploiting a recent implementation of the Expectation Propagation algorithm [35,36]. To highlight the phenomenology of the heterogeneous cell population, we establish a comparison to the chemostat model of Ref. [13] which considers a homogeneous population, but we note that the framework presented here can easily accommodate other models of continuous cell cultures and metabolic networks. Although some of the predictions of Ref. [13] remain valid in the presence of heterogeneity (such as the importance of the ratio between cell density and dilution rate as a control parameter), the heterogeneity introduced here leads to important qualitative differences that we discuss in detail below. These include changes in the kinds of metabolic transitions observed (e. g. , diverging concentrations of byproducts that should be zero if the population were homogeneous or depletion of different nutrients), higher cell densities in the presence of toxicity, and drastic modifications of the dynamical landscape of the system such as disappearance of a bistable regime.
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 every reaction in the metabolic network of the cell line under study and the metabolite uptake rates, respectively. Here and in what follows, the super-index h denotes an individual cell. If Nik denotes the stoichiometric coefficient of metabolite i in reaction k (Nik > 0 if metabolite i is produced in the reaction, Nik < 0 if it is consumed, and Nik = 0 if it does not participate in the reaction), then cell h produces metabolite i at an overall rate ∑ k N i k r k (h). This production must balance the cellular demands for metabolite i. In particular we consider a constant maintenance demand e i (h) which is independent of growth [37,38], and the requirements y i (h) of each metabolite for the synthesis of biomass [39,40]. If biomass is synthesized at a rate z (h), we obtain the following overall balance equation for each metabolite i: ∑ k N i k r k (h) + u i (h) = e i (h) + y i (h) z (h) (1) For simplicity we will assume that the overall macromolecular composition of the cell is constant, i. e. that e i (h) = e i and y i (h) = y i are independent of h. The net growth rate of a cell depends both on the rate of biomass synthesis, z, and on the concentration of metabolites in the media. We assume that some metabolites are toxic (such as lactate and ammonia) and affect the growth rate of cells. In order to accommodate 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 net growth rate, K (s _) is a growth inhibition factor, and σ (s _) a death rate. In addition to Eq 1, the fluxes r k (h), u i (h) must respect physico-chemical constrains arising from thermodynamics and molecular crowding. These constrains will also depend on the extracellular conditions of the culture, specifically the concentrations s _ of metabolites. They can be summarized as: ∑ k α k | r k (h) | ≤ C (3) lb k ≤ r k (h) ≤ ub k (4) - L i ≤ u i (h) ≤ U i (5) where αk are the enzymatic flux costs per unit of flux through reaction k, C is the maximum enzymatic cost available to the cell, lbk, ubk are the lower and upper bounds of reaction k, Li = 0 for metabolites that cannot be excreted from the cell and Li = ∞ otherwise, Ui is the maximum uptake rate of metabolite i, and X the total number of alive cells. The value of Ui will depend on the conditions of the culture and will be specified below. The reaction bounds lbk, ubk are −∞ and ∞ respectively for reversible reactions, while lbk = 0, ubk = ∞ for irreversible reactions. The constrains (1–5) define a convex polytope of feasible metabolic states [41,42] 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 can be stated as follows [28,32]: For β = 0 every point in the space P s _ is equally likely. In this case cells explore uniformly the space of allowed solutions defined by Eqs (1–5). For β = ∞ the distribution is a Dirac delta with infinite mass on the points of maximum growth rate. In the language of Statistical Mechanics this is the ground state of the system, also known in the Systems Biology literature as the Flux Balance Analysis (FBA) solution [41]. To gain some insight, we analyze first a simplified metabolic model that exhibits a switch from oxidative pathways at low growth rates to fermentative pathways at higher growth rates [43]. The model consists of the following stoichiometric constrains: 2 v g + v l - v o = 0 (6) 2 v g + 18 v o = v a t p (7) v g ≤ min { V g, c g D / X } (8) v l ≤ 0,0 ≤ v o ≤ R (9) where vg is the velocity of consumption of glucose, vl the excretion of lactate, vo the velocity of oxidative phosphorylation, and vatp the synthesis of ATP (see Fig 1). We only consider the ATP requirements of growth and maintenance and ignore additional biomass components. Therefore the rate of biomass synthesis (z) is an affine function of the rate of ATP synthesis, z = (vatp − e) /y, where y is the biomass yield and e is the constant maintenance demand for ATP [37]. The rate of oxidative phosphorylation is caped by the limited capacity of mitochondria, R [44]. This simplified model or close variations have been extensively discussed in the literature, e. g. [13,43,45,46]. We assume that the secreted lactate is toxic, so that the net growth rate is given by λ = z − σ, where the death rate σ is some increasing function of the extracellular 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 the slope of a more realistic description of the dependence of the death rate on the concentrations of toxic compounds, thus avoiding the introduction of too many parameters into this simplified model. Parameters are based on experimental values obtained for mammalian cells in the literature. The ATP maintenance demand e = 1mmol/gDW/h has been measured for mouse LS cells [37] (see [47] for similar values in other cell types). The value R = 0. 45mmol/gDW/h was calculated in [13], based on measurements of the glucose uptake threshold that triggers a switch to fermentation in mammalian cells [43]. The maximum velocity of glucose uptake, Vg = 0. 5mmol/gDW/h is based on a value measured for HeLa cells [48]. The value y = 348mmol/gDW was then adjusted so that the maximum growth rate matches typical duplication rates of mammalian cells of 1/day. From a linear approximation of the death rate dependence on lactate concentration measured in a mammalian culture [49] we obtain τ = 0. 0022h−1mM−1. The concentration of glucose was set at a value typical of mammalian cell culture media, cg = 15mM. Motivated by the fact that most therapeutic proteins requiring complex post-translational modifications 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 CHO metabolism available at the time of writing [50]. The network recapitulates experimental growth rates, essential enzyme requirements and cell line specific amino acid auxotrophies. In order to enforce (3), we complemented this network with a set of reaction costs of the form αi = MWi/ai, where MWi is the molecular weight of the enzyme catalyzing reaction i and ai the specific activity. These values were mined by T. Shlomi et. al [51] from public enzyme data repositories for use in a human metabolic network. Unfortunately there is no experimental data available to determine a set of reaction costs specific for CHO-K1, but those of a human network should give a good approximation of the cost coefficients of generic mammalian cells. For reactions where the corresponding value could not be found, we set its flux cost to the median of all the flux costs available. Finally C = 0. 078mg/mgDW is the mass fraction of metabolic enzymes in the dry weight of a typical mammalian cells, which can be estimated from protein abundance measurements [51,52]. The predictions of combining flux balance analysis with a crowding constrain obtained in this manner have been shown to exhibit significant correlation with enzyme mRNA expression levels [51]. A constant maintenance demand was added as an ATP hydrolysis flux at a rate 1mmol/gDW/h, based on measurements for mouse LS cells [37]. The maximum glucose uptake was set at Vg = 0. 5mmol/gDW/h, a value measured in HeLa cells [48]. Unfortunately kinetic parameters to estimate Vi for many other metabolites are not available. However, based on multiple reports in the literature [49,53,54] we estimated that amino acid uptakes are typically tenfold slower than sugar uptake, and therefore set Vi = Vg/10 for amino acids. In the simulations we used Iscove’s modified Dulbecco’s medium (IMDM) to set the values of ci, and set Vi = ∞ and infinite concentrations for water, protons and oxygen (see Supplementary Materials for full specification). Since lactate and ammonia are the most commonly recognized toxic byproducts in mammalian cell cultures, we set λ = K × z with K = (1 + snh4/ Knh4) - 1 (1 + s lac / K lac) - 1 (10) with Knh4 = 1. 05mM and Klac = 8mM, based on the values and functional form reported by [55]. In order to compute the steady state concentrations si, it is necessary to compute the expected value of ui under the MaxEnt distribution. Although in the simple model this poses no difficulty, in a genome-scale metabolic network such as the CHO-K1 considered in this paper the vector v _ has hundreds of components and an exact computation [56] becomes intractable. A hit and run Monte Carlo approach has been used for moderately sized networks [32,57]. Although these methods are guaranteed to converge to an uniform sample, this is only true in the asymptotic limit of an infinite number of steps. Unfortunately the geometrical shape of the metabolic flux space is highly elongated in some directions but very compressed in others [57]. In practice it becomes very hard to determine how long the Monte Carlo computation should run to achieve convergence, particularly so for large metabolic networks. A better approach is to use message-passing algorithms [58]. Recently, Expectation Propagation (EP) [35] has been successfully employed to compute a very good approximation of the marginal flux distributions in absence of selection (β = 0) [36]. In [36] the reader can find an exhaustive assessment of the quality of this approximation in a variety of metabolic networks. In Supplementary Materials we describe how the same method can be used to approximate the marginal flux distributions for non-zero β. All model simulations were carried out in Julia [59]. The CHO-K1 metabolic network [50] was loaded and setup using the COBRApy package [60]. The expectation propagation implementation was taken from https: //github. com/anna-pa-m/Metabolic-EP [36]. Since taking the absolute value in (3) is not a linear operation, we must replace reversible reactions in the model by two reactions, one in the forward and another in the backward direction. This ensures that all reaction fluxes are non-negative variables and (3) becomes a linear inequality. However, this almost doubles the number of reactions in the CHO-K1 model, which makes Expectation Propagation extremely slow to converge or in some cases it even fails to do so. In order to obtain a reduced tractable network, we first solved our model at β = ∞, which can be done with standard linear programming packages [61]. In cases where FBA produces multiple solutions maximizing biomass production [62], we select the one that minimizes the enzymatic cost in (3) [63]. This guarantees that the solution is uniquely determined, because the cost coefficients are sufficiently heterogeneous that the contours of constant cost are almost surely not parallel to the contours of constant biomass production. Moreover from the biological viewpoint it is reasonable that the cell minimizes the enzymatic mass required to support a certain metabolic goal [63–65]. Then we removed from 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 (provided as supplementary materials), where Expectation Propagation converges robustly and fast.
In order to determine P s _ we need to know the concentrations si, the function λ s _ (v _) giving the growth rate for each feasible metabolic state, and the average growth rate 〈λ〉. The later is easy in the chemostat: 〈λ〉 = D (the dilution rate) [15]. It is much more difficult to obtain experimental data for all the relevant si. In the worst case in which no information about the si is available, we notice that when the chemostat is in steady state, the input flux of metabolite i must balance its consumption by the cells and its output flow. The input and output fluxes per culture volume in the chemostat are given by Dci and Dsi, respectively, where ci is the concentration of metabolite i in the external feed. Since the rate of consumption of individual cells is u i (h), the total consumption by the population of cells can be estimated as: ∑ h u i (h) ≈ X ∫ 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 this equation becomes exact. In steady state, the concentration si is constant, which leads to the following mass-balance equation: D (c i - s i) - X ∫ u i (v _) P s _ (v _) d v _ = 0 (12) To obtain the steady state values of si, Eq (12) must be solved self-consistently together with the MaxEnt form of the distribution P s _ (v _). Moreover since the steady state value of si must be non-negative, we can extract from this equation an approximate form for the maximum uptake rate that further simplifies the computation (see Supplementary Materials for a derivation of this equation from Michaelis-Menten kinetics), U i = min { V i, c i X / D } = min { V i, c i ξ } (13) where Vi is the maximum uptake rate of metabolite i. The solution of the model in the case β = ∞ (FBA) is described in detail in Ref. [13], where it was argued that the ratio ξ = X/D determines the steady state of the culture and links steady states 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 a period 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 polytope of phenotypic states that cells can adopt. According to the maximum entropy principle, the distribution of cells within this polytope is of the form: P (v _) = e β ′ z (v _) ∫ P e β ′ z (v _) d v _ (14) where β ′ = β K (s _) and the parameter β quantifies the level of heterogeneity in the population of cells. It is a macroscopic representation of the underlying noisy processes 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 corresponds 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 〈ui〉 using the Expectation Propagation algorithm described in the appendix [36]. Next, the values of the metabolite concentrations are obtained from (12): s i = c i - ⟨ u i ⟩ ξ (15) Given s _, we compute K (s _), σ (s _) and then β = β ′ / K (s _). Similarly, 〈 λ 〉 = 〈 z 〉 K (s _) - σ (s _), which then determines the value of the dilution rate consistent with this solution, D = 〈λ〉. Finally, the total number of cells in the steady state is given by X = ξD. Within our framework the simple metabolic model admits an exact solution that provides 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 a t p, v g; ξ) = e - β v a t p / y / Z (ξ) (16) for (v a t p, v g) ∈ P ξ, where P ξ is the polytope defined by the constrains (6) – (9) (after eliminating vl, vo), and Z (ξ) = ∫ P ξ d v a t p d v g e - β v a t p / y (17) Notice that constant terms, including the additive death rate σ and the maintenance demands, cancel upon normalization. Due to the low-dimensionality of this model, the moments of (16) can be evaluated by numerical integration to any desired accuracy. Then the steady state concentrations of glucose and lactate can be calculated using (15). Fig 2 shows typical flux distributions of vatp and vglc in the population of cells, for different values of ξ and β. As β increases (more homogeneous populations), 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 optimal point. On the other hand, when β → 0 cells distribute uniformly in the space of feasible metabolic states, favoring states with low growth rate. This observation has been interpreted as implying that higher growth rates require active regulation in the cell [32] (see Supplementary Material for a study of an evolutionary model consistent with this explanation within our framework). Fig 3 shows the solution of the model as a function of ξ, for selected representative values of β. For comparison the β = ∞ solution is shown in purple. As the crossing curves in Fig 3a indicate, the effect of decreasing β is not simply to decrease the average growth rate of the population of cells. At first sight this seems to contradict the fact that the β = ∞ solution is where cells adopt a phenotype with maximum growth rate (the FBA solution). However, due to the accumulation of toxic byproducts in the culture, maximizing z may result in higher toxicities, and the net effect is to decrease the overall growth rate of the population. This translates into the fact that the heterogenous population may have higher cell numbers than the homogenous one (where the red curve is higher than the purple in Fig 3b). This explanation is confirmed in Fig 3d, which shows that 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 the strictly homogeneous limit (β = ∞) it goes to zero. This striking difference has important implications in the interpretations of bulk metabolic measurements in populations of cells. For example, the observation of an increasing concentration of lactate could be interpreted as the result of a selective pressure, pushing cells towards fermentation. On the contrary, our model shows that in a chemostat this is a natural consequence of the heterogeneity of the population. Moreover notice that, for high enough β, the curve of D versus ξ (Fig 3a) displays multistability. To the same value of D there may correspond more than one value of ξ. A theorem proved in Ref. [13] establishes that a steady state is stable if D (ξ) is decreasing at ξ. The theorem also holds in the present model, where the heterogeneity only redefines the function D (ξ) to which the theorem was originally applied. In this case, an important consequence of the heterogeneity is that it reshapes the dynamical landscape of the system, which is tightly connected to the decrease of the accumulation of lactate in the system as β decreases (cf. Fig 3d). For example, increasing the heterogeneity (decreasing β) abolishes the bistable regime. This is shown clearly in Fig 4, where the steady state values of X are plotted against the dilution rate. Heterogeneity also has the undesirable consequence of decreasing the medium depth (ξm), the maximum value of ξ with a non-zero steady state concentration of live cells. It is important to realize that even if the heterogeneity may help to increase the number of cells in the presence of toxicity (as discussed in the previous paragraphs), the medium depth never decreases with β. A plot of ξm as a function of β is shown in Fig 5 (continuous line). For highly heterogeneous systems (low values of β), ξm is almost insensitive to changes in β. Then there is a sharp slope change after which ξm steadily increases with β. Since the stability of the system depends on the level of heterogeneity, in Fig 5 we show how the values of ξ corresponding to unstable steady states depend on β (discontinuous lines). This way the ξ, β plane is divided into three regions: infeasible, where no steady state is possible, stable and unstable, according to the type of steady state. This confirms the dramatic effect of heterogeneity on the dynamic landscape of the system, as unstable steady states disappear below a certain threshold value of β, i. e. after a critical level of heterogeneity. For large β, the system becomes more robust in the sense that the range of values of ξ defining unstable states is constant, while the stable regime becomes wider. Fig 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 occur near unstable states. This should be interpreted as a word of caution, since trying to increase the number of cells in the industrial setting can have the undesired effect of washing the culture. Next, we study a reconstruction of the CHO-K1 cell line metabolic network [50]. Fig 6 shows the steady state concentrations of selected metabolites as functions of ξ and for certain representative values of β. There are several differences between the homogeneous (shown in purple in the plots) and heterogeneous regimes. Formate, a byproduct of mitochondrial oxidative metabolism secreted by normal tissues and especially cancer cells [66], stabilizes at a constant concentration for large ξ when β = ∞. The selective pressure for growth predicted by an FBA calculation only favors a mild secretion of this metabolite that is not enough to support its accumulation at larger ξ. In contrast, even mild levels of heterogeneity result in an increasing accumulation of formate at larger values of ξ. Ammonia shows a similar behavior. An heterogeneous population shows an accumulation of formate and ammonia that does not originate from the selective pressures for faster growth acting on individual cells. These differences resemble the behavior of lactate in the simple model. In the case of histidine and valine, we have the opposite situation, where the homogeneous model predicts non-vanishing concentrations at large ξ, but even mild levels of heterogeneity drive the concentrations of these metabolites to zero. Other metabolites show less significant differences, such as glutamate, serine and glycine, where the level of heterogeneity only controls the rate of depletion as ξ increases but does not seem to produce any qualitative differences. Glucose and glutamine are almost insensitive to variations in β. This means that the structure of the metabolic network itself favors maximal consumption of these metabolites, even in absence of active regulation. In Fig 7 (a) and 7 (b) we plot the dilution rate and the cell concentration in steady state as functions of ξ, for representative values of β. The crossing of curves corresponding to different values of β indicates that under certain conditions, heterogeneity may enable larger quantities of cells or faster dilution rates. This is surprising because larger values of β mean that the population of cells concentrate nearer the point of maximum growth rate. The explanation is that in this case larger β also leads to higher secretion of toxic byproducts. For example, in Fig 7 the red curve has a higher β than the black curve, but at ξ = 0. 1 the black curve has higher cell counts. Comparing with Fig 6, we see that the red curve at this point also has higher concentrations of the toxic byproducts ammonia and lactate, explaining the reduction in cell growth. This also confirms a prediction already made in the simpler model. Fig 7 (c) plots the steady state values of the cell density for each dilution rate, with unstable steady states shown in dashed lines. As in the simple metabolic model, we find that the system admits more than one steady state for some dilution rates. The dynamical landscape of the system depends on the value of β. In particular, if β is too low the unstable regime disappears, again recapitulating the behavior found in the simple metabolic model. Since many of the enzymatic flux costs used in these simulations (cf. Eq (3) ) are unknown or poorly annotated, we repeated these simulations after doing random perturbations on all these coefficients of up to 25% relative to the original values. Figures S3 and S4 in the Appendix show that all of the qualitative features discussed here are preserved under these perturbations.
In this work we have developed a framework to model cell heterogeneity based on the Maximum Entropy principle (MaxEnt) [28]. In contrast to previous applications of MaxEnt to model metabolic heterogeneity [32], our work has focused on continuous cell cultures in steady state. As we have shown, in this case it is possible to write down a self-consistent set of equations that determines the distribution of cells in phenotypic space as a function of the dilution rate. The dependence has a non-trivial character due to the nonlinear nature of the chemostat, in some cases including multistable regimes. We applied this framework in a simple metabolic network first, where all the computations can be carried out exactly. In this toy model we were able to discuss in detail many qualitative features of the model. We found that heterogeneity enables larger populations of cells because of reduced toxicity, that the dynamical landscape includes includes a multistable regime which shrinks or even disappears entirely as the level of heterogeneity increases, and that a byproduct predicted to be zero if heterogeneity is ignored exhibits increasing concentrations when heterogeneity is included in the model. In this work we have also shown how the Maximum Entropy approach to metabolic modeling can be applied to large metabolic networks, by taking advantage of a recent implementation of the Expectation Propagation algorithm [36]. By exploiting this algorithm we were able to obtain a numerical solution of our model, applied to a CHO-K1 genome-scale metabolic network reconstruction. Although this metabolic model is much more complicated than the toy model, it nonetheless 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. We chose to ignore cell-to-cell heterogeneity in mass composition (the parameters ei, yi in the notation of Methods). Indeed, although cells can vary widely in their metabolic arrangements, the mass fractions of major constituents such as total protein, lipids, and nucleotides, are approximately invariant for a given cell type. This approximation is not essential for the formal statement of the model, but it greatly simplifies calculations because then the exponent in the Maximum Entropy distribution is a linear function of metabolic fluxes. We have assumed that the availability of nutrients depended only on the total cell concentration. In a more realistic model cells compete for certain metabolites and their concentrations determines the uptake bounds. In the literature models of this type have been studied [67,68] and they have an interesting phenomenology of their own. Our approach can be seen as a “mean-field” approximation where each cell interacts with the entire population instead of with specific nutrients. We believe that our main conclusions will remain valid if this approximation is relaxed, but further work is needed in this direction. Finally, incorporating the crowding constrain entails replacing reversible reversible reactions by two reactions in the forward and backward direction, almost doubling the total number of reactions. We therefore applied a reduction where reactions found to carry no flux identically for all ξ in the β = ∞ regime were removed from the model. Obviously this modifies the solution space of the model. On the other hand, an extensive model such as the genome-scale CHO-K1 model [50] contains many reactions that are inactive under most conditions and therefore including them would only add noise to our analysis. Since FBA has been a successful tool in the detection of relevant metabolic pathways in multiple organisms [50,51,69–71], we believe that this reduction improves the quality of our predictions. Indeed, all the byproducts predicted by this model have 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 a consequence of negative feedback by toxicity accumulation, although extremely heterogeneous populations become monostable as shown in Figs 4,5 and 7. For a given combination of cell-line and feed media, the ratio of cell density to the dilution rate (ξ = X/D, inverse of the cell-specific perfusion rate [16]) determines the steady state of a continuous culture and therefore can be used to connect different modes of continuous cultures (chemostat and perfusion). This conclusion was obtained in [13] in the context of an homogeneous model (β = ∞). In the presence of heterogeneity the steady state does not consist of single flux values for every reaction in the network. Instead, it must be described by a global probability distribution representing the fraction of cells adopting each metabolic phenotype. The parameter β quantifies the spread of this distribution. Although following the standard practice of the MaxEnt principle, β should be determined to match experimental data on the average growth rate of a population, for the sake of generality, in this work β was treated as a free parameter. This is analogous to studying different temperatures in statistical physics. Therefore the selected values of β used in the figures carry no special significance, except that they are representative of the most salient aspects of the general phenomenology of the model. In the Appendix we show how β can be connected to a simple underlying model of the evolutionary dynamics of a cell population. Interestingly, toxicity also leads to the paradoxical result that more heterogeneous populations can achieve larger sizes. This happens because at finite β less cells are secreting the toxic byproduct at maximal rate. Thus, a prediction of our model is that inducing heterogeneity might be beneficial in an industrial setting where cell numbers are limited by the accumulation of toxic byproducts. We also showed that heterogeneity might be responsible for features of the bulk population not derivable from selective pressures on individual cells. The CHO-K1 reconstruction exhibits divergence of formate accumulation, while for the simplified metabolic model it is lactate that accumulates. This difference is not surprising because in the toy model lactate is the only allowed byproduct, while the CHO-K1 model has many possible byproducts. In both cases an homogeneous model would predict zero or constant concentrations. Compared to population balance models of heterogeneity, our approach has several advantages. First, by using MaxEnt, assumptions about the detailed mechanisms behind cellular heterogeneity are not required. As a consequence, the description of heterogeneity is simple and only one additional parameter β suffices for this purpose. This simplicity enables the study of genome scale metabolic networks such as the CHO-K1, which are not accessible to detailed population balance models.
While the advantages of continuous culture in the biotechnological industry have been widely advocated in the literature, its adoption over batch or fed-batch modes stalls due to the complexities of these systems. In particular, continuous cell cultures display hallmark nonlinear phenomena such as multi-stability, hysteresis, and sharp transitions between metabolic phenotypes. Moreover, the impact of the heterogeneity of a cell population on these features is not well understood. We use the maximum entropy principle to model the phenotypic distribution of an heterogeneous population of cells in a chemostat. Given the metabolic network and the dilution rate, we obtain a self-consistent solution for the stationary distribution of metabolic fluxes in cells. We apply the formalism in two examples: a simplified model where the exact solution is tractable, and a genome-scale metabolic network of the Chinese hamster ovary (CHO) cell line widely used in industry. We demonstrate that heterogeneity may be responsible for qualitative changes in the dynamical landscape of the system, like the disappearance of a bistable regime, the increase of concentrations of byproducts that vanish in the homogeneous system and larger number of cells. We explain the causes behind these phenomena.
cell physiology
carbohydrate metabolism
protein metabolism
biological cultures
metabolic networks
enzymology
cell metabolism
glucose metabolism
metabolites
network analysis
enzyme metabolism
cell cultures
thermodynamics
enzyme chemistry
research and analysis methods
computer and information sciences
entropy
physics
biochemistry
cell biology
biology and life sciences
physical sciences
metabolism
Maximum entropy and population heterogeneity in continuous cell cultures