Search tips
Search criteria 


Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
PLoS Comput Biol. 2017 November; 13(11): e1005835.
Published online 2017 November 13. doi:  10.1371/journal.pcbi.1005835
PMCID: PMC5703580

Characterizing steady states of genome-scale metabolic networks in continuous cell cultures

Jorge Fernandez-de-Cossio-Diaz, Conceptualization, Investigation, Writing – original draft, Writing – review & editing,1,2,* Kalet Leon, Conceptualization, Methodology, Supervision, Writing – review & editing,1 and Roberto Mulet, Conceptualization, Methodology, Supervision, Writing – review & editing2
Jens Nielsen, Editor


In the continuous mode of cell culture, a constant flow carrying fresh media replaces culture fluid, cells, nutrients and secreted metabolites. Here we present a model for continuous cell culture coupling intra-cellular metabolism to extracellular variables describing the state of the bioreactor, taking into account the growth capacity of the cell and the impact of toxic byproduct accumulation. We provide a method to determine the steady states of this system that is tractable for metabolic networks of arbitrary complexity. We demonstrate our approach in a toy model first, and then in a genome-scale metabolic network of the Chinese hamster ovary cell line, obtaining results that are in qualitative agreement with experimental observations. We derive a number of consequences from the model that are independent of parameter values. The ratio between cell density and dilution rate is an ideal control parameter to fix a steady state with desired metabolic properties. This conclusion is robust even in the presence of multi-stability, which is explained in our model by a negative feedback loop due to toxic byproduct accumulation. A complex landscape of steady states emerges from our simulations, including multiple metabolic switches, which also explain why cell-line and media benchmarks carried out in batch culture cannot be extrapolated to perfusion. On the other hand, we predict invariance laws between continuous cell cultures with different parameters. A practical consequence is that the chemostat is an ideal experimental model for large-scale high-density perfusion cultures, where the complex landscape of metabolic transitions is faithfully reproduced.

Author summary

While at present most biotechnology industrial facilities adopt batch or fed-batch processes, continuous processing has been vigorously defended in the literature and many predict its adoption in the near future. However, identical cultures may lead to distinct steady states and the lack of comprehension of this multiplicity has been a limiting factor for the widespread application of this kind of processes in the industry. In this work we try to remediate this providing a computationally tractable approach to determine the steady-states of genome-scale metabolic networks in continous cell cultures and show the existence of general invariance laws across different cultures. We represent a continuous cell culture as a metabolic model of a cell coupled to a dynamic environment that includes toxic by-products of metabolism and the cell capacity to grow. We show that the ratio between cell density and dilution rate is the control parameter fixing steady states with desired properties, and that this is invariant accross perfusion systems. The typical multi-stability of the steady-states of this kind of culture is explained by the negative feedback loop on cell growth due to toxic byproduct accumulation. Moreover, we present invariance laws connecting continuous cell cultures with different parameters that imply that the chemostat is the ideal experimental model to faithfully reproduce the complex landscape of metabolic transitions of a perfusion system.


Biotechnological products are obtained by treating cells as little factories that transform substrates into products of interest. There are three major modes of cell culture: batch, fed-batch and continuous. In batch, cells are grown with a fixed initial pool of nutrients until they starve, while in fed-batch the pool of nutrients is re-supplied at discrete time intervals. Cell cultures in the continuous mode are carried out with a constant flow carrying fresh medium replacing culture fluid, cells, unused nutrients and secreted metabolites, usually maintaining a constant culture volume. While at present most biotechnology industrial facilities adopt batch or fed-batch processes, the advantages of continuous processing have been vigorously defended in the literature [15], and currently some predict its widespread adoption in the near future [6].

A classical example of continuous cell culture is the chemostat, invented in 1950 independently by Aaron Novick and Leo Szilard [7] (who also coined the term chemostat) and by Jacques Monod [8]. In this system, microorganisms reside inside a vessel of constant volume, while sterile media, containing nutrients essential for cell growth, is delivered at a constant rate. Culture medium containing cells, remanent substrates and products secreted by the cells are removed at the same rate, maintaining a constant culture volume. The main dynamical variable in this system is the dilution rate (D), which is the rate at which culture fluid is replaced divided by the culture volume. In a well-stirred tank any entity (molecule or cell) has a probability per unit time D of leaving the vessel. In industrial settings, higher cell densities are achieved by attaching a cell retention device to the chemostat, but allowing a bleeding rate to remove cell debris [9]. Effectively only a fraction 0 ≤ ϕ ≤ 1 of cells are carried away by the output flow D. This variation of the continuous mode is known as perfusion culture.

By definition, a continuous cell culture ideally reaches a steady state when the macroscopic properties of the tank (cell density and metabolite concentrations) attain stationary values. Industrial applications place demands on the steady state, usually: high-cell density, minimum waste byproduct accumulation, and efficient nutrient use. However, identical external conditions (dilution rate, media formulation) may lead to distinct steady states with different metabolic properties (a phenomenon known in the literature as multi-stability or multiplicity of steady states) [1014]. Therefore, for the industry, it becomes fundamental to know in advance, given the cell of interest and the substrates to be used, which are the possible steady states of the system and how to reach them. Moreover, to satisfy production demands, it may be advantageous to extend the duration of a desired steady state indefinitely [6], implying that their stability properties are also of great interest.

Fortunately, in the last few years it has been possible to exploit an increasingly available amount of information about cellular metabolism at the stoichiometric level to build genome-scale metabolic networks [15, 16]. These networks have been modeled by different approaches [17, 18] but Flux Balance Analysis (FBA) has been particularly successful predicting cell metabolism in the growth phase [19]. FBA starts assuming a quasi-steady state of intra-cellular metabolite concentrations, which is easily translated into a linear system of balance equations to be satisfied by reaction fluxes. This system of equations is under-determined and a biologically motivated metabolic objective, such as biomass synthesis, is usually optimized to determine the complete distribution of fluxes through the solution of a Linear Programming problem [20]. This approach was first used to characterize the metabolism of bacterial growth [21], but later has been applied also to eukaryotic cells [22, 23]. Alternatively, given a set of under-determined linear equations, one can estimate the space of feasible solutions of the system and average values of the reaction fluxes [2426].

To consider the temporal evolution of a culture, FBA may be applied to successive points in time, coupling cell metabolism to the dynamics of extra-cellular concentrations. This is the approach of Dynamic Flux Balance Analysis (DFBA) [27] and has been applied prominently either to the modeling of batch/fed-batch cultures or to transient responses in continuous cultures, being particularly successful in predicting metabolic transitions in E. Coli and yeast [23, 27, 28]. However, to the best of our knowledge, the steady states of continuous cell cultures have not been investigated before. First, because DBFA for genome-scale metabolic networks may be a computational demanding task, particularly when the interest is to understand long-time behavior. Second, because it assumes knowledge of kinetic parameters describing metabolic exchanges between the cell and culture medium, that are usually unknown in realistic networks. Moreover, although the importance of toxic byproduct accumulation has been appreciated for decades [29, 30], its impact on steady states of continuous cultures has been studied mostly in simple metabolic models involving few substrates [31, 32], while it has been completely overlooked in DFBA of large metabolic networks. Lactate and ammonia are the most notable examples in this regard and have been widely studied in experiments in batch and continuous cultures [30, 3336].

Our goal in this work is to introduce a detailed characterization of the steady states of cell cultures in continuous mode, considering the impact of toxic byproduct accumulation on the culture, and employing a minimum number of essential kinetic parameters. To achieve this and inspired by the success of DFBA in other settings we couple macroscopic variables of the bioreactor (metabolite concentrations, cell density) to intracellular metabolism. However, we explain how to proceed directly to the determination of steady states, bypassing the necessity of solving the dynamical equations of the problem. This spares us from long simulation times and provides an informative overview of the dynamic landscape of the system. The approach, presented here for a toy model and for a genome-scale metabolic network of CHO-K1, but easily extensible to other systems, supports the idea that multi-stability, i.e., the coexistence of multiple steady states under identical external conditions, arises as a consequence of toxic byproduct accumulation in the culture. We find and characterize specific transitions, defined by simultaneous changes in the effective cell growth rate and metabolic states of the cell, and find a wide qualitative agreement with experimental results in the literature. Our analysis implies that batch cultures, typically used as benchmarks of cell-lines and culture media, are unable to characterize the landscape of metabolic transitions exhibited by perfusion systems. On the other hand, our results suggest a general scaling law that translates between the steady states of a chemostat and any perfusion system. Therefore, we predict that the chemostat is an ideal experimental model of high-cell density perfusion cultures, enabling a faithful characterization of the performance of a cell-line and media formulation truly valid in perfusion systems.

Materials and methods

Dynamical model of the perfusion system

We study an homogeneous population of cells growing inside a well-mixed bioreactor [37], where fresh medium continuously replaces culture fluid (Fig 1). The fundamental dynamical equations describing this system are:



where X denotes the density of cells in the bioreactor (units: gDW/L), μ the effective cell growth rate (units: 1/hr), ui the specific uptake of metabolite i (units: mmol/gDW/hr), and si the concentration of metabolite i in the culture (units: mM). The external parameters controlling the culture are the medium concentration of metabolite i, ci, the dilution rate, D (units: 1 / time), and the bleeding coefficient, ϕ (unitless), which in perfusion systems characterizes the fraction of cells that escape from the culture through a cell-retention device [9] or a bleeding rate. For convenience of notation, in what follows an underlined symbol like s_ will denote the vector with components {si}.

Fig 1
Cell culture in continuous mode.

Eq 1 describes the dynamics of the cell density as a balance between cell growth and dilution, while Eq 2 describes the dynamics of metabolite concentrations in the culture as a balance between cell consumption (or excretion if ui < 0) and dilution. One must notice that at variance with the standard formulation of DFBA, the terms involving the dilution rate in the right-hand side of both equations enable the existence of non-trivial steady states (with non-zero cell density) which are impossible in batch. These are the steady states that are relevant for industrial applications adopting the perfusion model and that we study in what follows.

Still, we require a functional connection between variables describing the macroscopic state of the tank (X, s_) and the average behavior of cells (u_, μ). We start assuming that metabolites inside the cell attain quasi-steady state concentrations [38], so that fluxes of intra-cellular metabolic reactions balance at each metabolite. 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 rk is the flux of reaction k, then the metabolic network produces a net output flux of metabolite i at a rate ∑k Nik rk, where Nik = 0 if metabolite i does not participate in reaction k. This output flux must balance the cellular demands for metabolite i. In particular we consider a constant maintenance demand at rate ei which is independent of growth [39, 40], as well as the requirements of each metabolite for the synthesis of biomass components. If yi units of metabolite i are needed per unit of biomass produced [41, 42], and biomass is synthesized at a rate z, we obtain the following overall balance equation for each metabolite i:


It is also well known that a cell has a limited enzymatic budget [17]. The synthesis of new enzymes, needed to catalyze many intracellular reactions, consumes limited resources, including amino acids, energy, cytosolic [22, 43, 44] or membrane space [45] (for enzymes located on membranes), ribosomes [46], all of which can be modeled as generic enzyme costs [17]. We split reversible reaction fluxes into negative and positive parts, rk=rk+-rk-, with rk±0, and quantify the total cost of a flux distribution in the simplest (approximate) linear form [17]:


where αk+,αk- are constant flux costs. The limited budget of the cell to support enzymatic reactions is modeled as a constraint αC, where C is a constant maximum cost. Thermodynamics places additional reversibility constrains on the flux directions of some intra-cellular reactions [47], which can be written as:

lbk ≤ rk ≤ ubk,  lbk, ubk ∈ { - , 0, }.

Additionally, some uptakes ui are limited by the availability of nutrients in the culture. We distinguish two regimes. If the cell density is low, nutrients will be in excess and uptakes are only bounded by the intrinsic kinetics of cellular transporters. In this case uiVi, where Vi is a constant maximum uptake rate determined by molecular details of the transport process. These will be the only kinetic parameters introduced in the model. When the cell density increases and the concentrations of limiting substrates reach very low levels, a new regime appears where cells compete for resources. In this regime the natural condition si ≥ 0 together with the mass balance equation (Eq 2 in steady state) imply that uici D/X. In summary,

Li ≤ ui ≤ min {ViciD/X}

where Li = 0 for metabolites that cannot be secreted, and Li = ∞ otherwise. Thus, an important approximation in our model is that low concentrations of limiting nutrients are replaced by an exact zero. The ratio D/X in Eq 6 establishes the desired coupling between internal metabolism and external variables of the bioreactor. The appendix contains an alternative derivation of Eq 6.

Next, we reason that, although cellular clones in biotechnology are artificially chosen according to various productivity-related criteria [48], the growth rate is typically under an implicit selective pressure. We will consider then that the flux distribution of metabolic reactions inside the cell maximizes the rate of biomass synthesis, z, subject to all the constrains enumerated above. Note that to carry out this optimization it is enough to solve a linear programming problem, for which efficient algorithms are available [49]. This formulation is closely related to Flux Balance Analysis (FBA) [20, 21, 50, 51], but some of the constrains imposed here might be unfamiliar. In particular, Eq 4 has been used before to explain switches between high-yield and high-rate metabolic modes under the name FBA with molecular crowding (FBAwMC) [43, 52, 53], while the right-hand side of Eq 6 is a novel constraint introduced in this work to model continuous cell cultures. If multiple metabolic flux distributions are consistent with a maximal biomass synthesis rate [54], the one with minimum cost α (cf. Eq 4) is selected [17]. Summarizing, from the complete solution of the linear program we obtain the optimal z, and the metabolic fluxes u_ feeding the synthesis of biomass.

Finally, the net growth rate of cells μ (see Eq 1) is essentially determined by the cellular capacity to synthesize biomass (rate z), but it may also be affected by environmental toxicity. In the examples presented below we considered that: μ=z-σ(s_) or μ=z×K(s_), corresponding to two different mechanisms explored in the literature [36, 55]. In the first case σ(s_) is easily interpreted as the death rate of the cell, while 01-K(s_)1 represents a fraction of biomass that must be expended on non-growth related activities, for example, due to increased maintenance demands on account of environmental toxicity (but see also Refs. [56, 57] and in particular B. Ben Yahia et al. [37] for a recent review of the subject). Both σ(s_) and K(s_) depend on the concentrations of toxic metabolites in the culture, such as lactate and ammonia.

Metabolic networks

Toy model

To gain insight into the kind of solutions expected, we examined first a simple metabolic network that admits an analytic solution. It is based on the simplified network studied by A. Vazquez et al. to explain the Warburg effect [58], and serves as a minimal model of metabolic transitions in the cell [5860].

A diagram of the network is shown in Fig 2. There are four metabolites: a primary nutrient S, an energetic currency E, an intermediate P, and a waste product W. Only S and W can be exchanged with the extracellular medium, and their concentrations in the tank will be denoted by s and w, respectively. The cell can consume S from the medium at a rate u ≥ 0. The nutrient is first processed into P, generating NF units of E per unit S processed. The intermediate can have two destinies: it can be excreted in the form of W (rate −v ≤ 0), or it can be further oxidized (rate r ≥ 0), generating NR units of E per unit P. These two pathways are reminiscent of fermentation and respiration. We assume that NF [double less-than sign] NR, which is consistent with the universally lower energy yield of fermentation versus respiration. Therefore, a maximization of energy output implies that the respiration mode is preferred. However, the enzymatic costs required to enable respiration are very high compared to fermentation. Therefore in Eq 4 only the costs of respiration are significant [58], which implies that this flux is bounded:

r ≤ rmax

Metabolic overflow occurs when the nutrient uptake is higher than the respiratory bound rmax. The remaining S must then be exported as waste, W. A balance constraint (cf. Eq 3) at the intermediate metabolite P requires that:

uvr = 0

where stoichiometric coefficients are set to 1 for simplicity. Another balance constraint at the internal energetic currency metabolite, E, leads to:

NFuNRreyz = 0

where e denotes an energetic maintenance demand. The currency E is a direct precursor of biomass, at a yield y. Finally, the waste byproduct W is considered toxic, inducing a death rate proportional to its concentration:

Fig 2
Diagram of a simple metabolic network.

The parameters were set as follows. The stoichiometric coefficients NF = 2, NR = 38 are the characteristic ATP yields of glycolysis and respiration, respectively [61]. Maintenance demand is modeled as a constant drain of ATP at a rate e = 1.0625 mmol/gDW/h, typical of mammalian cells [39]. The maximum respiratory capacity is computed as rmax = Fthr × Vol × DW = 0.45 mmol/g/h, where Fthr = 0.9 mM/min is a glucose uptake threshold (per cytoplasmic volume) beyond which mammalian cells secrete lactate [58], Vol = 3 pL and DW = 0.9 ng are the volume [62] and dry weight, respectively, of mammalian (HeLa) cells, the later estimated from the dry mass fraction (≈ 30%, [63, BNID100387]) and total weight (= 3 ng [64]) of one HeLa cell. The concentration of substrate in the medium was set c = 15 mM, which is a typical glucose concentration in mammalian cell culture media (for example, RPMI-1640 [65]). Next, V = 0.5 mmol/gDW/h, also measured for HeLa cells [66] (the measured flux is per protein weight, so we multiplied by 0.5 to obtain a flux per cell dry weight, since roughly half of a generic cell dry weight is protein [63, BNID101955]). The parameter y = 348 mmol/gDW was adjusted so that the maximum growth rate was ≈1 day−1, which is within the range of duplication rates in mammalian cells [67, 68]. Finally, the toxicity of waste was set as τ = 0.0022 h−1mM−1, obtained from linearizing the death rate dependence on lactate in a mammalian cell culture reported by S. Dhir et al. [68]. To convert between grams of dry weight and number of cells, we used the cellular dry weight estimated above for HeLa cells (0.9ng).

Genome-scale metabolic network of CHO-K1

Exploiting the increasingly available information about cellular metabolism at the stoichiometric level [15, 16], a metabolic network can be reconstructed containing the biochemical reactions occurring inside a cell of interest. These reconstructions typically contain data about stoichiometric coefficients (Nik, cf. Eq 3), thermodynamic bounds (lbk, ubk, Li, cf. Eqs 5 and 6), and a biomass synthesis pseudo-reaction (yi, cf. Eq 3) [18, 41, 69]. Motivated by the fact that most therapeutic proteins requiring complex post-transnational modifications in the biotechnological industry are produced in Chinese hamster ovary (CHO) cell lines [48], we analyzed the steady states of a genome scale model of the CHO-K1 line [70]. Based on the latest consensus reconstruction of CHO metabolism available at the time of writing, containing 1766 genes and 6663 reactions, a cell-line-specific model for CHO-K1 was built by Hefzi et al. [70], comprising 4723 reactions (including exchanges) and 2773 metabolites (with cellular compartmentalization). It accounts for biomass synthesis through a virtual reaction that contains the moles of each metabolite required to synthesize one gram of biomass. The network recapitulates experimental growth rates and cell-line-specific amino acid auxotrophies.

In order to enforce Eq 4, we complemented this network with a set of reaction costs. Following T. Shlomi et. al [22], we assigned costs as follows: αk±=MWk±/kcat,k±, where MWk± and kcat,k± are the molecular weight and catalytic rate of the enzyme catalyzing reaction k in the given direction. The parameters MWk±, kcat,k± were gathered by T. Shlomi et. al from public repositories of enzymatic data. Missing values are set to the median of available values. An estimate of the enzyme mass fraction C = 0.078mg/mgDW was obtained for mammalian cells by the same authors. A constant maintenance energetic demand (cf. term ei in Eq 3) was added in the form of an ATP hydrolysis drain at a flux rate 2.24868mmol/gDW/h [39] (the reported value is for mouse LS cells, which we converted by accounting for the dry weight of CHO cells [70]). The maximum uptake rate of glucose was set at Vglc = 0.5 mmol/gDW/h, from previous models of cultured CHO cells fitted to experimental data [71, 72] (which also closely matches the values obtained from kinetic measurements on other mammalian cell lines [66]). However, kinetic parameters needed to estimate Vi for most metabolites are not known at present. Based on data in the literature [33, 68, 73], we estimated that the uptake rates of amino acids is typically one order of magnitude slower than the uptake rate of glucose, accordingly we set Vi = Vglc/10 for amino acids. Other metabolites have an unbounded uptake (Vi = ∞). In the simulations we used Iscove’s modified Dulbecco’s medium (IMDM), and set infinite concentrations for water, protons and oxygen. We converted between grams of dry weight and cell number using a cellular dry weight of 350pg [70].

The two toxic byproducts most commonly studied in mammalian cell cultures are ammonia and lactate. Their toxicity is primarily attributed to their effects in osmolarity and pH [3336]. It has been suggested that the accumulated toxicity may result in increased maintenance demands [55, 74] and in reduced biomass yields [55]. Parameters describing these effects quantitatively vary over an order of magnitude [57, 75, 76] depending on culture conditions and cell-line. In our model we incorporate these effects through the factor K and for the sake of specificity in this example we use:

K = (1+snh4/Knh4)−1(1+slac/Klac)−1

with Knh4 = 1.05mM, Klac = 8mM [67], and set μ = K × z.

Additional details

Numerical simulations were carried out in Julia [77]. Linear programs were solved with Gurobi [78]. The CHO-K1 metabolic network [70] was read and setup with all relevant parameters using a script written in Python with the COBRApy package [7981]. All scripts (which also include parameter values) are freely available in a public Github repository [82].

Results and discussion

General properties of steady states

In this section we present the general procedure to determine the steady states of Eqs 1 and 2 and discuss some general results of our model that are independent of the specificities of the metabolic networks of interest. The first step is to set the time-dependence in Eqs 1 and 2 to zero,



Note that Eq 13 depends on X and D only through the ratio 1/ξ = D/X (known in the literature as cell-specific perfusion rate, or CSPR [83]), such that ξ is the number of cells sustained in the culture per unit of medium supplied per unit time (the units of ξ are cells × time / volume). If we recall that in our cellular model, ui is constrained by a term that also depends on X and D only through ξ (cf. Eq 6), it immediately follows that the values of the uptakes and metabolite concentrations in steady state must be functions of ξ, which we denote by ui*(ξ) and si*(ξ) respectively. To compute ui*(ξ), solve the linear program of maximizing the biomass synthesis rate (z) subject to Eqs 35, but replacing Eq 6 with:

Li ≤ ui ≤ min {Vici/ξ}.

The resulting optimal value of z will be denoted by z*(ξ). Moreover, once ui*(ξ) is known, the stationary metabolite concentrations in the culture follow from Eq 13:


Then, given z*(ξ), the effective growth rate in steady state can also be given as a function of ξ, μ*(ξ), by evaluating K or σ using the concentrations si*(ξ) from Eq 15. Next, Eq 12 implies that the dilution rate at which a steady state occurs must also be a function of ξ, that we denote by D*(ξ). Combining this with the relation ξ = X/D, we obtain the steady state cell density, X*(ξ), as well:

D*(ξ) = μ*(ξ)/ϕX*(ξ) = ξμ*(ξ)/ϕ.

Note that while Eqs 12 and 13 are satisfied by any D ≥ 0 when X = 0, the value Dmax = D*(0) given by Eq 16 is actually the washout dilution rate, i.e., the minimum dilution rate that washes the culture of cells. Clearly all steady states with non-zero cell density are required to satisfy D*(ξ) < Dmax.

From Eq 16 it is evident that the function μ*(ξ) encodes all the information needed to get the values of X in steady state at different dilution rates and for any value of the bleeding coefficient ϕ. On the other hand, if multiple values of ξ are consistent with the same dilution rate (i.e. if the function D*(ξ) is not one-to-one), the system is multi-stable (i.e., multiple steady states coexist under identical external conditions). A necessary condition multi-stability is that μ*(ξ) is not monotonously decreasing. Since the biomass production rate z*(ξ) is a non-increasing function of ξ (proved in the Appendix), a change in the monotonicity of μ*(ξ) must be a consequence of toxic byproduct accumulation, modeled through the terms K and σ.

A noteworthy consequence of Eq 16 is that a plot displaying the parametric curve (ϕD*(ξ), ϕX*(ξ)) as a function of ξ is invariant to changes in ϕ. This means that for a given cell line and medium formulation, this curve can be obtained from measurements in a chemostat (which corresponds to ϕ = 1), and the result will also apply to any perfusion system with an arbitrary value of ϕ. Moreover, since si*(ξ) and ui*(ξ) are independent of ϕ, cellular metabolism in steady states is equivalent in the chemostat and any perfusion system (with an arbitrary value of ϕ), provided that the values of ξ = X/D in both systems match.

Finally, we mention that generally a threshold value ξm exists, such that a steady state with ξ = X/D is feasible only if ξξm. When ξ > ξm, some of the constrains in Eqs 3, 5 and 6 cannot be met. In degenerate scenarios we could have ξm = ∞ (e.g., this could be the case if the maintenance demand in Eq 3 is neglected) or ξm ≤ 0 (e.g., if the medium is so poor that the maintenance demand cannot be met even with a vanishingly small cell density). The parameter ξm arising in this way in our model, coincides with the definition of medium depth given in the literature [84], and it quantifies for a given medium composition the maximum cell density attainable per unit of medium supplied per unit time.


To determine the stability of steady states we compute the Jacobian eigenvalues of the system of Eqs 1 and 2. If the real parts are all negative the state is stable, but if at least one eigenvalue has a positive real part, the state is unstable [85]. The critical case where all eigenvalues have non-negative real parts but at least one of them has a zero real part is dealt with using the Center Manifold Theorem [86, § 8.1]. The Appendix contains a detailed mathematical treatment. Briefly, a steady state is unstable if μ*(ξ) is increasing in a neighborhood, and stable otherwise. We stated above that steady states of a given cell line in continuous culture, using a fixed medium formulation, can be given as functions of ξ. The condition for stability stated here is also uniquely determined by ξ and, in particular, it is independent of ϕ. Therefore, a steady state in a chemostat (with ϕ = 1) is stable if and only if the same steady state in perfusion (with a matching value of ξ, but arbitrary ϕ) is also stable. Since our results are qualitatively invariant to changes in ϕ, we set ϕ = 1 in what follows.

Insight from the toy model

We first consider the small metabolic network depicted in Fig 2. In this example, maximization of growth sets the nutrient uptake (u) and respiratory flux (r) at the maximum rates allowed by their respective upper bounds, in Eqs 6 and 7. Employing Eq 8 to determine the waste secretion rate (v) from u, r, we obtain:

u = min {Vc/ξ},  r = min {urmax},  vru.

Thus the toy model admits simple analytical expressions giving the rates of metabolic fluxes in steady states as functions of ξ. A minimum nutrient uptake rate um is required to sustain the maintenance demand e. Since most cell types are able to grow under certain conditions without waste secretion, we make the biologically reasonable assumption that umrmax (which is satisfied by the parameters chosen in Materials and methods). It then follows that um = e/(NF + NR). There are three critical thresholds in ξ that correspond to important qualitative changes in the culture:


These transitions can be interpreted in the following way: 1) if ξξm the growth rate is zero because the maintenance demand cannot be met; 2) if ξsecξξm, cells grow without secreting waste; 3) if ξξsec, there is waste secretion; 4) for ξξ0 cells are competing for the substrate and growth is limited by nutrient availability (cf. discussion before Eq 6); 5) finally, if ξξ0 there is nutrient excess and cells are growing at the maximum rate allowed by intrinsic kinetic limitations. We emphasize that the threshold ξsec carries a special metabolic significance, because it controls the switch between two qualitatively distinct metabolic modes: if ξξsec, respiration is saturated and the intermediate P overflows in the form of secreted waste, with a lower energy yield; on the other hand, if ξξsec, the cell relies entirely on respiration to generate energy, with a higher yield (cf. Fig 3).

Fig 3
Metabolic modes of toy model.

The medium carries a concentration c of primary nutrient and zero waste content. Under these assumptions, Eq 15 has the following analytical solution for the steady state values of the metabolite concentrations, s*(ξ), w*(ξ):

s*(ξ) = c − min{Vξc}

w*(ξ) = max{0, c − s*(ξ)−rmaxξ}

Note that s*(ξ) is a decreasing function of ξ, while w*(ξ) has at most a single maximum. Eqs 810, 17, 19 and 20 can be used to define μ*(ξ). Then D*(ξ), X*(ξ) are given by Eq 16.

Fig 4 shows plots of μ*(ξ), X*(ξ), s*(ξ) and w*(ξ) for this model. Parameter values are given in Materials and Methods. As ξ ranges from ξ = 0 to ξ = ξm, stable and unstable steady states are drawn in continuous and discontinuous line, respectively. The system is stable in two regimes: ξ [less, similar] 1 × 106 cells · day/mL, with high toxicity, low biomass yield and low cell density; and ξξsec, with no toxicity, high biomass yield and high cell density that decays as ξ increases. The later states rely solely on respiration for energy generation (Fig 3a), while the former states exhibit overflow metabolism (Fig 3b). Waste concentration initially increases with ξ until a maximum value is reached. Then w decays during the unstable phase, all the way to zero at ξsec, where waste secretion stops and the system becomes stable again.

Fig 4
Steady states of toy model.

Intuitively, unstable states become stressed due to high levels of toxicity, which also makes the system very sensitive to perturbations. The typical behavior of nutrients and waste products (in particular glucose and lactate, respectively) in continuous cell cultures, as observed in experiments [84], is that as ξ increases, nutrient concentration in the culture decreases while waste initially accumulates [84] but eventually phases out as cells switch towards higher-yield metabolic pathways [10, 11]. This behavior is qualitatively reproduced by s and w in our toy model.

The function μ*(ξ) is not monotonically decreasing. As explained above, this implies a coexistence of multiple steady states under identical external conditions. This is readily apparent in a bifurcation diagram of the steady values of X versus D, as shown in Fig 5a. In a range of dilution rates (0.25 [less, similar] D [less, similar] 0.7, units: day−1), the system exhibits three steady states, one of which is unstable (discontinuous line in the figure), while the other two are stable (continuous line in the figure). Thus a stable state of high-cell density coexists with another of low cell density, over a range of dilution rates. Cellular metabolism in the former state is respiratory (Fig 3a), whereas cells in the later state exhibit an overflow metabolism (Fig 3b). The unstable state is an intermediate transition state lying between these two extremes.

Fig 5
Cell density versus dilution rate in the solvable model.

Multi-stability of continuous cultures has been repeatedly observed in experiments [1014]. In our model it is a direct consequence of toxicity induced by the accumulation of waste [87]. Small variations in the dilution rate near D ≈ 0.25 day−1 or D ≈ 0.7 day−1 result in discontinuous transitions where the cell density rises or drops abruptly, respectively. These jumps can be traced around an hysteresis loop, drawn in orange arrows in Fig 5a. More generally one may also expect that the system jumps from one state to the other due to random fluctuations. In particular, note that the basin of attraction of the high-cell density state decreases with D (since the discontinuous line of unstable states eventually intersects with the high cell density states). Therefore, our toy model exhibits a plausible mechanism through which increasing dilution rates translate into high cell density states with diminishing resilience to perturbations.

The role of toxicity becomes evident if we consider ideal cells resistant to waste accumulation (setting τ = 0). The resulting plot of X vs. D in this case (Fig 5b) reveals a single stable steady state for each value of the dilution rate. There is a discontinuous transition at the washout dilution rate (Dmax), where the cell density suddenly drops to zero. Away from this value, the system is resilient to perturbations since there are no multiple steady states between which jumps can occur.

Multi-stability implies that system dynamics are non-trivial, in the sense that different trajectories might lead to different steady states. Therefore it is important for industrial applications to understand how the system is driven to one or another state. We numerically solved Eqs 1 and 2 by performing the FBA optimization at each instant of time, in a manner analogous to DFBA [27]. With the parameter values given in Materials and Methods, we simulated the response of the system to three different profiles of the dilution in time. First, in Fig 6a, a constant dilution rate of D = 0.6 day−1 is used. Two possible stable steady states are consistent with this dilution rate, attaining different cell densities (cf. Fig 5a). Starting from a very low initial cell density, the system responds by settling at the steady state of lowest cell-density. As can be appreciated in the bottom row of Fig 6a, this state is characterized by an accumulation of toxic waste that prevents further cell growth. A smarter manipulation of the dilution rate makes the high-cell density state accessible. This is demonstrated in Fig 6b, where D starts from a lower value (D = 0.2 day−1), and is gradually increased until the final value (D = 0.6 day−1) is reached. The final cell density resulting from this smooth increase of the dilution rate is five-fold larger than the one obtained with a constant dilution rate. This state is also characterized by very low levels of waste accumulation (cf. last row of Fig 6b). We stress that external conditions in the final steady state (dilution rate and medium formulation) are the same in both cases.

Fig 6
Dynamic simulations of toy model.

Finally, Fig 6c shows how the dilution rate can be manipulated to switch from one steady state to another. Starting from the final state of the simulation in Fig 6a, the dilution rate is first decreased to a low value (D = 0.2 day−1), and then it is pushed back up to the starting value (D = 0.6 day−1). The system responds by switching from the state with low cell-density to the state with high-cell density. These simulations nicely reproduce the qualitative features of the experiment performed by B. Follstad et al. [11], where a continuous cell culture under the same steady external conditions (dilution rate and medium) switches between different steady states by transient manipulations of the dilution rate. The response of the cell density to transient manipulations of the dilution rate best illustrated in the X, D plane (cf. first row of Fig 6). Then it becomes obvious that the dilution rate must be pushed down to ≈ 0.2 day−1, otherwise the system will not leave the low cell density state.

Analysis of the CHO-K1 cell-line using a genome-scale metabolic network

We determined the steady states of a continuous cell culture of the CHO-K1 line. Cellular metabolism was modeled using the reconstruction given by Hefzi et al. [70], the most complete available at the time of writing. In the simulations we used Iscove’s modified Dulbecco’s growth media (IMDM), which is typically employed in mammalian systems. Similar to what we found in the toy model (cf. Fig 3), and in qualitative agreement with experimental observations [10, 14], cells exhibited several metabolic transitions between distinct flux modes as ξ was varied. However, in contrast to the the toy model, the CHO-K1 genome-scale metabolic network displays a rich multitude of transitions, as expected from its greater complexity. Because of their importance in the performance of the culture, we focused on metabolic changes that have an impact on macroscopic properties of the bioreactor, i.e., those that affect metabolite exchanges with the extracellular media (ui).

Although many classifications are possible, we organized our discussion by focusing on five qualitatively different modes based on the secretion of lactate and formate. Fig 7 shows cartoon diagrams of these phases in order of increasing ξ. On the top of each diagram we annotate the nutrients that became limiting for growth during a phase. Blue arrows indicate consumption and red arrows secretion. We focused on metabolites that changed their role between phases. In particular, NH4 was secreted in all phases and therefore was omitted from the figure to reduce clutter. A more detailed representation of our results is given in Fig 8, which shows metabolite concentrations (si) and uptakes (ui) in steady states as functions of ξ for a sub-set of selected metabolites. Red lines in these plots indicate the transitions depicted in Fig 7.

Fig 7
CHO-K1 metabolic exchange modes.
Fig 8
Uptake rates and concentrations of selected extracellular metabolites as functions of ξ.

For small values of ξ we found that glucose and almost all the amino acids available in the media were consumed, but none of them reached limiting concentrations. We call this the nutrient excess phase, where substrate uptake is limited only by intrinsic kinetic properties of cellular transporters. Remarkably lactate was not secreted at this stage, since pyruvate was converted instead to alanine [33] (although a small fraction of pyruvate was secreted as well [88]). The cell also produced succinate [89] and formate, the later being an overflow product of one-carbon metabolism of serine and glycine [53].

As ξ continues to increase, the first metabolite that becomes limiting is serine. This marks the end of the nutrient excess phase, coinciding also with the onset of lactate secretion. At this point pyruvate is no longer secreted into the culture. Remarkably, aspartate switches from being a secreted byproduct in the first phase [89], to consumption. Even more striking is that the specific uptake rate of aspartate and proline quickly increase until both reach limiting concentrations. A third phase is entered when succinate and formate production ceases, coinciding with a limitation of glycine. Histidine consumption rises steeply until it too reaches limiting concentrations. Other nutrients that limit growth include tyrosine, tryptophan, arginine, lysine and phenylalanine. This phase is also characterized by secretion of acetaldehyde. Remarkably, formate secretion is resumed in a later phase, where glucose, glutamine and asparagine also become limiting.

Finally, as ξ approaches the maximum value ξm, lactate and alanine secretion cease. This ideal state attains the highest possible biomass yield per unit of medium supplied per unit time. Note that the increase of ξ has brought an overall qualitative switch to a state of metabolic efficiency where the number of secreted byproducts has dropped significantly, compared to the nutrient excess phase. Notably, the cell-specific ammonia secretion was sustained even in the states of highest biomass yield, indicating a nitrogen imbalance. This behavior has been seen qualitatively in some experiments. For example, using a CHO-derived cell line [12], secretion of ammonia was sustained even after a transition to an efficient metabolic phenotype with low lactate secretion and high cellular yields. However this observation seems to be cell-line dependent, and in another experiment with an hybridoma, ammonia accumulation decreased with increasing ξ [84].

All of the secreted metabolites predicted by our model have been verified in experiments in mammalian cell cultures [33, 53, 88, 89], with the exception of acetaldehyde. Although some of these experiments do not match the cell line and media used in our simulations (CHO and IMDM), these byproducts have been observed in mammalian cell cultures in a variety of conditions, suggesting that they are not restricted to a specific cell line or media. For acetaldehyde our search in the literature did not reveal any experimental evidence refuting or validating its presence in mammalian cell cultures. A possible explanation is the high reactivity of this metabolite. Acetaldehyde binds covalently to glutathione and proteins, forming adducts that are subsequently detoxified [96, 97]. However, since the metabolic network does not account for these interactions [70], the model predicts excretion of pure acetaldehyde instead. We thank an anonymous reviewer for bringing this fact to our attention.

The performance of cell-lines and media are typically evaluated by measurements performed in batch experiments [90]. Measurements performed in the exponential phase of batch only reveal the behavior of continuous cultures at very low ξ, in conditions of nutrient excess. The existence of a rich multitude of qualitatively distinct metabolic behaviors at higher values of ξ is missed by these experiments and therefore the assessment should not be extrapolated to high-density perfusion systems. As our analysis reveals, several nutrients may switch from basal rates of consumption to growth limiting at later values of ξ, while others go from secreted byproducts to consumption [91]. These examples indicate that nutrients could be in excess in a batch experiment but need not be so in the ideal regime of high-cell density perfusion cultures, at high ξ. Our model suggests that a better characterization of a cell-line and media formulation can be obtained in a chemostat, since the full spectrum of values of ξ can be explored in this device and it faithfully reproduces all the metabolic transitions found in perfusion.

The effects of toxic byproduct accumulation are explored in Fig 9. We considered the toxic effects of the most commonly studied metabolites in this regard: lactate and ammonia, although the model can easily accommodate the effects of additional toxic compounds if necessary. In Fig 9a we plot the effective growth rate, μ as function of ξ. Stable states are drawn in continuous line, unstable states are dashed and the red dots indicate the metabolic transitions depicted in Fig 7. Note that μ*(ξ) is not monotonous. In particular, metabolic transitions resulting in lactate and ammonia secretion peaks produce a sink in the curve μ*(ξ). On the other hand, metabolic transitions associated to the secretion of other non-toxic byproducts do not imply changes in the monotonicity of μ*(ξ).

Fig 9
Steady states of CHO-K1.

The non-monotonicity of μ*(ξ) results in multiple stable states coexisting at the same dilution rate, as evident in the bifurcation diagram Fig 9b. This resonates with the results obtained in the simpler model considered above, and is also consistent with many experimental observations of bi-stability in the literature [1014]. The regime with high-cell density corresponds to a higher value of ξ and exhibits a lower accumulation of toxic byproducts (lactate and ammonia). Metabolism in this regime is also more efficient, with less byproducts secreted (cf. Fig 7). On the other hand, low cell density states are wasteful, with high levels of environmental toxicity preventing further cell growth. Again, bi-stability implies the existence of an hysteresis loop (orange arrows in the figure), where the system may suffer abrupt transitions between high and low cell densities.

Concluding remarks

In this work we have presented a model of cellular metabolism in continuous cell culture. Although similar in spirit to DFBA, our dynamic equations include terms accounting for the continuous medium exchange that enables steady states in this system. We presented a simple method to compute the steady states of the culture as a function of the ratio between cell density and dilution rate (ξ = X/D), scalable to metabolic networks of arbitrary complexity. In the literature 1/ξ is known as the cell-specific perfusion rate (CSPR), introduced by S. Ozturk [83] who already made the empirical observation that control of the CSPR can be used to maintain a constant cell environment independent of cell growth [83]. Our model theoretically supports this idea and leads to a stronger conclusion: that for a given cell line and medium formulation, the steady state values of the macroscopic variables of the bioreactor are all unequivocally determined by ξ. Therefore, ξ is an ideal control parameter to fix a desired steady state in a continuous cell culture.

The model is consistent with multi-stability, a phenomenon repeatedly observed in experiments in continuous cell cultures where multiple steady states coexist under identical external conditions. Moreover, our model accounts for metabolic switches between flux modes, experimentally observed in continuous cell culture in response to variations in the dilution rate [92]. These transitions affect the consumption or secretion of metabolites and the set of nutrients limiting growth. As a consequence, the metabolic landscape of steady states in perfusion cell cultures is complex and cannot be reproduced in batch cultures. This has the practical implication that assesments of medium quality and cell line performance carried out in batch [90] should not be extrapolated to perfusion, since they might be missleading in this setting.

However, our analysis reveals a simple scaling law between steady states in the chemostat and any perfusion system. The landscape of metabolic transitions in the later system can be faithfully reproduced in the chemostat. Thus, for a fixed cell-line and medium formulation, the diagram displaying the values of ϕX versus ϕD in steady state is invariant across perfusion systems with any bleeding ratio (ϕ), cf. Eq 16, while metabolism is equivalent if the ratio ξ = X/D is the same. The practical consequence is that the chemostat is an ideal experimental model where cell-lines and medium formulations can be benchmarked for their performance in high-cell density industrial continuous cultures.

Further, the model predicts that multi-stability is a consequence of negative feedback on cell growth due to accumulation of toxic byproducts in the culture. The qualitative complexity of the ϕX versus ϕD diagram depends only on the behavior of toxic metabolites. Moreover, multi-stability implies that the system is sensitive to initial conditions and transient manipulations of external parameters. In practice, the dilution rate must be manipulated carefully to bring the system to a desired state. Thus, starting from a seed of low cell density, sharp increases of the dilution may land the system on a steady state of high toxicity and low biomass. On the other hand, slowly increasing the dilution rate will surely lead towards high-cell density states.

The conclusions stated above rely on the validity of our assumptions. In particular, we have considered a homogeneous cell population in a well-mixed bioreactor. Both assumptions are behind many models published in the field and provide reasonable fits to experimental data [37]. Mechanical stirring of the culture typically achieves a well-mixed solution, but care must be taken to prevent mechanical damage to the cells [93] (but see Ref [94]). Moreover, that the cell population can be treated attending only to its average properties is justified by the large number of cells in a typical culture (~106 − 108 cells/mL), although in some settings cell-to-cell heterogeneity might become relevant [95]. Next, to develop a specific model of cellular metabolism, we adopted a flux-balance approach [38], where cells are assumed to optimize their metabolism towards growth rate maximization. Although this framework is well supported in the literature [20], it is worth noting that we did not consider the kinetics of intracellular metabolites or additional regulatory mechanisms that may also control metabolic fluxes. Additionally, the quantitative predictions of the model rely on the accuracy of parameters found in the literature and databases. Among these, the flux cost coefficients (αi, Eq 4) are not available for many enzymes. If too many of these parameters are absent, calculations from FBA might be degenerate [17, 54]. Another important omission from the present model is that we did not consider explicitly the exchange between the culture and a gaseous phase. In particular, this includes oxygen exchange. Therefore our approach is only valid if this exchange does not become limiting to cellular growth. Despite these limitations, we have shown that the model predictions are in qualitative agreement with experimental data. More importantly, the conclusions stated above are independent of the values of model parameters.

Finally, we discuss briefly the relevance of our work in tissues. Constant perfusion flows resulting from the blood stream imply that tissues can reach a steady state in principle similar to the chemostat. Although many complications must be considered in this case, including circadian oscillations, heterogeneous cell populations and spatial structure (Ref. [60] presents an attempt to deal with the later two), some of the conclusions derived in this manuscript could still be relevant in this context. For example, it would be very interesting to explore if the metabolic profile of cells from different tissues can be correlated to the ratio between cell density and average perfusion rate experienced at the tissue site (i.e, the ξ parameter defined in this work). Moreover, the majority of in vitro experiments studying the metabolism of human cells, are conducted in batch due to ease of operation. Our work points out that the metabolic profile observed in this setting might be very different from the state in a perfused tissue.

Supporting information

S1 Appendix

Supplementary note.

Contains mathematical derivations.


S1 Data

Source code.

Contains the source code of scripts used in the manuscript.



The authors warmly thank Tamy Boggiano and Ernesto Chico for many helpful discussions and for reviewing this manuscript.

Funding Statement

This project has received funding from the European Union’s Horizon 2020 research and innovation programme MSCA-RISE-2016 under grant agreement No. 734439 INFERNET. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

Data Availability

All relevant data are within the paper and its Supporting Information files.


1. Werner RG, Walz F, Noé W, Konrad A. Safety and economic aspects of continuous mammalian cell culture. Journal of Biotechnology. 1992;22(1-2):51–68. doi: 10.1016/0168-1656(92)90132-S [PubMed]
2. Griffiths JB. Animal cell culture processes—batch or continuous? Journal of Biotechnology. 1992;22(1-2):21–30. doi: 10.1016/0168-1656(92)90129-W [PubMed]
3. Kadouri A, Spier RE. Some myths and messages concerning the batch and continuous culture of animal cells. Cytotechnology. 1997;24(2):89–98. doi: 10.1023/A:1007932614011 [PMC free article] [PubMed]
4. Werner RG, Noe W. Letter to the Editor. Cytotechnology. 1998;26(2):81–82. doi: 10.1023/A:1007985828899 [PMC free article] [PubMed]
5. Croughan MS, Konstantinov KB, Cooney C. The future of industrial bioprocessing: Batch or continuous? Biotechnology and Bioengineering. 2015;112(4):648–651. doi: 10.1002/bit.25529 [PubMed]
6. Konstantinov KB, Cooney CL. White Paper on Continuous Bioprocessing May 20–21 2014 Continuous Manufacturing Symposium. Journal of Pharmaceutical Sciences. 2015;104(3):813–820. doi: 10.1002/jps.24268 [PubMed]
7. Novick A, Szilard L. Description of the chemostat. Science. 1950;112:716 doi: 10.1126/science.112.2920.715 [PubMed]
8. Monod J. La technique de culture continue, théorie et applications. Ann Inst Pasteur. 1950;79:390–410.
9. Castilho L, Medronho R. Cell Retention Devices for Suspended-Cell Perfusion Cultures. Tools and Applications of Biochemical Engineering Science. 2002;74:129–169. doi: 10.1007/3-540-45736-4_7 [PubMed]
10. Europa AF, Gambhir A, Fu PC, Hu WS. Multiple steady states with distinct cellular metabolism in continuous culture of mammalian cells. Biotechnology and Bioengineering. 2000;67(1):25–34. doi: 10.1002/(SICI)1097-0290(20000105)67:1%3C25::AID-BIT4%3E3.0.CO;2-K [PubMed]
11. Follstad BD, Balcarcel RR, Stephanopoulos G, Wang DIC. Metabolic flux analysis of hybridoma continuous culture steady state multiplicity. Biotechnology and Bioengineering. 1999;63(6):675–683. doi: 10.1002/(SICI)1097-0290(19990620)63:6%3C675::AID-BIT5%3E3.0.CO;2-R [PubMed]
12. Altamirano C, Illanes A, Casablancas A, Gámez X, Cairó JJ, Gòdia C. Analysis of CHO cells metabolic redistribution in a glutamate-based defined medium in continuous culture. Biotechnology Progress. 2001;17(6):1032–1041. doi: 10.1021/bp0100981 [PubMed]
13. Hayter PM, Curling EMa, Baines AJ, Jenkins N, Salmon I, Strange PG, et al. Glucose-Limited Chemostat Culture of Chinese Hamster Ovary Cells Producing Recombinant Human Interferon-γ. Biotechnology and Bioengineering. 1992;39:327–335. doi: 10.1002/bit.260390311 [PubMed]
14. Gambhir A, Korke R, Lee J, Fu PC, Europa AF, Hu WS. Analysis of cellular metabolism of hybridoma cells at distinct physiological states. Journal of bioscience and bioengineering. 2003;95(4):317–27. doi: 10.1016/S1389-1723(03)80062-2 [PubMed]
15. Kanehisa M, Goto S, Sato Y, Kawashima M, Furumichi M, Tanabe M. Data, information, knowledge and principle: Back to metabolism in KEGG. Nucleic Acids Research. 2014;42(D1):199–205. doi: 10.1093/nar/gkt1076 [PMC free article] [PubMed]
16. Caspi R, Altman T, Billington R, Dreher K, Foerster H, Fulcher CA, et al. The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of Pathway/Genome Databases. Nucleic Acids Research. 2014;42(D1):D459–D471. doi: 10.1093/nar/gkt1103 [PMC free article] [PubMed]
17. Noor E, Flamholz A, Bar-Even A, Davidi D, Milo R, Liebermeister W. The Protein Cost of Metabolic Fluxes: Prediction from Enzymatic Rate Laws and Cost Minimization. PLoS Computational Biology. 2016;12(11). doi: 10.1371/journal.pcbi.1005167 [PMC free article] [PubMed]
18. Lewis NE, Nagarajan H, Palsson BØ. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nature Reviews Microbiology. 2012;10(4):291–305. doi: 10.1038/nrmicro2737 [PMC free article] [PubMed]
19. Ibarra RU, Edwards JS, Palsson BØ. Escherichia coli K-12 undergoes adaptive evolution to achieve in silico predicted optimal growth. Nature. 2002;420(6912):186–189. doi: 10.1038/nature01149 [PubMed]
20. Palsson BØ. Systems Biology: Properties of Reconstructed Networks. Cambridge University Press; 2006.
21. Edwards JS, Ibarra RU, Palsson BØ. In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. Nature biotechnology. 2001;19(2):125–30. doi: 10.1038/84379 [PubMed]
22. Shlomi T, Benyamini T, Gottlieb E, Sharan R, Ruppin E. Genome-Scale Metabolic Modeling Elucidates the Role of Proliferative Adaptation in Causing the Warburg Effect. PLoS Computational Biology. 2011;7(3):e1002018 doi: 10.1371/journal.pcbi.1002018 [PMC free article] [PubMed]
23. Jouhten P, Wiebe M, Penttilä M. Dynamic flux balance analysis of the metabolism of Saccharomyces cerevisiae during the shift from fully respirative or respirofermentative metabolic states to anaerobiosis. FEBS Journal. 2012;279(18):3338–3354. doi: 10.1111/j.1742-4658.2012.08649.x [PubMed]
24. Braunstein A, Muntoni AP, Pagnani A. An analytic approximation of the feasible space of metabolic networks. Nature Communications. 2017;8:14915 doi: 10.1038/ncomms14915 [PMC free article] [PubMed]
25. Braunstein A, Mulet R, Pagnani A. Estimating the size of the solution space of metabolic networks. BMC Bioinformatics. 2008;9:240 doi: 10.1186/1471-2105-9-240 [PMC free article] [PubMed]
26. Fernández-de Cossio-Díaz J, Mulet R. Fast inference of ill-posed problems within a convex space. Journal of Statistical Mechanics: Theory and Experiment. 2016;(7):073207.
27. Mahadevan R, Edwards JS, Doyle FJ. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophysical Journal. 2002;83(3):1331–40. doi: 10.1016/S0006-3495(02)73903-9 [PubMed]
28. Meadows AL, Karnik R, Lam H, Forestell S, Snedecor B. Application of dynamic flux balance analysis to an industrial Escherichia coli fermentation. Metabolic Engineering. 2010;12(2):150–160. doi: 10.1016/j.ymben.2009.07.006 [PubMed]
29. Visek WJ. Some aspects of ammonia toxicity in animal cells. Journal of dairy science. 1968;51(2):286–295. doi: 10.3168/jds.S0022-0302(68)86976-0 [PubMed]
30. Hassell T, Gleave S, Butler M. Growth inhibition in animal cell culture. Applied Biochemistry and Biotechnology. 1991;30(1):29–41. doi: 10.1007/BF02922022 [PubMed]
31. Guthke R, Knorre WA. Bistability in a model of microbial product formation. Zeitschrift für allgemeine Mikrobiologie. 1980;20(7):441–447. doi: 10.1002/jobm.19800200703 [PubMed]
32. Hegewald E, Wolleschensky B, Guthke R, Neubert M, Knorre WA. Instabilities of product formation in a fed-batch culture of Penicillium chrysogenum. Biotechnology and Bioengineering. 1981;23(7):1563–1572. doi: 10.1002/bit.260230715
33. Ozturk SS, Riley MR, Palsson BØ. Effects of ammonia and lactate on hybridoma growth, metabolism, and antibody production. Biotechnology and Bioengineering. 1992;39(4):418–431. doi: 10.1002/bit.260390408 [PubMed]
34. Bakker WAM, Schäfer T, Beeftink HH, Tramper J, De Gooijer CD. Hybridomas in a bioreactor cascade: modeling and determination of growth and death kinetics. Cytotechnology. 1996;21(3):263–277. doi: 10.1007/BF00365349 [PubMed]
35. Hu WS, Aunins JG. Large-scale mammalian cell culture. Current Opinion in Biotechnology. 1997;8(2):148–153. doi: 10.1016/S0958-1669(97)80093-6 [PubMed]
36. Schneider M, Marison IW, von Stockar U. The importance of ammonia in mammalian cell culture. Journal of biotechnology. 1996;46(3):161–85. doi: 10.1016/0168-1656(95)00196-4 [PubMed]
37. Ben Yahia B, Malphettes L, Heinzle E. Macroscopic modeling of mammalian cell growth and metabolism. Applied Microbiology and Biotechnology. 2015;99(17):7009–7024. doi: 10.1007/s00253-015-6743-6 [PMC free article] [PubMed]
38. Edwards JS, Covert M, Palsson BØ. Metabolic modelling of microbes: the flux-balance approach. Environmental Microbiology. 2002;4(3):133–140. doi: 10.1046/j.1462-2920.2002.00282.x [PubMed]
39. Kilburn DG, Lilly MD, Webb FC. The energetics of mammalian cell growth. Journal of cell science. 1969;4(3):645–654. [PubMed]
40. Sheikh K, Förster J, Nielsen LK. Modeling Hybridoma Cell Metabolism Using a Generic Genome-Scale Metabolic Model of Mus musculus. Biotechnology Progress. 2008;21(1):112–121. doi: 10.1021/bp0498138 [PubMed]
41. Feist AM, Palsson BØ. The biomass objective function. Current Opinion in Microbiology. 2010;13(3):344–349. doi: 10.1016/j.mib.2010.03.003 [PMC free article] [PubMed]
42. Feist AM, Palsson BØ. What do cells actually want? Genome Biology. 2016;17(1):110 doi: 10.1186/s13059-016-0983-3 [PMC free article] [PubMed]
43. Beg QK, Vazquez A, Ernst J, de Menezes MA, Bar-Joseph Z, Barabási ALAL, et al. Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity. Proceedings of the National Academy of Sciences. 2007;104(31):12663–12668. doi: 10.1073/pnas.0609845104 [PubMed]
44. Vazquez A, Oltvai ZN. Macromolecular crowding explains overflow metabolism in cells. Scientific Reports. 2016;6:31007 doi: 10.1038/srep31007 [PMC free article] [PubMed]
45. Zhuang K, Vemuri GN, Mahadevan R. Economics of membrane occupancy and respiro-fermentation. Molecular Systems Biology. 2014;7(1):500–500. doi: 10.1038/msb.2011.34 [PMC free article] [PubMed]
46. Basan M, Hui S, Okano H, Zhang Z, Shen Y, Williamson JR, et al. Overflow metabolism in Escherichia coli results from efficient proteome allocation. Nature. 2015;528(7580):99–104. doi: 10.1038/nature15765 [PMC free article] [PubMed]
47. Beard DA, Qian H. Chemical Biophysics—Quantitative Analysis of Cellular Systems. Cambridge University Press; 2008.
48. Wurm FM. Production of recombinant protein therapeutics in cultivated mammalian cells. Nature Biotechnology. 2004;22(11):1393–1398. doi: 10.1038/nbt1026 [PubMed]
49. Vanderbei RJ. Linear Programming Foundations and Extensions. 4th ed Springer; 2014.
50. Varma A, Palsson BØ. Metabolic Flux Balancing: Basic Concepts, Scientific and Practical Use. Bio/Technology. 1994;12(10):994–998. doi: 10.1038/nbt1094-994
51. Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nature Publishing Group. 2010;28(3):245–248. [PMC free article] [PubMed]
52. Vazquez A. Optimal cytoplasmatic density and flux balance model under macromolecular crowding effects. Journal of Theoretical Biology. 2010;264(2):356–359. doi: 10.1016/j.jtbi.2010.02.024 [PubMed]
53. Meiser J, Tumanov S, Maddocks O, Labuschagne CF, Athineos D, Van Den Broek N, et al. Serine one-carbon catabolism with formate overflow. Science Advances. 2016;2(10):e1601273–e1601273. doi: 10.1126/sciadv.1601273 [PMC free article] [PubMed]
54. Mahadevan R, Schilling CH. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metabolic Engineering. 2003;5(4):264–276. doi: 10.1016/j.ymben.2003.09.002 [PubMed]
55. Martinelle K, Häggström L. Mechanisms of ammonia and ammonium ion toxicity in animal cells: Transport across cell membranes. Journal of Biotechnology. 1993;30(3):339–350. doi: 10.1016/0168-1656(93)90148-G [PubMed]
56. Bertolazzi E. A combination formula of Michaelis-Menten-Monod type. Computers and Mathematics with Applications. 2005;50(1-2):201–215. doi: 10.1016/j.camwa.2004.10.045
57. Pörtner R, Schäfer T. Modelling hybridoma cell growth and metabolism—a comparison of selected models and data. Journal of Biotechnology. 1996;49(1-3):119–135. doi: 10.1016/0168-1656(96)01535-0 [PubMed]
58. Vazquez A, Liu J, Zhou Y, Oltvai ZN. Catabolic efficiency of aerobic glycolysis: the Warburg effect revisited. BMC systems biology. 2010;4:58 doi: 10.1186/1752-0509-4-58 [PMC free article] [PubMed]
59. Capuani F, De Martino D, Marinari E, De Martino A. Quantitative constraint-based computational model of tumor-to-stroma coupling via lactate shuttle. Scientific Reports. 2015;5:11880 doi: 10.1038/srep11880 [PMC free article] [PubMed]
60. Fernandez-de Cossio-Diaz J, De Martino A, Mulet R. Microenvironmental cooperation promotes early spread and bistability of a Warburg-like phenotype. Scientific Reports. 2017;7:3103 doi: 10.1038/s41598-017-03342-3 [PMC free article] [PubMed]
61. Alberts B, Johnson A, Lewis Julian, Morgan D, Raff M, Roberts K, et al. Molecular Biology of the Cell. 6th ed; 2014.
62. Zhao L, Kroenke CD, Song J, Piwnica-Worms D, Ackerman JJH, Neil JJ. Intracellular water-specific MR of microbead-adherent cells: the HeLa cell intracellular water exchange lifetime. NMR in Biomedicine. 2008;21(2):159–164. doi: 10.1002/nbm.1173 [PMC free article] [PubMed]
63. Milo R, Jorgensen P, Moran U, Weber G, Springer M. BioNumbers—the database of key numbers in molecular and cell biology. Nucleic Acids Research. 2010;38(Database):D750–D753. doi: 10.1093/nar/gkp889 [PMC free article] [PubMed]
64. Park K, Jang J, Irimia D, Sturgis J, Lee J, Robinson JP, et al. ’Living cantilever arrays’ for characterization of mass of single live cells in fluids. Lab on a chip. 2008;8(7):1034–1041. doi: 10.1039/b803601b [PMC free article] [PubMed]
65. Moore GE, Gerner RE, Franklin HA. Culture of normal human leukocytes. Journal of the American Medical Association. 1967;199(8):519–524. doi: 10.1001/jama.1967.03120080053007 [PubMed]
66. Rodríguez-Enríquez S, Marín-Hernández A, Gallardo-Pérez JC, Moreno-Sánchez R. Kinetics of transport and phosphorylation of glucose in cancer cells. Journal of Cellular Physiology. 2009;221(3):552–559. doi: 10.1002/jcp.21885 [PubMed]
67. Bree MA, Dhurjati P, Geoghegan RF, Robnett B. Kinetic modelling of hybridoma cell growth and immunoglobulin production in a large-scale suspension culture. Biotechnology and Bioengineering. 1988;32(8):1067–1072. doi: 10.1002/bit.260320814 [PubMed]
68. Dhir S, Morrow KJ, Rhinehart RR, Wiesner T. Dynamic optimization of hybridoma growth in a fed-batch bioreactor. Biotechnology and Bioengineering. 2000;67(2):197–205. doi: 10.1002/(SICI)1097-0290(20000120)67:2%3C197::AID-BIT9%3E3.0.CO;2-W [PubMed]
69. Schellenberger J. Monte Carlo simulation in systems biology; 2010.
70. Hefzi H, Ang KS, Hanscho M, Bordbar A, Ruckerbauer DE, Lakshmanan M, et al. A Consensus Genome-scale Reconstruction of Chinese Hamster Ovary Cell Metabolism. Cell Systems. 2016;3(5):434–443.e8. doi: 10.1016/j.cels.2016.10.020 [PMC free article] [PubMed]
71. Nolan RP, Lee K. Dynamic model of CHO cell metabolism. Metabolic Engineering. 2011;13(1):108–124. doi: 10.1016/j.ymben.2010.09.003 [PubMed]
72. Kiparissides A, Koutinas M, Kontoravdi C, Mantalaris A, Pistikopoulos EN. ‘Closing the loop’ in biological systems modeling—From the in silico to the in vitro. Automatica. 2011;47(6):1147–1155. doi: 10.1016/j.automatica.2011.01.013
73. Altamirano C, Paredes C, Cairó JJ, Godia F. Improvement of CHO Cell Culture Medium Formulation: Simultaneous Substitution of Glucose and Glutamine. Biotechnology Progress. 2000;16(1):69–75. doi: 10.1021/bp990124j [PubMed]
74. Lao MS, Toth D. Effects of Ammonium and Lactate on Growth and Metabolism of a Recombinant Chinese Hamster Ovary Cell Culture. Biotechnology Progress. 1997;13(5):688–691. doi: 10.1021/bp9602360 [PubMed]
75. Glacken MW, Adema E, Sinskey AJ. Mathematical Descriptions of Hybidoma Culture Kinetics: I. Initial Metabolic Rates. Biotechnology and Bioengineering. 1988;32:491–506. doi: 10.1002/bit.260320412 [PubMed]
76. Gaertner JG, Dhurjati P. Fractional factorial study of hybridoma behavior. 1. Kinetics of growth and antibody production. Biotechnology Progress. 1993;9(3):298–308. doi: 10.1021/bp00021a009 [PubMed]
77. Bezanson J, Edelman A, Karpinski S, Shah VB. Julia: A Fresh Approach to Numerical Computing. SIAM Review. 2017;59(1):65–98. doi: 10.1137/141000671
78. Gurobi Optimization, Inc. Gurobi Optimizer Reference Manual; 2016. Available from:
79. Ebrahim A, Lerman JA, Palsson BØ, Hyduke DR. COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC systems biology. 2013;7(1):74 doi: 10.1186/1752-0509-7-74 [PMC free article] [PubMed]
80. Becker SA, Feist AM, Mo ML, Hannum G, Palsson BØ, Herrgård MJ. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox. Nature protocols. 2007;2(3):727–38. doi: 10.1038/nprot.2007.99 [PubMed]
81. Schellenberger J, Que R, Fleming RMT, Thiele I, Orth JD, Feist AM, et al. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox v2.0. Nature Protocols. 2011;6(9):1290–1307. doi: 10.1038/nprot.2011.308 [PMC free article] [PubMed]
82. Fernandez-de Cossio-Diaz J. Simulation scripts used in this work; 2017. Available from:
83. Ozturk SS. Engineering challenges in high density cell culture systems. Cytotechnology. 1996;22(1-3):3–16. doi: 10.1007/BF00353919 [PubMed]
84. Konstantinov KB, Goudar CT, Ng M, Meneses R, Thrift J, Chuppa S, et al. The “push-to-low” approach for optimization of high-density perfusion cultures of animal cells. Advances in Biochemical Engineering/Biotechnology. 2006;101(July):75–98. doi: 10.1007/10_016 [PubMed]
85. Strogatz SH. Nonlinear Dynamics and Chaos. Westview Press; 1994.
86. Khalil HK. Nonlinear systems. 3rd ed Prentice-Hall; 2002.
87. Xiu ZL, Zeng AP, Deckwer WD. Multiplicity and Stability Analysis of Microorganisms in Continuous Culture: Growth Inhibition. Biotechnology and bioengineering. 1998;57:251–261. [PubMed]
88. O’Donnell-Tormey J, Nathan CF, Lanks K, DeBoer CJ, de la Harpe J. Secretion of pyruvate. An antioxidant defense of mammalian cells. The Journal of experimental medicine. 1987;165(2):500–14. doi: 10.1084/jem.165.2.500 [PMC free article] [PubMed]
89. Duarte TM, Carinhas N, Barreiro LC, Carrondo MJT, Alves PM, Teixeira AP. Metabolic responses of CHO cells to limitation of key amino acids. Biotechnology and Bioengineering. 2014;111(10):2095–2106. doi: 10.1002/bit.25266 [PubMed]
90. Reinhart D, Damjanovic L, Kaisermayer C, Kunert R. Benchmarking of commercially available CHO cell culture media for antibody production. Applied Microbiology and Biotechnology. 2015;99(11):4645–4657. doi: 10.1007/s00253-015-6514-4 [PMC free article] [PubMed]
91. Martínez VS, Dietmair S, Quek LE, Hodson MP, Gray P, Nielsen LK. Flux balance analysis of CHO cells before and after a metabolic switch from lactate production to consumption. Biotechnology and bioengineering. 2013;110(2):660–6. doi: 10.1002/bit.24728 [PubMed]
92. Hayter PM, Curling EMA, Gould ML, Baines AJ, Jenkins N, Salmon I, et al. The effect of the dilution rate on CHO cell physiology and recombinant interferon-gamma production in glucose-limited chemostat culture. Biotechnology and Bioengineering. 1993;42(9):1077–1085. doi: 10.1002/bit.260420909 [PubMed]
93. Hu W, Berdugo C, Chalmers JJ. The potential of hydrodynamic damage to animal cells of industrial relevance: Current understanding. Cytotechnology. 2011;63(5):445–460. doi: 10.1007/s10616-011-9368-3 [PMC free article] [PubMed]
94. Nienow AW. Reactor engineering in large scale animal cell culture. Cytotechnology. 2006;50(1-3):9–33. doi: 10.1007/s10616-006-9005-8 [PMC free article] [PubMed]
95. Delvigne F, Zune Q, Lara AR, Al-Soud W, Sørensen SJ. Metabolic variability in bioprocessing: Implications of microbial phenotypic heterogeneity. Trends in Biotechnology. 2014;32(12):608–616. doi: 10.1016/j.tibtech.2014.10.002 [PubMed]
96. Sorrell MF, Tuma DJ. The functional implications of acetaldehyde binding to cell constituents. Annals of the New York Academy of Sciences. 1987;492(1):50–62. doi: 10.1111/j.1749-6632.1987.tb48652.x [PubMed]
97. Lieber CS. Metabolic effects of acetaldehyde. Biochemical Society Transactions. 1988;16(3):241–247. doi: 10.1042/bst0160241 [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science