Computers and Chemical Engineering 154 (2021) 107471

Contents lists available at ScienceDirect

Computers and Chemical Engineering

journal homepage: www.elsevier.com/locate/compchemeng

Synergising stoichiometric modelling with artificial neural networks to

predict antibody glycosylation patterns in Chinese hamster ovary cells
Athanasios Antonakoudis a,1, Benjamin Strain a,1, Rodrigo Barbosa a, Ioscani Jimenez del
Val b, Cleo Kontoravdi a,∗
a
Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College London, London SW7 2AZ, United Kingdom
b
School of Chemical & Bioprocess Engineering, University College Dublin, Belfield D04 V1W8, Ireland

a r t i c l e i n f o a b s t r a c t

Article history: In-process quality control of biotherapeutics, such as monoclonal antibodies, requires computationally ef-
Received 16 February 2021 ficient process models that use readily measured process variables to compute product quality. Existing
Revised 18 July 2021
kinetic cell culture models can effectively describe the underlying mechanisms but require considerable
Accepted 31 July 2021
development and parameterisation effort. Stoichiometric models, on the other hand, provide a generic,
Available online 1 August 2021
parameter-free means for describing metabolic behaviour but do not extend to product quality predic-
Keywords: tion. We have overcome this limitation by integrating a stoichiometric model of Chinese hamster ovary
Constraint-based modelling (CHO) cell metabolism with an artificial neural network that uses the fluxes of nucleotide sugar donor
Flux balance analysis synthesis to compute the profile of Fc N-glycosylation, a critical quality attribute of antibody therapeu-
Artificial neural network tics. We demonstrate that this hybrid framework accurately computes glycan distribution profiles using a
Monoclonal antibody set of easy-to-obtain experimental data, thus providing a platform for process control applications.
Protein glycosylation
© 2021 Elsevier Ltd. All rights reserved.
Hybrid modelling

Article-at-a-glance files over four-time intervals, in 3 independent fed-batch experi-

ments, to within 0.27% error using a set of easy-to-obtain experi-
The glycosylation of monoclonal antibodies, a critical quality at- mental data. The presented framework provides a platform for pro-
tribute known to impact therapeutic efficacy and safety, can vary cess control, cell line selection and potential metabolic engineering
significantly during antibody manufacturing processes. As a result, applications.
monitoring and controlling the glycan distribution of therapeutic
antibodies is vital for meeting stringent product quality require- 1. Introduction
ments. Significantly, wet-lab methods to determine glycan distri-
butions on therapeutic proteins are time-consuming and resource- Therapeutic proteins have revolutionized the way we treat life-
intensive, meaning model-based solutions for the accurate predic- threatening diseases such as cancer thanks to their specificity, rel-
tion of glycan distributions represent significant industrial value. ative safety and high therapeutic efficacy. However, biological drug
Here, using evidence generated during an in-depth stochiomet- development and manufacture suffers from lack of integrative ap-
ric model comparison study, we present a kinetic parameter-free, proaches, resulting in costly development campaigns and often
computationally inexpensive hybrid modelling framework for the suboptimal processes in terms of flexibility, productivity and qual-
prediction of antibody glycosylation using commonly generated ex- ity (Croughan et al., 2015). These limitations, in turn, translate into
tracellular measurements. This framework feeds nucleotide sugar high production costs. The development of industrial-scale pro-
donor fluxes computed by CHOmpact, a small-scale stochiomet- cesses is typically based on well-established mAb production plat-
ric model of central Chinese Hamster Ovary (CHO) cell metabolism form operations that give a good starting point for high titres but
into an artificial neural network, which calculates glycan distribu- often require optimisation to meet product quality requirements.
tions on the protein of interest. We demonstrate that this hybrid The latter results in pre-designed (i.e., open loop) operational (pri-
framework accurately computes antibody glycan distribution pro- marily reactor feeding) strategies applicable to the cell line/product
in question. Often, the methodology for developing such optimised
∗ processes is based on extensive and, thus, time-consuming high-
Corresponding author at: Department of Chemical Engineering, Imperial College
London, London SW7 2AZ, UK. throughput experimentation. Furthermore, in absence of in-process
E-mail address: cleo.kontoravdi@imperial.ac.uk (C. Kontoravdi). control, manufacturing processes still exhibit batch-to-batch vari-
1
AA and BS contributed equally ability (Planinc et al., 2017; Schiestl et al., 2011). We and oth-

https://doi.org/10.1016/j.compchemeng.2021.107471
0098-1354/© 2021 Elsevier Ltd. All rights reserved.
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

ers have previously shown that adopting a kinetic model-based metabolite uptake and secretion rates and computes the fluxes of
approach, we can map a wide range of process conditions and nucleotide sugar donors (NSDs), the co-substrates of glycosylation.
rapidly explore and assess or design different operational strategies The latter is a critical quality attribute of mAbs and other glyco-
(Kotidis et al., 2019a, 2019b; Yang and Ierapetritou, 2021). Such a protein therapeutics such as erythropoietin, interferons and certain
knowledge-based platform is expected to be transferrable across fusion proteins. These fluxes are then fed to the ANN model of pro-
different cell lines and process conditions to enable fully rational tein N-glycosylation, which calculates the resulting glycan distribu-
decision-making during bioprocess development (as reviewed in tion on the protein of interest.
Kyriakopoulos et al., 2018). The use of stoichiometric modelling imparts several benefits:
Nonetheless, the need for re-parameterisation due to metabolic stoichiometric metabolic models are parameter-free, linear and
rewiring during cell culture (Ahn and Antoniewicz, 2011) limits the typically less computationally expensive than the respective ki-
effective use of fully kinetic models for online applications such netic metabolic models. Although they rely on the system being
as control, even when considering the relatively slow evolution of at quasi-steady state conditions, several advances have been made
most relevant biological changes. An alternative is the use of statis- towards constructing dynamic metabolic flux analysis models in-
tical models such as sequential multivariate tools with partial least cluding in CHO cell systems (e.g., Martínez et al., 2015). There
squares regression (PLSR), which have been combined with genetic is a vast choice of stoichiometric models for CHO cells ranging
algorithms to predict antibody titre and various quality indicators from small-scale representations of central carbon metabolic path-
in Chinese hamster ovary (CHO) cell cultures (Sokolov et al., 2018). ways (Carinhas et al., 2013) to genome-scale models (GeMs) of cell
Similarly, PLSR has been used to capture the information content metabolism, iCHO1766 (Hefzi et al., 2016), and, most recently, the
of a comprehensive metabolic dataset for four different fed-batch secretory pathway, iCHO2048s (Gutierrez et al., 2020). In this study
CHO cell culture processes and subsequently used the derived sta- we therefore sought to compare the suitability of each approach
tistical model to predict antibody glycosylation (Zürcher et al., before incorporation in the hybrid framework.
2020). Interestingly, the study concluded that only 25% of metabo- Overall, the proposed approach is easily adaptable to new
lites measured were sufficient to describe variations in glycosyla- systems by constraining the stoichiometric model and re-
tion, pointing to a significant reduction potential in terms of ex- parameterising the ANN component with the data at hand. We
perimental effort. suggest that both tasks are significantly less cumbersome than
Alternatively, machine learning techniques, such as artificial re-parameterising kinetic models, which would potentially require
neural networks (ANN), support vector machines, Kringing etc., can re-formulation to reflect the behaviour of new cell lines or pro-
be used to build data-driven models. For instance, antibody glyco- cesses. Towards the goal of process control, it can be readily in-
sylation profiles have been successfully simulated using an ANN formed by metabolite data generated inline using Raman spec-
trained on experimentally determined intracellular concentrations troscopy (Bhatia et al., 2018; Ryder, 2018). Raman can also rapidly
of nucleotides and nucleotide sugar donors (NSD). This model was provide feedback measurements for mAb titre and quality in the
able to accurately predict the time evolution of product glycoform form of glycosylation (André et al., 2015; Li et al., 2018; Zavala-
distribution, as well as the impact of NSD precursor molecules Ortiz et al., 2020). Together with novel process analytical tech-
and manganese supplementation (Kotidis and Kontoravdi, 2020). nologies, the proposed hybrid framework will therefore enable in-
Interestingly, and despite not usually being reliable for extrapo- process quality control.
lation, the ANN successfully predicted the glycoform distribution
for process conditions that were distinctly different from the train- 2. Materials and methods
ing dataset, demonstrating the potential of such an approach for
model-directed process optimisation. 2.1. Experimental dataset
Both statistical and machine learning approaches are less com-
putationally expensive than fully kinetic descriptions, can be de- The experimental data used herein were generated with the
veloped more readily using off-the-shelf software and, importantly, GS46 cell line expressing glutamine synthetase and the cB72.3
can describe processes and relationships the mechanisms of which chimeric IgG4 antibody (kindly donated by Lonza Biologics, Slough,
are poorly understood. However, a significant drawback of purely UK), as described in Kyriakopoulos and Kontoravdi (2014). Briefly,
data-driven approaches is that they can require large amounts the dataset includes three experimental conditions each studied in
of data produced under a wide range of experimental conditions triplicate 50 mL cultures growing on CD CHO medium (Life Tech-
to obtain reliable models, which, in turn, involves a substantial nologies, Paisley, UK), supplemented with 25 μM methionine sul-
amount of experimentation. A hybrid approach that couples a foximine (MSX, Sigma–Aldrich, Dorset, UK), maintained at 36.5°C,
mechanistic model with a data-driven algorithm therefore imparts 140 rpm in 8% CO2 humidified air. The cultures were supplemented
benefits in terms of computational efficiency while retaining some with 10% v/v of feed every second day starting on day 2. Three
mechanistic components that support transferability across exper- feeds were evaluated: commercially available CD EfficientFeedTM
imental systems. For instance, work has successfully utilised a hy- C AGTTM (Invitrogen, Paisley, UK) and two in-house feed, Feed U
brid model consisting of a genome scale stoichiometric model cou- and Feed U40 (composition detailed in Kyriakopoulos and Kon-
pled to a data-driven regression model to predict amino acid con- toravdi, 2014). Viable and dead cell densities were determined us-
sumption rates (Schinn et al., 2020). The predictions generated ing the trypan blue dye exclusion method and light microscopy.
with this hybrid approach were significantly improved compared Amino acid analysis was carried out using the PicoTag method
to the stoichiometric model alone and were able to accurately (Waters, Hertfordshire, UK) using phenyl-isothiocyanate for pre-
predict the consumption of amino acids highly abundant in re- column derivatization of all amino acids apart from cysteine.
combinant antibodies, allowing the control of nutrient feeding to Quantitation was performed on the HPLC Alliance system (Wa-
avoid premature nutrient depletion. A review of recent advances ters, UK) with UV detection at 254 nm. The concentrations of
in hybrid modelling applied to bioprocesses and biological systems lactate, ammonia and glucose were determined with a Biopro-
more broadly is presented in Antonakoudis et al. (2020). file 400 analyzer (NOVA Biomedical, Waltham, MA). Pyruvate was
Herein, we propose a hybrid modelling framework in which quantified using an enzymatic assay kit (Abcam, Cambridge, UK).
mechanistic information is organised in the form of a stoichiomet- mAb titre was determined using the BLItz® system (Pall ForteBio,
ric model while mAb quality is modelled using an ANN, as illus- Portsmouth, UK). The glycan profile of the mAb Fc was performed
trated in Fig. 1. The stoichiometric model uses data in the form of at the National Institute for Bioprocessing Research and Train-

2
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Fig. 1. Hybrid modelling strategy: step 1 involves the evaluation of small-scale and genome-scale stoichiometric models. Step 2 uses the metabolic fluxes of nucleotide sugar
donors (NSD) computed by the stoichiometric model as inputs to an artificial neural network, which computes the glycan distribution of the produced antibody.

ing (Dublin, Ireland) with an automated sample preparation work- consumption/production in each culture period for each flask.
flow using ultra-performance hydrophilic interaction chromatogra- These bounds are outlined in Supplementary Table 1.
phy (Stöckmann et al., 2013).
2.3. Small-scale stoichiometric model
2.2. Data processing
The small-scale stoichiometric model, CHOmpact, uses the
Cell culture data was processed for use in stochiometric models metabolic reaction network proposed by Carinhas and co-workers
as follows: for GS-CHO cells (Carinhas et al., 2013). This has been augmented
with a detailed description of the aspartate-malate (Asp-Mal) shut-
a) IVCD, representing the cumulative amount of biomass and the
tle, the urea cycle, de novo serine synthesis from glycolytic in-
time that cells are viable per system volume (L), was calculated
termediates and nucleotide sugar synthesis to give a network
as following:
comprising 101 metabolites and 144 reactions. The stoichiometric
IV C Dt = IV C Dt−1 + 0.5 × (V C Dt + V C Dt−1 ) t biomass composition is based on Hefzi et al. (2016). Two optimi-
sation strategies have been used to solve the FBA problem: in the
Where IV C Dt is the integral of viable cell density ( cel lLs·hr ) up to
first, the objective function involves the maximisation of biomass
time t (h ), V C Dt is the viable cell density ( celLl s ) at time t and t
synthesis, while the second maximises recombinant protein syn-
is the corresponding time interval.
thesis. All optimisations were performed using the nonlinear pro-
a) The concentration of components following a feeding event was gramming sequential quadratic programming (NLPSQP) solver built
calculated as follows: into gPROMS ModelBuilder v5.1.1 (Process Systems Enterprise).
FIN × Ciin + Ci t−1 × Vt−1
Ci = 2.4. Product-specific genome scale model construction
Vt
Where Ci is concentration of component i (mM ), FIN is the feed An expanded iCHO2048s product-specific genome scale model
volume (L ), Ciin is the concentration of component i in FIN , Cit−1 (GeM) was constructed following the methodology outlined by
is the concentration (mM ) of component i and Vt−1 is the system Gutierrez et al., (2020). In brief, information regarding the secreted
volume (L ) before each feeding event, respectively. Vt is the volume product was obtained (Kyriakopoulos, 2014) including amino acid
of the system after the feeding event FIN . composition, presence of a signal peptide, number of disulphide
a) Specific rates of consumption/production were calculated as fol- bonds, number of core N-linked glycans, and molecular weight.
lowing: This information was used to add the appropriate secretory path-
way reactions to the model.
Ci
qi =
IV CD 2.5. Nitrogen constrained FBA
Where qi is the specific consumption/production of compo-
nent i ( mmol
cel l .h
). We augmented the carbon constraint FBA (ccFBA) algorithm
For input into the genome scale model, units were converted (Lularevic et al., 2019) with nitrogen atoms, resulting in nitrogen
to mMol/gDCW/hr, as this is the most widely used unit for FBA. constraint FBA (ncFBA). The algorithm is the same as ccFBA with
The dry cell weight (DCW) of the CHO cell was assumed to the difference being that instead of using the total uptake of car-
be 0.315 mg/106 cells (Kyriakopoulos and Kontoravdi, 2014). Up- bons to set the lower and upper bounds of the reactions, we use
per bounds were calculated as the maximum specific consump- the total sum of the nitrogen atoms entering the cell. The algo-
tion/production rate and lower bounds as the minimum specific rithm first calculates the total nitrogen uptake from the cell (based

3
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Fig. 2. Comparison of upper bounds with the ccFBA and ncFBA algorithms for nine enzymatic reactions that showed the greatest level of reduction after imposing the
nitrogen constraint.

Table 1 method that constrained each reaction the most was selected on
Error between ANN results and experimental values for
a reaction-by-reaction basis.
each of the four cell culture intervals. Training error is the
average of both the Feed C and Feed U. The prediction er-
ror refers to the Feed U40 dataset. 2.7. Flux sampling
Days 0-8 8-10 10-12 12+
Artificial centering hit-and-run (ACHR) Monte Carlo sampling
Training Error (%) 0.14 0.4 0.27 0.07
Prediction error (%) 0.11 0.22 0.13 0.03
(Schellenberger and Palsson, 2009) was used to sample the
metabolic flux solution space of iCHO1766 and iCHO2048 using the
sample function in the flux analysis submodule of COBRApy. The
solution space was sampled 50,0 0 0,0 0 0 times, of which 50 0 0 were
on the lower bounds of the exchange reactions) and then each re-
stored and the rest were discarded. Both models were constrained
action’s upper bound is the ratio of the total nitrogen uptake and
using the same dataset used during the FBA.
the nitrogen atoms produced by the reaction:

Nitrogen uptake
vupper
j
= 2.8. Development of ANN model
produced nitrogen
Nitrogen is abundant in amino acids and, thus, constraining Based on the work by Kotidis and Kontoravdi (2020), we de-
intracellular reactions based on the total nitrogen uptake led to veloped an ANN which considers the NSD fluxes computed by the
a number of reactions having their bounds reduced even further stoichiometric model as inputs and the glycan distribution of the
when compared to ccFBA, as can be seen in Fig. 2. protein of interest, IgG in our case, as outputs. We used the NSD
flux and glycan data from the Feed C and Feed U datasets to train
2.6. FBA solution of iCHO2048s and iCHO1766 the model and the data from Feed U40, which is an independent
dataset, for testing the ANN model’s prediction accuracy. The ex-
Both the product-specific GeM (iCHO2048) and original CHO perimental glycan data are in the form of glycoform distribution
cell GeM (iCHO1766) were constrained for each time-period us- (%). For the ANN development, they converted to glycan secretion
ing the minimum and maximum experimentally measured uptake rates across specific time intervals to match the respective NSD
rates of nutrients/metabolites for the Feed C dataset (Table 1). fluxes. To do so we took into account the specific IgG productiv-
When biomass growth rate was used as the objective function, ex- ity and difference in glycan distribution for each time interval as
perimentally measured specific productivity rates were added as proposed by del Val et al., (2016). The culture period was split into
constraints. When the mAb productivity was defined as the ob- four intervals: exponential phase 0-8 days, production phase 8-10
jective function, experimentally measured biomass growth rates days, late production phase 10-12 days and death phase 12+ days.
were added as constraints. All GeM FBA problems were solved This generated 16 input-output samples for training and eight sam-
using the COBRApy (Ebrahim et al., 2013) in Python 3.8.6 and ples for prediction.
the glpk solver. Parsimonious FBA (pFBA) (Lewis et al., 2010) was The formulation of an ANN model typically consists of four
solved using the inbuilt function within COBRApy. Carbon con- main steps: the initialization where random values are assigned to
straint FBA (ccFBA) was calculated following the methodology out- weights, forward propagation where we are propagating through
lined by Lularevic et al. (2019), where the maximum flux through the neural network to calculate the output values using the acti-
each intracellular reaction is calculated based on the total amount vation function, the calculation of the loss function where in each
of carbon available from uptake rates and the number of carbon step we determine if the difference of input to the output values
atoms participating in each reaction. When calculating nitrogen is greater than our objective and backward propagation where the
constraint FBA (ncFBA), this same methodology was applied to ni- weight parameters are updated using the derivative of the error
trogen uptake rates and number of nitrogen atoms participating in function based on the calculated error. These steps are repeated
each reaction. When combining ncFBA with ccFBA (nccFBA), the iteratively until the loss function value is acceptable. For the con-

4
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Table 2
Summary of different FBA methods applied to the GeM.

Method Description Reference

Flux balance analysis (FBA) A linear programming method that finds the flux distributions (Orth, Jeffrey D., Ines Thiele, 2010)
that optimise a selected objective function, subject to constrains.
Parsimonious enzyme usage FBA (pFBA) Bi-level optimisation problem that minimises enzyme-associated (Lewis et al., 2010)
fluxes, subject to an optimal objective function, theoretically
improving the biological relevance of predictions.
Carbon constrained FBA (ccFBA) Refines flux range predictions by constraining the flux of carbon (Lularevic et al., 2019)
through different reactions, based on the maximum amount of
carbon taken up by the cells.
Nitrogen/Carbon constrained FBA (nccFBA) A novel extension of ccFBA that follows the same rationale but This work
also considers nitrogen uptake by the cell, constraining each
reaction by either carbon or nitrogen availability depending on
which is most conservative.

struction, training, and validation of the ANN, we used the Tensor- timisation problem that minimizes enzyme-associated fluxes, sub-
Flow (v2.4) library in Python 3.8. The first step was to normalize ject to optimal biomass, thereby predicting the most stoichiomet-
the input and output data, for which we used the following equa- rically efficient pathways and thus theoretically improving the bio-
tion: logical relevance of predictions (Lewis et al., 2010).
  Carbon constrained FBA (ccFBA) was evaluated as an alternative
xi − min xi
xinorm =     methodology to obtain flux distributions across the reaction net-
max xi − min xi work. This technique refines flux range predictions by constrain-
ing the amount of carbon able to flow through reactions based
The ANN is sequential and made of 3 layers. First, the input
on the maximum carbon uptake by the cell. This approach helps
layer is a Gaussian noise layer, for which the function is included
force flux predictions to physiologically meaningful ranges and has
in TensorFlow. Gaussian noise helps mitigate overfitting by adding
been shown to substantially improve the accuracy of predicted flux
noise to the inputs of the ANN and creating artificial points on a
values compared with standard FBA (Lularevic et al., 2019). Here,
given error distribution. Then, there is a hidden layer with 50 neu-
we also apply a novel extension of ccFBA by introducing nitro-
rons, which, given the inputs, uses a Sigmoid activation function to
gen/carbon constrained FBA (nccFBA). This approach follows the
calculate the outputs. We did not use more than one intermediate
same methodology as ccFBA but also considers nitrogen uptake by
hidden layer to avoid overfitting. Finally, the output layer, which
the cell, constraining each reaction by whichever element limits
consists of eight outputs, provides the percentage of each glycan
the reaction the most. By closing the elemental balance between
on the antibody product. The number of training iterations was set
the secretion and uptake of these elements, these methods should
to 200, which was deemed sufficient for achieving the minimum
theoretically refine physiological relevance of flux range predic-
error. As a loss function for our model, we made use of the mean
tions and improve model performance.
absolute error function in TensorFlow, which computes the mean
of absolute difference between predicted and experimentally deter- 3.1.2. Biomass accumulation optimisation
mined glycan distribution values. Due to the small number of sam- Where biomass growth rate was the objective function, CHOm-
ples and to avoid overfitting, the ANN was chosen to have only one pact was the most accurate in predicting growth and productiv-
hidden layer with 50 neurons. The number of neurons and epochs ity rates at all time points except the decline phase, where it
was calculated iteratively until the error was minimised. After 200 failed to capture the decline phase of cellular growth (Fig. 3A).
epochs the model results exhibited negligible training and predic- In late stage fed-batch culture, cell growth is typically not limited
tion error as shown in Table 1. by nutrient availability but rather inhibited by the accumulation of
metabolic by-products, such as ammonia (Kyriakopoulos and Kon-
3. Results toravdi, 2014). Growth inhibition is not accounted for in stoichio-
metric models. It is therefore unsurprising that, given that the cells
3.1. Comparison of constraint-based modelling approaches continue to uptake nutrients, the model overpredicts cell growth.
At all time periods, both genome-scale models provided com-
3.1.1. Comparison of phenotype predictive performance parable predictions, with iCHO2048 displaying marginally lower
In order to assess which constraint-based modelling approach predictions in all instances. This was expected as models are
to progress for hybrid modelling, FBA optimising either maxi- highly similar, with the added secretory reactions within iCHO2048
mum biomass accumulation or product formation was performed model providing a larger drain on cellular resources compared to
on the CHOmpact stochiometric reconstruction of central carbon the pseudo-IgG production reaction within iCHO1766, thus reduc-
metabolism, the original CHO cell GeM, iCHO1766 (Hefzi et al., ing biomass predictions. FBA and pFBA provided near identical pre-
2016), and a product-specific CHO GeM, iCHO2048 (Gutierrez et al., dictions. The biomass maximisation objective is thought to best
2020). Models were constrained using the same industrially rel- describe the exponential phase of the culture. However, in both
evant Feed C dataset of measured metabolites and compared to cases predictions were significantly higher than experimental val-
the experimentally determined values. In an attempt to overcome ues for that phase. This overprediction of growth rate has been
the limitations of basic FBA applied to GeMs, which include over- previously reported in the literature (Hefzi et al., 2016), and is
estimation of fluxes, a wide range of optimal solutions and fu- likely attributable to the high nutrient uptake rates that result in
tile metabolic cycles, several different FBA methodologies (sum- a large solution space. In contrast, both ccFBA and nccFBA heavily
marised in Table 2) were applied to the GeMs. The first method- constrained predictions across all time periods. As a result, these
ology, termed parsimonious enzyme usage FBA (pFBA), refines flux methods underpredicted biomass accumulation during the expo-
range predictions by taking into account the selective pressure that nential phase of cellular growth but provided more biologically rel-
exists within cell cultures for fast-growing cell lines with low en- evant predictions in later intervals. This may be partially due to the
zyme usage. This is represented mathematically via a bi-level op- value of dry cell weight used in the model being the higher of the

5
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Fig. 3. Comparison of results for A) biomass growth rate and B) specific antibody productivity for the small-scale metabolic model (CHOmpact) and the GeMs (iCHO1766
and iCHO2048) using FBA algorithms presented in Table 2. Error bars represent the experimental standard deviation.

two values determined by Kyriakopoulos and Kontoravdi (2014) (a improve predictions, they may in fact over-constrain the solu-
lower value of 219 mg/106 cells was experimentally determined for tion space. Again, much like biomass optimisation, the CHOmpact
the exponential growth phase). Significantly, the CHOmpact model metabolic model consistently provided the most accurate solutions
outperformed all methodologies applied to the two GeMs in terms across all time periods.
of quantitative biomass growth rate prediction.
3.1.4. Comparison of flux distributions
3.1.3. Product demand optimisation As our proposed hybrid framework involves the use of pre-
For the antibody productivity optimisation, CHOmpact accu- dicted intracellular fluxes, it was deemed necessary to further in-
rately predicted specific productivity rates across all time periods vestigate differences between these stoichiometric modelling ap-
(Fig. 3B). GeM predictions followed a similar trend to biomass opti- proaches by examining the predicted flux values for key shared in-
misation, where iCHO2048 showed slightly reduced, and thus gen- tracellular reactions. This comparison was carried out for the case
erally improved, predictions compared to iCHO1766. Additionally, of biomass optimisation during the exponential phase of growth.
FBA and pFBA predictions were significantly higher than experi- As shown in Fig. 4A, each modelling methodology gave substan-
mental values during exponential phase but showed improved rel- tially different predicted flux results. While the small-scale model
evance during later culture phases. This poor performance during displayed relatively consistent flux values through reactions, GeM
exponential phase was particularly surprising for iCHO2048 where approaches displayed apparently sporadic and non-sensical flux
it was expected that the addition of product-specific reactions into values. Of particular interest was the fact that the basic FBA ap-
the GeM would improve the biological relevance of productivity plied to iCHO1766 predicted zero flux for the HEX1 reaction. This
predictions by being able to better capture the cost of protein reaction, involving the phosphorylation of glucose to glucose-6-
secretion. The results highlight that the addition of further bio- phosphate, is critical in allowing glucose to enter central carbon
chemical reactions to a GeM does not necessarily guarantee bio- metabolism. Further investigation revealed the FBA was introduc-
logically relevant predictions in the absence of sufficient experi- ing glucose into central carbon metabolism via the hydrogenation
mental measurements for constraining the solution space. During of glucose to D-sorbitol, dehydrogenation of D-sorbitol to fruc-
exponential growth and stationary phases where a biomass con- tose and subsequent phosphorylation of fructose to fructose-6-
straint was present, both ccFBA and nccFBA predicted zero flux phosphate (Fig. 4B), representing a ‘non-typical’ flux distribution.
due to infeasibility as a result of an over-constrained solution Significantly, for these same reactions pFBA displayed vastly differ-
space. Due to the complete closure of the carbon/ nitrogen bal- ent flux predictions that were more in line with previous studies
ance, this likely indicates that not all carbon or nitrogen sources (Ahn and Antoniewicz, 2011), following the expected pathway via
were considered by experimentally measured uptake rates, mean- glucose-6-phosphate (Fig. 4B). This behaviour was likely thanks to
ing there was not enough carbon or nitrogen available for the the pFBA forcing the model to take more stoichiometrically effi-
model to allocate to antibody synthesis. As such, following the cient pathways.
method outlined by Lularevic et al. (2019), a relaxation factor was To evaluate if the inconsistencies in flux predictions for key re-
introduced that increased the amount of available carbon/nitrogen. actions between GeMs were also a global phenomenon, a mul-
To achieve a feasible solution, large relaxation factors of six and tivariate analysis was performed to assess whether the flux dis-
nine had to be applied to ccFBA and nccFBA, respectively, indicat- tributions of 6463 shared reactions correlated between iCHO2048
ing how over-constrained the problem was. This highlights how and iCHO1766 when using a variety of analysis technique (Fig. 4C).
resource-intensive antibody production is and further emphasises Considering the high level of similarity in terms of both model
that, while theoretically ccFBA and nccFBA methodologies should structure, phenotypic predictions, and flux solution space for key

6
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Fig. 4. A) Heatmap of predicted single point intracellular fluxes for shared key intracellular reactions computed with each FBA approach constrained using the same dataset.
B) Visual representation of the differing pathways FBA and pFBA analysis method take to generate fructose-6-phospate from glucose within iCHO1766. Flux value is repre-
sented by colour and size of arrow, where grey arrows indicate 0 flux. C) Multivariate analysis showing correlations of 6463 shared reaction fluxes between iCHO2048 and
iCHO1766 for each FBA approach.

reactions between iCHO2048 and iCHO1766, it was expected that vestigate this, flux sampling may be used to explore the feasible
when using the same analysis technique, model-wide flux distri- solution space by generating probability distributions of steady-
butions should strongly correlation between models. Significantly, state reaction fluxes. As demonstrated in Fig. 5, there is signifi-
however, this was not true in many cases, with only weak to mod- cant overlap of the solution space for iCHO1766 and iCHO2048 for
erate correlations reported for several of the analysis techniques. key reactions such as HEX1 (Fig. 5). As demonstrated above how-
For instance, flux corelations between models for FBA, ccFBA and ever, this does not translate into similar single point flux predic-
nccFBA techniques were not strong, with R values of 0.67, 0.39 and tions when using FBA (Fig. 4A). This arguably makes FBA applied
0.25 respectively. Conversely, pFBA correlated strongly between to GeMs less attractive for hybrid modelling efforts in which pre-
models (R= 0.99), which again may be attributed to the bilevel dicted intracellular fluxes are used as inputs for data-driven mod-
problem involved in pFBA forcing fluxes to a more stoichiometri- els, where a lack of intra-cellular flux prediction consistency may
cally efficient state, resulting in consistency across models. significantly impact the machine learning component of the hybrid
These intracellular flux prediction irregularities highlight a key modelling framework. On the other hand, pFBA generated more
issue with the use of FBA applied to large underdetermined promising results as, while it does not strictly produce a unique
metabolic networks, where it is possible for a vast number of solution (Kim et al., 2016), it does appear to improve the biological
flux distribution combinations to lead to the same optimal objec- relevance and consistency of intracellular flux predictions. While
tive function prediction (Orth, Jeffrey D., Ines Thiele, 2010). To in- other techniques such as geometric FBA (Smallbone and Simeoni-

7
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Fig. 5. Flux sampling distribution of key intracellular reactions for iCHO1766 (blue) and iCHO2048 (red). (For interpretation of the references to color in this figure legend,
the reader is referred to the web version of this article)

dis, 2009) have been developed that provide a true unique solu- constrain the large GeMs. While enzyme constrains could have
tion and hence flux distribution consistency, these are associated been added to the flux balance analysis framework (ecFBA) to try
with slow computational times, negating some of the benefits of and improve GeM performance (Yeo et al., 2020), this would have
the proposed hybrid framework. As a result, if progressing the CHO required significant amounts of kinetic data and so was deemed
cell GeM for future hybrid modelling efforts, using pFBA or flux impractical. Fundamentally, this highlights that, despite the large
sampling would be advisable. level of attention that CHO GeMs have received in recent years,
An alternative would be to constrain reactions to thermody- smaller scale metabolic models may provide a more appropriate
namically feasible ranges by integrating Gibbs free energy data into solution when limited data is available. As such, considering the
the model. This would remove thermodynamically infeasible reac- consistently accurate performance of CHOmpact, it was selected to
tions, improving the biological realism of flux predictions. While be progressed for hybrid-modelling studies.
this has been successfully achieved within smaller scale GeMs
(Fleming et al., 2009; Henry et al., 2007; Rui et al., 2015; Salvy and
Hatzimanikatis, 2020), it has yet to be effectively realised for the
3.2. Hybrid stoichiometric-ANN model
CHO cell GeM, due to the vast effort required for incorporating
Gibbs free energy data. The above observations also highlight the
NSDs are synthesised intracellularly using key nutrients sup-
requirement of additional measurements, such as oxygen uptake
plied in the cell culture media and are then transported to the
rate and CO2 production rate, to constrain the flux more accurately
ER and Golgi apparatus where a sequence of enzymatic reac-
through central carbon metabolism. In particular, CO2 production
tions takes place to yield the glycan chain attached to the crys-
rate may be used as a proxy of the flux through the TCA cycle.
tallisable fragment (Fc) of the antibody product (summarised in
Fig. 6). Due to the promiscuous nature of the glycosylation en-
3.1.5. Model selection for progression zymes, the reaction network in the Golgi apparatus is initially di-
Despite utilising multiple alternative FBA methodologies to im- vergent (Jimenez del Val et al., 2011), giving rise to a distribution of
prove GeM performance, the small-scale CHOmpact model consis- final glycan structures on the secreted antibody product. The inher-
tently outperformed the GeMs in terms of quantitative biomass ent complexity of the glycosylation process previously prompted
and productivity predictions. Moreover, flux distributions from this us to investigate the potential of an ANN to describe it without
small model were more tractable, with intra-cellular flux predic- the need for extensive kinetic parameter estimation required for
tions in a range that was perceived to be more biologically rel- mechanistic models. Based on this earlier work by Kotidis and Kon-
evant. This disparity in performance was an indication that the toravdi (2020), herein we developed an ANN that takes the NSD
metabolite consumption dataset we used, that is typical of datasets fluxes calculated using the CHOmpact as inputs and calculates the
generated in industrial settings, was not sufficient to effectively glycan distribution of the secreted antibody product.

8
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Fig. 6. Summary of how cell culture components are metabolised to produce NSDs, which, in turn are transported to the ER and Golgi apparatus where they react with the
growing oligosaccharide chain on the crystallisable fragment (Fc) of the antibody product. Created with biorender.com.

Fig. 7. Comparison of ANN model simulation results (pattern bars) with the training dataset that comprises data from the Feed C and Feed U experiments (solid bars) across
four cell culture time intervals (Kyriakopoulos and Kontoravdi, 2014).

9
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Fig. 8. Comparison of ANN model predictions (pattern bars) to independent dataset Feed U40 (solid bars) across four cell culture intervals (Kyriakopoulos and Kon-
toravdi, 2014).

The ANN closely matched the training dataset (from experi- Declaration of Competing Interest
ments Feed C and Feed U) across all four cell culture intervals con-
sidered as seen in Fig. 7. We then compared the ANN predictions The authors have no conflict of interest to declare.
against an independent dataset (Feed U40). As shown in Fig. 8, the
model correctly predicted the mAb glycoform distribution based CRediT authorship contribution statement
on the calculated NSD fluxes for the majority of structures. The
model exhibits a reasonable level of extrapolation ability, but dis- Athanasios Antonakoudis: Conceptualization, Methodology, Vi-
played discrepancy in the percentages of a1,6G1F and a1,3GF gly- sualization, Writing – original draft. Benjamin Strain: Concep-
cans, which were significantly different from the training dataset tualization, Data curation, Methodology, Visualization, Writing –
during the intervals of days 0-8 and 12+. This could be due to original draft. Rodrigo Barbosa: Conceptualization, Methodology.
these structures being in higher abundance in earlier time inter- Ioscani Jimenez del Val: Methodology, Writing – review & edit-
vals in the training dataset. Overall, despite the small sample num- ing. Cleo Kontoravdi: Conceptualization, Methodology, Supervi-
ber, the ANN performs yields quantitatively reliable results. Given sion, Writing – review & editing.
a richer training dataset, its predictive capability is expected to in-
crease further.
Acknowledgments

AA thanks the UK Engineering and Physical Sciences Research

Council for his studentship. BS and RB thank the UK Biotechnology
4. Conclusions and Biological Sciences Research Council and GlaxoSmithKline for
their studentships. IJV acknowledges the financial support of Sci-
The presented hybrid modelling framework predicts the glycan ence Foundation Ireland’s Solid State Pharmaceutical Cluster (SSPC)
distribution of the secreted antibody product given a dataset of Research Centre (Grant 12/RC/2275_P2). The authors wish to thank
common experimentally measured concentrations of key cell cul- Simone Albrecht, Henning Stockmann and Pauline M Rudd for per-
ture metabolites in the extracellular environment. This is achieved forming the glycomic analysis and Pavlos Kotidis for his help with
through a novel combination of a stoichiometric metabolic model the development of the ANN.
of CHO cells and an artificial neural network describing the glyco-
sylation process of the recombinant product. An in-depth compar- References
ison of small-scale and genome-scale metabolic models and FBA
methodologies revealed that the small-scale network is better able Ahn, W.S., Antoniewicz, M.R., 2011. Metabolic flux analysis of CHO cells at growth
and non-growth phases using isotopic tracers and mass spectrometry. Metab.
to accurately describe the experimental system across cell culture
Eng. 13, 598–609. doi:10.1016/j.ymben.2011.07.002.
duration. Overall, the framework considerably reduces the compu- André, S., Cristau, L.Saint, Gaillard, S., Devos, O., Calvosa, É., Duponchel, L., 2015. In-
tational expense of model parameterisation and execution com- line and real-time prediction of recombinant antibody titer by in situ Raman
pared to fully mechanistic alternatives. Additionally, it can be eas- spectroscopy. Anal. Chim. Acta 892, 148–152. doi:10.1016/j.aca.2015.08.050.
Antonakoudis, A., Barbosa, R., Kotidis, P., Kontoravdi, C., 2020. The era of big data:
ily trained for new conditions and systems, even in-process, setting genome-scale modelling meets machine learning. Comput. Struct. Biotechnol. J.
the foundations for online control applications. doi:10.1016/j.csbj.2020.10.011.

10
A. Antonakoudis, B. Strain, R. Barbosa et al. Computers and Chemical Engineering 154 (2021) 107471

Bhatia, H., Mehdizadeh, H., Drapeau, D., Yoon, S., 2018. In-line monitoring of amino consistent with computed optimal growth from genome-scale models. Mol.
acids in mammalian cell cultures using raman spectroscopy and multivariate Syst. Biol. 6. doi:10.1038/msb.2010.47.
chemometrics models. Eng. Life Sci. 18, 55–61. doi:10.10 02/elsc.20170 0 084. Li, M.-Y., Ebel, B., Paris, C., Chauchard, F., Guedon, E., Marc, A., 2018. Real-time
Carinhas, N., Duarte, T.M., Barreiro, L.C., Carrondo, M.J.T., Alves, P.M., Teixeira, A.P., monitoring of antibody glycosylation site occupancy by in situ Raman spec-
2013. Metabolic signatures of GS-CHO cell clones associated with butyrate troscopy during bioreactor CHO cell cultures. Biotechnol. Prog. 34, 486–493.
treatment and culture phase transition. Biotechnol. Bioeng. 110, 3244–3257. doi:10.1002/btpr.2604.
doi:10.1002/bit.24983. Lularevic, M., Racher, A.J., Jaques, C., Kiparissides, A., 2019. Improving the accu-
Croughan, M.S., Konstantinov, K.B., Cooney, C., 2015. The future of industrial biopro- racy of flux balance analysis through the implementation of carbon availabil-
cessing: batch or continuous? Biotechnol. Bioeng. 112, 648–651. doi:10.1002/bit. ity constraints for intracellular reactions. Biotechnol. Bioeng. 116, 2339–2352.
25529. doi:10.1002/bit.27025.
Del Val, I.J., Polizzi, K.M., Kontoravdi, C., 2016. A theoretical estimate for nucleotide Martínez, V.S., Buchsteiner, M., Gray, P., Nielsen, L.K., Quek, L.-E., 2015. Dynamic
sugar demand towards Chinese Hamster Ovary cellular glycosylation. Sci. Rep. metabolic flux analysis using B-splines to study the effects of temperature shift
6, 1–15. doi:10.1038/srep28547. on CHO cell metabolism. Metab. Eng. Commun. 2, 46–57. doi:10.1016/j.meteno.
Ebrahim, A., Lerman, J.A., Palsson, B.O., Hyduke, D.R., 2013. COBRApy: constraints- 2015.06.001.
based reconstruction and analysis for python. BMC Syst. Biol. 7. doi:10.1186/ Orth, Jeffrey D., Ines Thiele, B.Ø.P., 2010. What is flux balance? Nat. Biotechnol. 28,
1752- 0509- 7- 74. 245–248. doi:10.1038/nbt.1614.What.
Fleming, R.M.T., Thiele, I., Nasheuer, H.P., 2009. Quantitative assignment of reac- Planinc, A., Dejaegher, B., Vander Heyden, Y., Viaene, J., Van Praet, S., Rappez, F.,
tion directionality in constraint-based models of metabolism: Application to Es- Van Antwerpen, P., Delporte, C., 2017. Batch-to-batch N-glycosylation study of
cherichia coli. Biophys. Chem. 145, 47–56. doi:10.1016/j.bpc.20 09.08.0 07. infliximab, trastuzumab and bevacizumab, and stability study of bevacizumab.
Gutierrez, J.M., Feizi, A., Li, S., Kallehauge, T.B., Hefzi, H., Grav, L.M., Ley, D., Baycin Eur. J. Hosp. Pharm. 24. doi:10.1136/ejhpharm- 2016- 001022, 286 LP –292.
Hizal, D., Betenbaugh, M.J., Voldborg, B., Faustrup Kildegaard, H., Min Lee, G., Rui, B., Yi, Y., Shen, T., Zheng, M., Zhou, W., Du, H., Fan, Y., Wang, Y., Zhang, Z.,
Palsson, B.O., Nielsen, J., Lewis, N.E., 2020. Genome-scale reconstructions of the Xu, S., Liu, Z., Wen, H., Xie, X., 2015. Elementary flux mode analysis revealed
mammalian secretory pathway predict metabolic costs and limitations of pro- cyclization pathway as a powerful way for NADPH regeneration of central car-
tein secretion. Nat. Commun. 11, 1–10. doi:10.1038/s41467- 019- 13867- y. bon metabolism. PLoS One 10, 1–16. doi:10.1371/journal.pone.0129837.
Hefzi, H., Ang, K.S., Hanscho, M., Bordbar, A., Ruckerbauer, D., Lakshmanan, M., Ryder, A.G., 2018. Cell culture media analysis using rapid spectroscopic methods.
Orellana, C.A., Baycin-Hizal, D., Huang, Y., Ley, D., Martinez, V.S., Kyriakopou- Curr. Opin. Chem. Eng. 22, 11–17. doi:10.1016/j.coche.2018.08.008.
los, S., Jiménez, N.E., Zielinski, D.C., Quek, L.E., Wulff, T., Arnsdorf, J., Li, S., Salvy, P., Hatzimanikatis, V., 2020. The ETFL formulation allows multi-omics inte-
Lee, J.S., Paglia, G., Loira, N., Spahn, P.N., Pedersen, L.E., Gutierrez, J.M., King, Z.A., gration in thermodynamics-compliant metabolism and expression models. Nat.
Lund, A.M., Nagarajan, H., Thomas, A., Abdel-Haleem, A.M., Zanghellini, J., Kilde- Commun. 11. doi:10.1038/s41467- 019- 13818- 7.
gaard, H.F., Voldborg, B.G., Gerdtzen, Z.P., Betenbaugh, M.J., Palsson, B.O., An- Schellenberger, J., Palsson, B., 2009. Use of randomized sampling for analysis of
dersen, M.R., Nielsen, L.K., Borth, N., Lee, D.Y., Lewis, N.E., 2016. A consensus metabolic networks. J. Biol. Chem. 284, 5457–5461. doi:10.1074/jbc.R80 0 04820 0.
genome-scale reconstruction of chinese hamster ovary cell metabolism. Cell Schiestl, M., Stangler, T., Torella, C., Čepeljnik, T., Toll, H., Grau, R., 2011. Acceptable
Syst 3, 434–443. doi:10.1016/j.cels.2016.10.020, .e8.. changes in quality attributes of glycosylated biopharmaceuticals. Nat. Biotech-
Henry, C.S., Broadbelt, L.J., Hatzimanikatis, V., 2007. Thermodynamics-based nol. 29, 310–312. doi:10.1038/nbt.1839.
metabolic flux analysis. Biophys. J. 92, 1792–1805. doi:10.1529/biophysj.106. Schinn, S.-M., Morrison, C., Wei, W., Zhang, L., Lewis, N.E., 2020. A genome-s-
093138. cale metabolic network model synergizes with statistical learning to predict
Jimenez del Val, I., Nagy, J.M., Kontoravdi, C., 2011. A dynamic mathematical model amino acid concentrations in Chinese hamster ovary cell cultures. bioRxiv
for monoclonal antibody N-linked glycosylation and nucleotide sugar donor 2020.09.02.279687.
transport within a maturing Golgi apparatus. Biotechnol. Prog. 27, 1730–1743. Smallbone, K., Simeonidis, E., 2009. Flux balance analysis: A geometric perspective.
doi:10.1002/btpr.688. J. Theor. Biol. 258, 311–315. doi:10.1016/j.jtbi.2009.01.027.
Kim, M.K., Lane, A., Kelley, J.J., Lun, D.S., 2016. E-Flux2 and sPOT: Validated methods Sokolov, M., Morbidelli, M., Butté, A., Souquet, J., Broly, H., 2018. Sequential mul-
for inferring intracellular metabolic flux distributions from transcriptomic data. tivariate cell culture modeling at multiple scales supports systematic shaping
PLoS One 11, 1–22. doi:10.1371/journal.pone.0157101. of a monoclonal antibody toward a quality target. Biotechnol. J. 13, e1700461.
Kotidis, P., Demis, P., Goey, C.H., Correa, E., McIntosh, C., Trepekli, S., Shah, N., doi:10.10 02/biot.20170 0461.
Klymenko, O.V, Kontoravdi, C., 2019a. Constrained global sensitivity analysis Stöckmann, H., Adamczyk, B., Hayes, J., Rudd, P.M., 2013. Automated, high-
for bioprocess design space identification. Comput. Chem. Eng. 125, 558–568. throughput igg-antibody glycoprofiling platform. Anal. Chem. 85, 8841–8849.
doi:10.1016/j.compchemeng.2019.01.022. doi:10.1021/ac402068r.
Kotidis, P., Jedrzejewski, P., Sou, S.N., Sellick, C., Polizzi, K., del Val, I.J., Kontoravdi, C., Yang, O., Ierapetritou, M., 2021. mAb production modeling and design space evalu-
2019b. Model-based optimization of antibody galactosylation in CHO cell cul- ation including glycosylation process. Process doi:10.3390/pr9020324.
ture. Biotechnol. Bioeng. 116, 1612–1626. doi:10.1002/bit.26960. Yeo, H.C., Hong, J., Lakshmanan, M., Lee, D.Y., 2020. Enzyme capacity-based genome
Kotidis, P., Kontoravdi, C., 2020. Harnessing the potential of artificial neural net- scale modelling of CHO cells. Metab. Eng. 60, 138–147. doi:10.1016/j.ymben.
works for predicting protein glycosylation. Metab. Eng. Commun. 10, e00131. 2020.04.005.
doi:10.1016/j.mec.2020.e00131. Zavala-Ortiz, D.A., Ebel, B., Li, M.-Y., Barradas-Dermitz, D.M., Hayward-Jones, P.M.,
Kyriakopoulos, S., 2014. Amino acid metabolism in Chinese hamster ovary cell cul- Aguilar-Uscanga, M.G., Marc, A., Guedon, E., 2020. Interest of locally weighted
ture 247. regression to overcome nonlinear effects during in situ NIR monitoring of CHO
Kyriakopoulos, S., Ang, K.S., Lakshmanan, M., Huang, Z., Yoon, S., Gunawan, R., cell culture parameters and antibody glycosylation. Biotechnol. Prog. 36, e2924.
Lee, D.Y., 2018. Kinetic modeling of mammalian cell culture bioprocessing: the doi:10.1002/btpr.2924.
quest to advance biomanufacturing. Biotechnol. J. 13, 1–11. doi:10.1002/biot. Zürcher, P., Sokolov, M., Brühlmann, D., Ducommun, R., Stettler, M., Souquet, J.,
201700229. Jordan, M., Broly, H., Morbidelli, M., Butté, A., 2020. Cell culture process
Kyriakopoulos, S., Kontoravdi, C., 2014. A framework for the systematic design of metabolomics together with multivariate data analysis tools opens new routes
fed-batch strategies in mammalian cell culture. Biotechnol. Bioeng. 111, 2466– for bioprocess development and glycosylation prediction. Biotechnol. Prog. 36,
2476. doi:10.1002/bit.25319. e3012. doi:10.1002/btpr.3012.
Lewis, N.E., Hixson, K.K., Conrad, T.M., Lerman, J.A., Charusanti, P., Polpitiya, A.D.,
Adkins, J.N., Schramm, G., Purvine, S.O., Lopez-Ferrer, D., Weitz, K.K., Eils, R.,
König, R., Smith, R.D., Palsson, B., 2010. Omic data from evolved E. coli are

11