CHOmpact: A reduced metabolic model of Chinese hamster ovary cells with enhanced interpretability

Abstract Metabolic modeling has emerged as a key tool for the characterization of biopharmaceutical cell culture processes. Metabolic models have also been instrumental in identifying genetic engineering targets and developing feeding strategies that optimize the growth and productivity of Chinese hamster ovary (CHO) cells. Despite their success, metabolic models of CHO cells still present considerable challenges. Genome‐scale metabolic models (GeMs) of CHO cells are very large (>6000 reactions) and are difficult to constrain to yield physiologically consistent flux distributions. The large scale of GeMs also makes the interpretation of their outputs difficult. To address these challenges, we have developed CHOmpact, a reduced metabolic network that encompasses 101 metabolites linked through 144 reactions. Our compact reaction network allows us to deploy robust, nonlinear optimization and ensure that the computed flux distributions are physiologically consistent. Furthermore, our CHOmpact model delivers enhanced interpretability of simulation results and has allowed us to identify the mechanisms governing shifts in the anaplerotic consumption of asparagine and glutamate as well as an important mechanism of ammonia detoxification within mitochondria. CHOmpact, thus, addresses key challenges of large‐scale metabolic models and will serve as a platform to develop dynamic metabolic models for the control and optimization of biopharmaceutical cell culture processes.


| INTRODUCTION
Production of recombinant proteins is known to compete with biomass synthesis for externally provided nutrients. This is particularly true for mammalian cell lines, such as Chinese hamster ovary (CHO) cells, which are the dominant host for the industrial production of therapeutic proteins (O'Flaherty et al., 2020). Metabolic modeling has become an essential tool for understanding resource allocation and, when coupled with advances in genome editing, for designing rational cell engineering strategies. Publication of the CHO-K1 genome, and the omics analyses this enabled, laid the foundation for systems-level understanding of this host. This knowledge has been reconstructed mathematically in a community-curated genomescale metabolic model (GeM) of CHO cells termed iCHO1766 (Hefzi et al., 2016). Crucially, iCHO1766 organized knowledge of all biochemical conversions, transport, and exchange reactions to create a large, interlinked network of metabolites and their associated reactions.
The inclusion of gene-protein associations provided a direct link between genes and metabolic reactions. Since then, significant expansions and improvements to iCHO1766 have been achieved, such as gap-filling studies that also removed dead-end reactions (Fouladiha et al., 2021), and the integration of a core protein secretory pathway, iCHO2048, enabling the computation of energetic costs and machinery demands of each secreted protein . iCHO2048 was subsequently used to direct host cell protein knockout studies, which resulted in increased recombinant protein productivity and a cleaner feedstock for downstream purification (Kol et al., 2020), highlighting the power of these models in identifying cellular engineering strategies.
The solution of GeMs, and any underdetermined metabolic model, relies on constraint-based methods, such as flux balance analysis (FBA), to predict steady-state intracellular flux distributions (Orth et al., 2010). Although FBA offers the advantage of not requiring detailed knowledge of enzymatic kinetic parameters, it does not return a unique set of intracellular flux values. In addition, the larger the metabolic network considered, the more difficult it becomes to interpret such predictions (Gardner & Boyle, 2017).
GeMs, therefore, require large datasets, preferably across different omics levels (e.g., metabolomic, transcriptomic) to increase confidence in results. This is also true for curating GeMs for specific cell lines, raising the need for characterization that goes beyond what is typically conducted in industrial settings.
Several algorithms have been developed to improve the predictive performance of GeMs, for example, by constraining the amount of carbon able to flowthrough reaction fluxes, based on the maximum amount carbon uptake by the cell (ccFBA) (Lularevic et al., 2019), taking into account the selective pressure that exists within cell cultures for fast-growing cell lines with a low enzyme usage (Lewis et al., 2010), or introducing enzyme capacity constraints (Yeo et al., 2020). Despite these advances, difficulties in handling the models and interpreting their intracellular flux predictions remain unaddressed challenges. An additional limitation is the computational difficulty in creating dynamic versions of GeMs that would reflect the nature of cell culture processes, although recent efforts coupling a CHO GeM with statistical models have yielded promising results in predicting the time evolution of extracellular (EC) amino acid concentrations (Martínez et al., 2015).
In this work, we introduce a reduced-scale metabolic model, CHOmpact, where the reaction network is based on the work by Carinhas et al. (2013) and has been augmented with a detailed description the aspartate-malate (Asp-Mal) shuttle, the urea cycle, de novo serine synthesis from glycolytic intermediates, and nucleotide sugar donor (NSD) biosynthesis. The resulting network comprises 101 metabolites and 144 reactions, which due to its compact nature, significantly enhances the interpretability of simulation results. The reduced scale of the reaction network also allows for more complex, nonlinear formulations of the objective function to be incorporated, compared to biomass maximization that is often employed in FBA of GeMs. Our proposed optimization framework allows solution across all phases of cell culture and provides insight into the dynamics of cellular metabolism. We envisage that the advantages presented by CHOmpact will enable the development of dynamic flux balance models that can serve as digital twins for the control and optimization of biopharmaceutical cell culture processes.
Triplicate cultures for each feeding regime were performed in orbitally shaken (140 rpm) 250 mL vented conical flasks (Corning) with a 50 mL working volume. The flasks were placed in a humidified incubator with 8% CO 2 and a temperature of 36.5°C. The basal medium for all cultures was CD CHO (Life Technologies) supplemented with 25 μM methionine sulfoximine (Sigma-Aldrich). All feeding regimes consisted in adding 10% v/v every 48 h of culture.
Feed C used commercial CD EfficientFeed™ C AGT™ (Invitrogen), whereas U and U40 provided amino acids beyond the amounts available in Feed C. Glucose and amino acid concentrations in the feeds are detailed in Kyriakopoulos and Kontoravdi (2014).

| Analytical methods
Viable and dead cell density was determined using the trypan blue dye exclusion method. mAb titer was determined using the BLItz ® system (Pall ForteBio). Time profiles for glucose, lactate, and ammonia were generated using the Bioprofile 400 analyzer (NOVA Biomedical). Residual amino acid profiles were quantified with the PicoTag method (Waters) on an Alliance HPLC instrument (Waters). EC pyruvate concentrations were determined with an enzyme assay kit (Abcam). mAb Fc glycoprofiling was performed with an automated sample preparation workflow (Stockmann et al., 2013) where mAb samples were affinity-purified from the cell culture supernatant with a 96-well protein G immunoglobulin G purification plate (Thermo Fisher Scientific). Glycans were released from the mAb through PNGase (Prozyme) digestion and labeled with 2-amino benzamide (Ludger). Labeled glycans were separated using hydrophilic interaction ultra performance liquid chromatography and quantified with fluorescence detection (Stöckmann et al., 2013). Glycans were initially assigned by comparing their Glucose Unit retention times with those available in the national institute for bioprocessing research & training GlycoBase 3.2 structural N-glycan library (Campbell et al., 2008). Glycan assignment was confirmed through weak anion exchange chromatography and quadrupole time-of-flight mass spectrometry on exoglycosidase-digested samples (Albrecht et al., 2014).

| Dry cell weight measurement
Duplicate cultures were harvested at Days 4 (mid-exponential) and 10 (stationary) for dry cell weight measurements. First, viable cell density was determined using trypan blue dye exclusion. Immediately after cell counting, 40 mL of the cultures was harvested and centrifuged at 1000g for 1 min in preweighed 50 mL falcon tubes.
After discarding the supernatant, cell pellets were washed once with 40 mL 0.9% w/v NaCl (Sigma-Aldrich) and centrifuged at 1000g for 1 min. The wash was discarded, and the cell pellet was left to dry in a nonhumidified incubator at 37°C until no changes in weight were observed. Tubes were weighed within 1 mg accuracy (ACCULAB; Sartorius).

| Data processing and analysis
The cell-specific rates for nutrient consumption and metabolite/ product secretion, q t ( ) i n , were calculated with Equation 1 (Sauer et al., 2000) is the consumed/produced amount of component i up to time t n and IVC t ( ) n is the integral of viable cells up to time t n .  (Fan et al., 2015;del Val, Fan, et al., 2016;).

| Flux balance model development
The CHOmpact flux balance model ( Figure 1) is based on previous work by Carinhas et al. (2013) and has been expanded to include the Asp-Mal shuttle (Mulukutla et al., 2012;Nolan & Lee, 2011), the urea cycle (Zamorano et al., 2010), de novo serine synthesis from glycolytic intermediates, and NSD biosynthesis (Kremkow & Lee, 2018). Details on the Asp-Mal shuttle were included with the aim of gaining further insight into glutamate and aspartate anaplerosis and cataplerosis. Reactions of the urea cycle, catalyzed by enzymes whose genes are present in CHO cells, were added to provide additional avenues for ammonia detoxification and better insight into arginine metabolism.
De novo serine synthesis was included to ensure that this prototrophic amino acid in CHO cells is not growth limiting.
NSD biosynthesis was included to gain insight into the metabolic burden of cellular and recombinant product glycosylation on CHO cell metabolism.
Manual curation of our FBA model was performed using the kyoto encyclopedia of genes and genomes (KEGG) database (Kanehisa et al., 2017(Kanehisa et al., , 2019 and the reference CHO-K1 and Cricetulus griseus genome annotations (Kremkow et al., 2015;Lewis et al., 2013;Rupp et al., 2018). Sequential reactions were combined into single reaction fluxes to reduce degrees of freedom within the model (Nolan & Lee, 2011). CHOmpact is comprised of material balances for 101 species, an additional equation that defines the consumption of ATP toward active amino acid transport, and 144 fluxes ( Figure 1 and Supporting Information: The amino acid composition of cellular proteins was computed from CHO cell proteomic data (Baycin-Hizal et al., 2012), as reported by del Val, Polizzi, et al. (2016). The glycan content of biomass has been included by using stoichiometric coefficients for cellular glycolipid and N-and O-linked protein glycosylation (del Val, Polizzi, et al., 2016). The stoichiometric equation for the cB72.3 humanized IgG4κ mAb was computed using the amino acid sequences for human IgG4 Fc (Heilig et al., 2003), the human kappa light chain constant fragment (Brady et al., 1991), and the variable heavy and light chain fragments for cB72.3 (Xiang et al., 1999), as shown in Supporting Information: Table 4. mAb glycoprofiling at three culture timepoints (192, 240, and 288 h) allowed us to calculate stoichiometric coefficients for NSD consumption toward mAb glycosylation across three culture intervals: 0-192 , 192-240 , and 240-288 h. These calculations were made with Equation 4, and the obtained stoichiometric values are presented in Supporting Information: Table 5.
F I G U R E 1 CHOmpact reaction network CHOmpact considers 101 species linked through 144 fluxes. Different colors indicate particular metabolic pathways: glycolysis, tricarboxylic acid, and oxidative phosphorylation (green), nucleotide sugar donor metabolism (red), aspartatemalate shuttle (orange), urea cycle (purple), amino acid and nucleotide metabolism (dark blue), and cycle fluxes that must be constrained and/or estimated during optimization (light blue).

| FBA solution: Nonlinear optimisation
As with most FBA models, no intracellular accumulation of species has been assumed in the material balances generated from our stoichiometric matrix, leading to a problem of the form S F × = 0, where S is the stoichiometric coefficient matrix presented in Supporting Information: Table 2 and F is the vector of unknown fluxes. Because the model contains more unknown fluxes (144) than equations (102), it is solved using constraint-based optimization (Banga, 2008).
Two constraint-based optimization strategies were used to solve CHOmpact (Table 1). The first maximizes the rate of biomass synthesis while maintaining the transport flux for all nutrients, metabolites, and product set to experimentally determined values.
Supporting Information: Table 6 presents the experimental data used for optimization. Reaction reversibility constraints, based on enzyme data available in KEGG (Kanehisa et al., 2017(Kanehisa et al., , 2019 and BRENDA (Jeske et al., 2019), were included as indicated in Supporting Information: Tables 1 and 2.
The optimization strategy proposed herein simultaneously maximizes the fluxes where ATP is synthesized, while minimizing the sum of squared intracellular fluxes. This objective function represents maximum energetic efficiency by the cells (Schuetz et al., 2007) and was used to ensure consistent directionality of central carbon metabolism fluxes.
Alongside maximizing the energetic efficiency of the cells, the sum of squared differences between measured and computed fluxes was minimized to ensure consistency between model results and experimental measurements. To avoid flux F 15 (Mal mit + NAD + ↔ OAA mit + NADH mit ) being bypassed by F 17 (Mal mit + NAD + → Pyr mit + CO 2 + NADH mit ), the F 17 /F 14 ratio was constrained to values within [0, 1] and minimized. If this constraint is not present, F 15 , may be reversed and possibly lead to stalling the tricarboxylic acid (TCA) cycle by the absence of carbon flowing through F 9 (AcCoA + OAA mit → Cit + CoASH). The secretion fluxes of nonmeasured by-products (except CO 2 ) were minimized, and flux reversibility constraints were also included within the CHOmpact optimization strategy. Our rationale for minimizing the secretion of nonmeasured fluxes intends to ensure that most of the carbon is destined toward energy production (correct TCA directionality) and that amino acids are mainly consumed toward biomass synthesis. It is worth noting that the growth rate maximization strategy that is prevalently used to solve flux balance models also implies by-product secretion minimization (i.e., it pushes nutrients toward biomass synthesis  LaNoue & Tischler, 1974). F 16 was constrained to be 1% of the glucose uptake flux (F 105 ), as previously measured for CHO cells (Ahn & Antoniewicz, 2013).
Both optimization strategies are outlined in Table 1, where F j are all intracellular fluxes, ATP synth. is the squared sum of ATP synthesis reaction fluxes, SSE is the sum of square errors (SSEs) between the T A B L E 1 Optimization strategies for flux balance model solution.   μ g optimizations were performed on CHOmpact using the biomass composition assumed by Hefzi et al. (2016) to discern whether it was the cause for deviations in predictive capabilities.

BM maximization Nonlinear optimization
These optimizations yielded deviations from the experimental data that are indistinguishable from those obtained with our assumed biomass composition (data not shown). These results are expected when considering that the differences between our assumed biomass composition and the one used in the GeM involve prototrophic amino acids only (Ala, Gln, Gly) (Supporting Information: Figure 3). 3.2 | CHOmpact identifies the source of μ g prediction inaccuracies Our proposed nonlinear optimization strategy enables the identification of nutrient uptake rates that lead to differences between calculated and experimental-specific growth rates. This is achieved by fixing the biomass growth rate and the uptake/secretion rates of Glc, The above strategy yields the results presented in Table 2, where the percent increases in specific uptake rates required to match the experimental values of μ g and q p are shown (positive/red values denote the percent increase in uptake fluxes required to match μ g and q p ). Table 2 shows that, for Car LP + NaBu, Selv Early Exp, Selv Late Exp, and Mart Cold 2, small or no increases in amino acid uptake rates are required to match μ g and q p . This is expected, considering the results of Figure 2, where the predicted values for μ g are matched or exceeded for these datasets.
Across all other datasets, substantial increases in amino acid uptake rates are required to match the experimental values for μ g and q p . Except F I G U R E 2 Comparison of experimentally determined and predicted specific growth rates (μ g ) for different Chinese hamster ovary cell lines and culture conditions. The dark blue bars present the experimentally determined μ g values reported by Carinhas et al. (2013), Martínez et al. (2015), and Selvarasu et al. (2012). The medium blue bars present predictions reported by Hefzi et al. (2016), and the light blue bars show the μ g values predicted with CHOmpact.
for both Mart Warm datasets, the amino acids which limit growth are auxotrophic (shown in Table 2) and, in all but one case (Car HP NaBu), multiple auxotrophic amino acids limit growth rate. Sharp deviations across multiple amino acids are observed for the Car LP and both Mart Warm datasets.
Underestimations for μ g arise from two possible sources: (i) either the assumed stoichiometric coefficients for auxotrophic amino acids in biomass are too high or (ii) the measured uptake rates for auxotrophic amino acids are underestimated. Errors in uptake rate measurements above 30% are unlikely (see Table 2); therefore, uncertainty in biomass weight and composition are likely responsible, especially when considering they are seldom measured for flux balance studies.  Table 2 and would improve the predictive capability of CHOmpact and Hefzi's GeM.
The results of Figure 2 and Table 2

| CHOmpact objective function robustness
The CHOmpact objective function is nonlinear, and when using a local nonlinear solver, such as gPROMS NLPSQP, there is uncertainty T A B L E 2 Percent increases in amino acid uptake rates required to match μ g and q p .

| CHOmpact facilitates easier interpretation of flux distributions
Our optimization strategy provides enhanced insight into metabolic dynamics by enabling the calculation of flux distributions during cell culture phases beyond exponential growth. In addition, our reduced reaction network simplifies model output interpretation to better relate the obtained flux distributions with cellular physiology. Conversely, fluxes through the TCA pathway increase with culture time, which occurs because reduced lactate production allows for more glycolysis-derived pyruvate to reach mitochondria. The high rates of glycolysis and lactate production observed during the initial growth phases of culture, often referred to as the Warburg effect (Vander Heiden et al., 2009), have been widely reported in CHO cells (Buchsteiner et al., 2018;Kelly et al., 2018).

| Flux distribution dynamics
The calculated fluxes through the Asp-Mal shuttle are defined by our imposed constraint on F 20 , which limits its value to the flux of

| Glutamate anaplerosis and cataplerosis
Glutamate can either be consumed toward TCA and energy production (anaplerosis) or produced from TCA metabolites for subsequent use in biomass synthesis (cataplerosis). In the context of our flux balance model, Glu anaplerosis is observed when more of this amino acid is transported into mitochondria (F 20 + F 30 ) than what is transported out, in the form of αKG. Net Glu cataplerosis is observed when less of it is transported into mitochondria than the αKG transported out. Importantly, however, anaplerosis and cataplerosis are driven by the balance of species along the TCA cycle, not solely amino acid metabolism.
Glu transport into/out of mitochondria is a consequence of carbon shortage/excess within the TCA (Owen et al., 2002). Figure 4 shows that, in Feed C and Feed U, net Glu cataplerosis is observed during the exponential growth intervals and is reversed during stationary phase, when net anaplerosis occurs. Feed U40 contrast with Feed C and Feed U by presenting Glu anaplerosis during early phases of culture (Figure 4-bottom).
Feed U presents considerable anaplerosis during the early stationary interval, where it accounts for 92.5% of total consumed Glu (Figure 4-bottom). This high anaplerosis is associated with increased Glu production from Gln (37.3% of total) and, to a lesser extent, from His (18.5% of total) (Figure 4-top). Interestingly, this is the only culture condition and interval where Gln is consumed from the extracellular environment (Supporting Information: Figure 1). In In Feed C, Glu cataplerosis is driven by the high Asn uptake rates observed during exponential growth, where F 31 provides a path for Asn overflow toward pyruvate via oxaloacetate (F 29 ) and malate (F 24 and F 21 ). An alternative cataplerotic pathway for Glu would be its direct synthesis from cytosolic αKG via a cytosol-localized EC 1.4.1.3 enzyme. Such an enzyme would consume NH 4 + and produce cytosolic NAD + (Yang et al., 2005) and could, thereby, reduce lactate production (Freund & Croughan, 2018). These properties would make a cytosolic version of EC 1.4.1.3 an interesting target for metabolic engineering; however, the endogenous version of this enzyme is not localized in the cytosol of CHO cells.
Prior flux balance work on standard (non-GS) CHO cells commonly reports net Glu anaplerosis during the exponential growth phase of cells, where high lactate production results in low TCA fluxes (Ahn & Antoniewicz, 2013). Glu anaplerosis is thought to be used by cells to replenish flux through TCA, and is also known to be a major source of ammonia production because much of this Glu anaplerosis is commonly described as the mitochondrial production of αKG from Glu (via F 23 - Figure 1) (Ahn & Antoniewicz, 2013;Mulukutla et al., 2012;Nicolae et al., 2014).
Crucially, the implicit preceding step is Glu transport into mitochondria, which occurs either through the Asp-Mal shuttle Glu/Asp Aralar1/Aralar2 antiporter (F 20 ) or the Asp-Mal shuttleindependent GCH1 Glu/H + antiporter (F 30  Across all feeds, the maximum proportion of Asn/Asp consumed toward biomass and mAb product is 28.8% for Feed U40 during early exponential phase ( Figure 5-bottom). These results show that Asn/ Asp are fed well beyond stoichiometric requirements for growth and product formation, and that anaplerosis provides an overflow pathway for when these nutrients are fed in excess. These results are consistent with prior work where excess Asn/Asp feeding has been found to increase Ala and Lac secretion by CHO cells (Calmels et al., 2019;Selvarasu et al., 2012).
The anaplerotic overflow pathway for Asn/Asp results in rapid uptake of Asn and a concomitant reduction in its concentration in the culture medium (Supporting Information: Figure 1). The low residual Asn concentrations observed in the Feed C culture, where Asn and Asp are fed at the lowest levels, could be interpreted as being growth limiting. However, our FBA results demonstrate the contrary: increasing Asn feeding was found to slow cell growth (Supporting Information: Figure 1), likely due to excess ammonia production. Interestingly, these flux values correlate closely with those of total Glu production (6.9, 5.5, and 6.1 nmol/10 6 cells/h for the corresponding culture intervals), indicating that Glu availability may regulate Asn/Asp anaplerosis.

| NH 4 + sources and sinks
NH 4 + is a key determinant of CHO cell culture performance because it is known to impact cell growth (Synoground et al., 2021;Wahrheit et al., 2014) and product quality (Borys et al., 1994;Hong et al., 2010).
In standard CHO cells, NH 4 + is mainly generated as a by-product of Gln anaplerosis (glutaminolysis) (Dean & Reddy, 2013;Hong et al., 2010;Wahrheit et al., 2014). Glutamine synthase (GS) cells satisfy their Gln requirements by producing it from Glu via ectopic GS expression. Despite considerable reductions in NH 4 + accumulation, GS-CHO cells still produce ammonia to levels that may impact product glycosylation (Borys et al., 1994;Hong et al., 2010), so it is therefore important to characterize the major sources and sinks of this key metabolite.
The dark blue bars in the top half of Figure 6 show that the majority (>50%) of ammonia is produced from asparagine (F 51 Figure 1). Additional sources of NH 4 + include Ser (F 36 ), Thr (F 37 ), and His (F 54 ), although to much lower levels. These results indicate that, in GS-CHO cells, NH 4 + production is mainly caused by Asn/Asp anaplerosis, which is consistent with previous work with GS-CHO cells (Calmels et al., 2019;Carinhas et al., 2013).
The total production rate of NH 4 + further confirms its link with Asn/Asp anaplerosis. The top half of Figure 6 shows that the total amount of NH 4 + produced by the cells increases with higher levels of Asn feeding. After the Mid Exponential interval, cells cultured with basal Asn feeding levels (Feed C) have NH 4 + production rates below 11.5 nmol/10 6 cells/h, whereas the Feed U and Feed U40 cultures (increasingly higher levels of Asn feeding) produce above 20 nmol/ 10 6 cells/h of NH 4 + .
The bottom half of Figure 6 (Heffner et al., 2020). It is worth noting that both enzymes are known to be expressed at low levels in CHO K1 cells, so their ability to contribute to ammonia detoxification remains to be confirmed.

| CONCLUDING REMARKS
We have presented a compact reaction network to describe the metabolism of mAb-producing CHO cells. Our reduced metabolic network (144 reactions) performs comparably with the iCHO1766 GeM (>6000 reactions) in predicting the growth rates of different CHO cell lines. Our FBA framework also allowed us to identify the absence of cellular weight and composition measurements as the most likely cause of inaccuracies in predicting growth rates. We have also presented a comprehensive optimization strategy that constrains the solution space to yield physiologically consistent flux distributions across all phases of cell culture.
When coupled with nonlinear optimization, CHOmpact greatly enhances the interpretability of metabolic flux distributions across different phases of cell culture. Our results provide insights into the mechanisms underlying Glu anaplerosis and its dependence on the uptake of Asn/Asp. We have also identified Asn/Asp as the key anaplerotic nutrients of GS-CHO cells, where they act as a considerable source of lactate during the early stages of culture.
Our results also show that Asn is the predominant source of NH 4 + across all culture conditions and that a major sink for this key metabolite is consumption within mitochondria. The presence of Asp within mitochondria determines whether this organelle is a source or Moving forward, the enhanced understanding of metabolic dynamics afforded by our CHOmpact reaction network and nonlinear optimization framework can be used to define feeding strategies that optimize cell culture performance. Furthermore, the compact size of our reaction network will also facilitate the creation of hybrid dynamic FBA/culture dynamics models, which can be used as digital twins for dynamic optimization and control of cell culture bioprocesses.

CONFLICT OF INTEREST STATEMENT
The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.