Многомасштабная разработка и оптимизация производства антител с помощью гибких сетей.

Introduction

Monoclonal antibodies are becoming increasingly important therapeutic agents in medicine. Their usefulness has been demonstrated in various applications, such as in the treatment of tumors [1] , [2] and autoimmune diseases [3] , [4] , as well as in biosensors [5] . For an antibody drug to be effective, it must have low immunogenicity and high efficiency in its intended application [6] . To achieve these properties, antibodies require multiple post-translational modifications, proper folding, and optimal cleavage [7] . Mammalian cells are the most popular and promising production host system to effectively meet these requirements [8] .

Here, we focus on Chinese Hamster Ovary (CHO) cells. CHO cells are capable of producing antibodies with glycosylation patterns similar to human serum antibodies. Although there may be some differences due to the absence of certain enzymes in CHO cells compared to the human system [9] , these antibodies closely resemble endogenous human antibodies. This similarity results in low immunogenicity and increased effectiveness, making these antibodies safer for therapeutic use. Moreover, CHO cells have a lower risk of viral contamination and spreading compared to other cell lines. For these reasons, CHO cells are the preferred system for producing large amounts of monoclonal antibody drugs while maintaining high product quality [8] .

Advances in technology have led to a 100-200-fold increase in antibody production in CHO cell cultures in the last 30 years [9] . However, large-scale cultures remain non-optimized due to their high cost relative to the total production of antibody drugs. The Design of Experiments (DoE) can optimize chemical production by identifying key process parameters, but it has limitations in capturing metabolic complexity and cellular responses over time. Since it relies on statistical correlations rather than mechanistic insights, in addition to being experimentally intensive, DoE can be less effective in predicting metabolic shifts [10] , [11] . Computational models of bioprocesses and cell metabolism offer a promising solution by guiding researchers in experimental design. By shifting from trial-and-error to rational bioprocess and cell line optimization, this approach can reduce investment costs and improve antibody production efficiency [12] , [13] , [14] , [15] .

The modeling of metabolic networks can be approached with several key methods. Constraint-based model methods, such as Flux Balance Analysis (FBA) [16] , use stoichiometry and optimization to predict cell behavior under steady-state assumptions, making them effective for genome-scale applications. Kinetic models [17] incorporate enzyme kinetics to capture transient time-dependent changes in metabolic fluxes, but often require extensive data. Metabolic Flux Analysis (MFA) [18] quantifies fluxes through measured uptake and secretion rates in simpler systems. Hybrid models, like dynamic FBA (dFBA) [19] , combine constraint-based models with dynamic elements to simulate changes over time. On the other hand, Petri nets [20] offer a graphical approach that enables the representation of regulatory and probabilistic elements.

Genome-scale metabolic models (GEMs) [21] are mathematical depictions of an organism's metabolism. They encode a comprehensive network of mass-balanced connections among the metabolites of the organism, based on gene-protein-reaction associations derived from genome annotation and experimental data. Although these models are not entirely accurate, many have been proven reliable enough to aid in decision making for wet lab experiments [22] , [23] . GEMs are inherently limited in their ability to capture multiscale variables, such as those arising from interactions between a bioreactor environment and a cell culture, because they focus primarily on intracellular metabolic reactions defined by stoichiometric constraints. GEMs operate at the cellular level, typically assuming a static, controlled environment without accounting for external variables like nutrient supply rate, pH shifts, product toxicity, biomass generation or other dynamic conditions that occur in a bioreactor [21] . These external factors can significantly impact cellular metabolism by altering reaction kinetics, product formation, and metabolic flux distributions.

Although the aforementioned approaches are valuable tools for modeling many real systems, they lack sufficient modeling power or analytical methods to address large multi-scale systems with uncertainties and nonlinear dynamics. Thus, alternative modeling frameworks and computational methods are necessary.

Flexible Nets (FNs) [24] is a modeling formalism that, in addition to representing the complete set of intracellular fluxes within a GEM, can incorporate easily bioreactor dynamics, including macroscopic variables such as dilution rate, biomass growth, and concentration of metabolites in the medium and in the tank. FNs offer a straightforward and intuitive graphical representation, and enable the modeling of regulation processes such as reaction activation or inhibition based on metabolite concentration [25] . Additionally, FNs can approximate nonlinear dynamics via piecewise linear functions and model uncertainty, enhancing their versatility as a modeling approach.

To demonstrate that FNs is a suitable formalism to model and optimize production by CHO cells, we model antibody production in a bioreactor operated in continuous culture mode, i.e., nutrients are continuously supplied to the tank at a specific dilution rate. Simultaneously, cells, cell products, and waste products are removed from the tank via the effluent, also regulated by the dilution rate. Although bioreactors often operate in batch and fed-batch modes for many biotechnological processes, there is a growing expectation that continuous cultures will become more prevalent [26] , [27] . Thus, our current modeling assumes a static exchange profile (given by measurements) of the cell with its environment, rather than dynamic changes over time. This allows us to investigate the trade-offs between growth and production using a novel methodology that has not been applied in this context before. Moreover, FN provides a flexible foundation for future work, where time-resolved exchange fluxes or other data can be seamlessly integrated, highlighting its advantages over dFBA for scalable and data-rich analyses.

In addition to validating and evaluating the model's performance, we also considered relevant aspects related to the antibody production optimization. One of these aspects is media design. Since not all amino acids in the medium have the same cost, and some may be present in excess beyond the requirements for antibody production, optimizing amino acid concentrations is essential.

In metabolic modeling, a phenotype represents a specific metabolic behavior defined by the flux distribution across reactions under given constraints. Different conditions, such as nutrient availability or gene knockouts, shape distinct phenotypes by altering metabolic fluxes. In our case, when producing a particular antibody, different phenotypes can be reached depending on the dilution rate and the biomass concentration in a bioreactor, each with distinct amino acid consumption patterns. Identifying the minimum amino acid requirements for each phenotype would help avoid resource waste and reduce production costs.

There are multiple approaches to develop an optimized medium, though many are based on empirical observations [28] , [29] . In this study, we utilize FNs to predict the minimum concentration of each amino acid in the medium needed to achieve a specific IgG concentration, or, in other words, a desired phenotype. While our product of interest is the antibody IgG, this approach can be applied to the production of other metabolites.

Materials and methods

Flexible nets

This section introduces Flexible Nets (FNs) [24] , a modeling formalism inspired by Petri Nets [30] , [31] . FNs can model effectively the interplay between the state and the processes of dynamic systems through two interconnected nets: the event net and the intensity net. While the event net defines how processes influence state variables, the intensity net determines how these variables regulate processes' rates [32] .

Both the event and the intensity nets are tripartite graphs, comprising places (associated with metabolites) depicted as circles (set P ), transitions (linked with reactions) depicted as rectangles (set T ), and handlers depicted as dots. There are two types of handlers: event handlers (set V ), which in the context of metabolic networks, trigger changes in metabolite concentrations when a reaction occurs, and intensity handlers (set S ) which determine the influence of metabolite concentrations on reaction rates.

Let us consider a simple system to exemplify some modeling capabilities of FNs. Hemoglobin is a tetrameric protein responsible for oxygen transport in the blood, consisting of two alpha and two beta subunits. It can bind a variable number of oxygen molecules, ranging from 1 to 4, depending on oxygen availability. Each subunit can independently bind an oxygen molecule. Although hemoglobin is most efficient when fully saturated with four oxygen molecules, it remains functional across intermediate states, allowing it to release oxygen to tissues based on their specific needs [33] .

If one molecule of hemoglobin binds two molecules of oxygen, the process results in two molecules of oxygen bound to one molecule of hemoglobin. Such a reaction can be represented as: R 1 : M a + 2 M b → M c where R 1 has two reactants, M a (hemoglobin) and M b (oxygen), and one product, M c (oxygen-bound hemoglobin).

The event net in Fig. 1 a, which models reaction R 1 , has three places, P = { p a , p b , p c } which are associated with M a , M b and M c and one event handler V = { v 1 } . The labels of the arcs, a , b , and c , represent the number of hemoglobin molecules consumed, the number of oxygen molecules consumed, and the number of oxygen-bound hemoglobin molecules produced by reaction R1, respectively. Since one hemoglobin molecule binds to two oxygen molecules, b must be twice a , i.e. b = 2 a . Additionally, because one oxygen-bound hemoglobin molecule is produced for each hemoglobin molecule, c must be equal to a , i.e. c = a . This way, the equations v 1 : b = 2 a ; c = a associated with event handler v 1 define the relationships between the number of molecules consumed and produced by a reaction. Fig. 1 (a) Event net modeling the stoichiometry of reaction R 1 : A  + 2 B  →  C . (b) Event net modeling the reaction R 1 : A  + 2 B  →  C with rate λ 0 [ t 1 ] constrained in the interval [ lb , ub ]. (c) Event net modeling the reaction R 1 : A  +  nB  →  C where n is the uncertain stoichiometric coefficient of p b . (d) Flexible Net modeling the reaction R 1 : A  +  nB  →  C with rate equal to λ 0 [ t 1 ] plus 3 times the concentration of M a . Here, the event net is represented by black arcs and edges, and the intensity net with blue arcs and edges. Fig. 1

(a) Event net modeling the stoichiometry of reaction R 1 : A  + 2 B  →  C . (b) Event net modeling the reaction R 1 : A  + 2 B  →  C with rate λ 0 [ t 1 ] constrained in the interval [ lb , ub ]. (c) Event net modeling the reaction R 1 : A  +  nB  →  C where n is the uncertain stoichiometric coefficient of p b . (d) Flexible Net modeling the reaction R 1 : A  +  nB  →  C with rate equal to λ 0 [ t 1 ] plus 3 times the concentration of M a . Here, the event net is represented by black arcs and edges, and the intensity net with blue arcs and edges.

Notice that the net in Fig. 1 a does not contain any dynamic information, i.e. the reaction rate is undefined. Such a rate can be modeled by adding the transition t 1 , see Fig. 1 b. Transition t 1 is connected to v 1 through an edge with label v which represents the number of occurrences of the reaction. Thus, the equations a = c = v and b = 2 v associated with the event handler v 1 imply that one occurrence of R 1 ( v = 1 ) consumes one molecule of M a ( a = 1 ), two molecules of M b ( b = 2 ) and produces one molecule of M c ( c = 1 ). The reaction rate of R 1 is given by the variable λ 0 [ t 1 ] which is associated with t 1 . In this particular example, λ 0 [ t 1 ] can be any real value in the interval [ l b , u b ] (see the inequalities associated with t 1 ).

The net in Fig. 1 c accounts for the fact that the number of oxygen molecules that can bind one molecule of hemoglobin is not a fixed value, but lies in the interval [ 1 , 4 ] . In this case, the reaction can be expressed as: (1) R 1 : M a + n M b → M c where n is an uncertain stoichiometric coefficient, that is, n units of M b in the interval [ 1 , 4 ] are consumed per each unit of M a to produce one unit of M c . Thus, the number of molecules of M b consumed per each occurrence of R 1 is an uncertain value in the interval [ 1 , 4 ] . In terms of fluxes (or rates), Equation (1) implies that the flux at which M b is consumed is n times the flux at which M a is consumed, and the flux at which M a is consumed is equal to the flux at which M c is produced.

Intensity nets can be used to model how reaction rates depend on the concentration of metabolites. In Fig. 1 d, the intensity net contains an intensity handler s 1 , an edge from p a to s 1 , and an arc from s 1 to t 1 . The equation t = 3 r associated with s 1 means that the rate, t , produced by M a in R 1 is equal to 3 times the concentration of M a , r . Thus, R 1 is a first order reaction. This rate is added to λ 0 [ t 1 ] , which is in [ l b , u b ] , to obtain the final reaction rate λ [ t ] . Similarly to event handlers, inequalities can be associated with intensity handlers to model uncertainties, e.g. 2 r ≤ t ≤ 4 r would mean that the rate produced is between 2 and 4 times the concentration of M a .

The reaction rate of R 1 in Fig. 1 d is linear with respect to the concentration of M a . In order to account for nonlinear dynamics, piecewise linear functions can be associated with the intensity handlers. For instance, if the piecewise linear function: (2) s 1 : { t = 3 r  if  r ≤ 5 t = 15  if  r > 5 is associated with s 1 , then the rate produced by M a will be equal to 3 times its concentration if such a concentration is equal or lower than a given threshold 5, otherwise the rate produced will not be proportional to the concentration but constant and equal to 15.

In most GEMs, reactions are described in terms of their stoichiometry (without uncertain coefficients) and flux bounds that are independent on metabolite concentrations. Thus, such reactions can be modeled straightforwardly through event nets similar to the one in Fig. 1 c.

Chinese Hamster Ovary cell model (iCHOv1)

iCHOv1 features

The iCHOv1 model from BiGG Models [34] is a comprehensive representation of the metabolic network of Chinese Hamster Ovary (CHO) cells, widely used in biopharmaceutical production [13] . It accounts for 4456 metabolites, 6663 reactions, and 1766 genes, and integrates extensive biochemical data, covering central carbon metabolism, amino acid and nucleotide biosynthesis, lipid metabolism, and energy production. The model includes Gene-Protein-Reaction (GPR) associations, linking specific genes to the enzymes and metabolic reactions they encode. It also has an artificial biomass composition reaction, precise reaction stoichiometry for mass and energy balance, and accounts for the following cellular compartments: cytosol, mitochondria, endoplasmic reticulum, and extracellular space.

IgG antibody synthesis reaction

Immunoglobulin G (IgG) is a versatile and abundant antibody crucial for immune response. Structurally, IgG is Y-shaped, with two heavy chains (Hc) and two light chains (Lc), and has two antigen-binding sites at the Y-tips for recognizing pathogens [35] , [36] .

The iCHOv1 model does not contain the reactions required to produce IgG. To enable antibody synthesis, we added these reactions to the model. Given that IgG is a protein, amino acids are sequentially integrated in the polypeptide chain until one IgG subunit, I g G s u b u n i t , is completed (see reaction R 1 below), as defined in [37] . Afterwards, two antibody subunits must assemble to produce the functional IgG (see R 2 ). To facilitate antibody secretion, we also included an exchange reaction (see R 3 ). The complete list of added reactions is as follows: R 1 : 39 A l a + 14 A r g + . . . + 1324 G T P + 1326 H 2 O → I g G s u b u n i t + . . . + 1326 H R 2 : 2 I g G s u b u n i t → I g G R 3 : I g G → ∅

An FN graphical representation of these reactions is shown in Fig. 2 . The labels a, b, c, etc. model the stoichiometry of each reaction. For instance, a = 39 r in reaction R 1 means that 39 moles of the amino acid Ala react to produce a mole of the I g G s u b u n i t . The complete reaction containing all the stoichiometric coefficients is reported in Table 3 (see Annex). Fig. 2 FN modeling the reactions that produce IgG in the iCHOv1 model. Note that in reaction R 1 , only some of the involved amino acids are represented (see Table 3 for the comprehensive list of stoichiometric coefficients). Fig. 2

FN modeling the reactions that produce IgG in the iCHOv1 model. Note that in reaction R 1 , only some of the involved amino acids are represented (see Table 3 for the comprehensive list of stoichiometric coefficients).

Modeling a continuous culture with FNs

The bioreactor used for the continuous culture consists of three compartments: the reservoir, the tank and the effluent. Fig. 3 sketches an FN model that integrates the microscopic variables of the model, i.e., the fluxes of the reactions in the metabolic network of a CHO cell, with the macroscopic variables of the bioreactor, see Table 1 . Fig. 3 Bioreactor consisting of three compartments: reservoir, tank and effluent. The blue rectangle “Cell” at the bottom contains the whole metabolic network of iCHOv1. Fig. 3 Table 1 Macroscopic variables involved in the bioreactor dynamics. Table 1 Macroscopic variables Units Description D h −1 Dilution rate s m mM Nutrient concentration in medium s t mM Nutrient concentration in tank X gDWL −1 Biomass concentration a t mM Antibody concentration in tank

Bioreactor consisting of three compartments: reservoir, tank and effluent. The blue rectangle “Cell” at the bottom contains the whole metabolic network of iCHOv1.

Table 1 Macroscopic variables involved in the bioreactor dynamics. Table 1 Macroscopic variables Units Description D h −1 Dilution rate s m mM Nutrient concentration in medium s t mM Nutrient concentration in tank X gDWL −1 Biomass concentration a t mM Antibody concentration in tank

Macroscopic variables involved in the bioreactor dynamics.

For a better understanding of Fig. 3 , the step-by-step model development process is discussed below. 1. Add the necessary reactions to the GEM to produce IgG: The GEM iCHOv1 is extended to accommodate the reactions detailed in Section 2.2.2 for IgG production. The cell compartment in Fig. 3 contains the whole iCHOv1 GEM together with these reactions (for clarity, only reaction R 2 for antibody production is depicted). Before discussing the next steps, it must be noted that the evolution of the concentration of a given nutrient, s , in a bioreactor is determined by: (3) d s t d t = ( s m − s t ) D − u s X where s m is the concentration of nutrient s in the medium, s t is the concentration of the same nutrient in the tank, D is the dilution rate, u s is the uptake rate of the nutrient by the cell, and X is the biomass concentration (see Table 1 for a description of the macroscopic variables). Note that the terms in Equation (3) have a negative contribution if the metabolite is removed from the tank, and a positive contribution if the metabolite goes from a certain compartment into the tank. Steps 2, 3 and 4 explain how each of these terms is modeled by FNs. 2. Model nutrients' fluxes from the medium to the tank ( D s m ): A given nutrient s , with concentration s m in the medium, is supplied from the reservoir to the tank at a dilution rate D , i.e. the flux of s from the reservoir to the tank is D s m (see Equation (3) ). This flux is modeled by the transition t s i n in Fig. 4 (a) whose intensity is λ 0 [ t s i n ] = D s m . Note that, although only one nutrient is depicted in Fig. 3 , the model accounts for all the nutrients in Table 5 . Fig. 4 (a) Subnet modeling the relationship between the reservoir, the nutrient in the tank and the effluent in the bioreactor. (b) Subnet modeling nutrient uptake. Fig. 4 3. Model the nutrient uptake from the tank ( − u s X ): Once a nutrient s is in the tank (such a nutrient is modeled by place S ), it can be either transported into the cell or expelled to the effluent. The former process corresponds with the term − u s X and is modeled by Fig. 4 (b). Such transport of s from the tank to a single cell is modeled by transition t s , which matches with the exchange reaction of s in the iCHOv1 model. This exchange reaction will have a metabolic flux u s indicating the rate of consumption of metabolite s by the cell. The key element that gathers the microscopic (intracellular fluxes of iCHOv1) and macroscopic (bioreactor dynamics) variables and allows for the development of a multi-scale model is the intensity handler h s . The equation associated with h s scales the exchange flux of metabolite s by the biomass concentration X , i.e. s t = u s X (see Equation (3) ), to obtain the overall flux of s , λ [ t s t ] , consumed by the cell culture. 4. Model nutrients' fluxes from the bioreactor to the effluent ( − D s t ): As mentioned in the previous step, a certain amount of the metabolite modeled by place S will be carried by the effluent, preventing it from being consumed by cells. This process is depicted in Fig. 4 (a), specifically by t s o u t and the elements connected to it ( v s o u t and s s o u t ). The rate at which this process occurs is defined by the intensity handler s s o u t (see Equation (3) ). 5. Model the biomass in the bioreactor: The biomass evolution in the bioreactor is given by: (4) d X d t = ( μ − D ) X where μ is the specific growth rate of the cell, which coincides with the flux of the biomass reaction of the iCHOv1 model, D is the dilution rate, and X is the biomass concentration in the bioreactor. The first part of Equation (4) , μX , is modeled with FNs in the same way as in step 3, with a slight difference. In this case, the term is positive because the growth rate produces biomass in the tank. In this part of the model, the intensity handler h x acts as the interface between the microscopic and macroscopic parts of the model ( Fig. 5 (a)). Fig. 5 (a) Subnet modeling biomass production. (b) Subnet modeling the output flux of X . Fig. 5 The second part of the Equation (4) , − D X , models how cells are removed from the tank towards the effluent ( t x o u t ) at the rate indicated by Equation (4) (see Fig. 5 (b)). 6. Add the IgG antibody production: The concentration of antibodies in the bioreactor is determined by: (5) d a t d t = u a X − D a t Equation (5) is closely related to Equation (3) . The first part consists of u a X , where u a is the flux of the antibody exchange reaction from the modified iCHOv1 model ( R 3 in Fig. 2 ). This cellular-level reaction is connected to the macroscopic variables through the intensity handler h a ( Fig. 6 (a)). Note that the term is positive, since the antibody is produced and secreted from cell to the tank. Fig. 6 (a) Subnet modeling the production of antibodies. (b) Subnet modeling the flux of antibodies from the bioreactor to the effluent. Fig. 6 The second term of Equation (5) , D a t , is the rate at which the antibody leaves the tank through the effluent as captured in s a o u t ( Fig. 6 (b)).

Add the necessary reactions to the GEM to produce IgG: The GEM iCHOv1 is extended to accommodate the reactions detailed in Section 2.2.2 for IgG production. The cell compartment in Fig. 3 contains the whole iCHOv1 GEM together with these reactions (for clarity, only reaction R 2 for antibody production is depicted).

Before discussing the next steps, it must be noted that the evolution of the concentration of a given nutrient, s , in a bioreactor is determined by: (3) d s t d t = ( s m − s t ) D − u s X where s m is the concentration of nutrient s in the medium, s t is the concentration of the same nutrient in the tank, D is the dilution rate, u s is the uptake rate of the nutrient by the cell, and X is the biomass concentration (see Table 1 for a description of the macroscopic variables).

Note that the terms in Equation (3) have a negative contribution if the metabolite is removed from the tank, and a positive contribution if the metabolite goes from a certain compartment into the tank. Steps 2, 3 and 4 explain how each of these terms is modeled by FNs.

Model nutrients' fluxes from the medium to the tank ( D s m ): A given nutrient s , with concentration s m in the medium, is supplied from the reservoir to the tank at a dilution rate D , i.e. the flux of s from the reservoir to the tank is D s m (see Equation (3) ). This flux is modeled by the transition t s i n in Fig. 4 (a) whose intensity is λ 0 [ t s i n ] = D s m . Note that, although only one nutrient is depicted in Fig. 3 , the model accounts for all the nutrients in Table 5 . Fig. 4 (a) Subnet modeling the relationship between the reservoir, the nutrient in the tank and the effluent in the bioreactor. (b) Subnet modeling nutrient uptake. Fig. 4

(a) Subnet modeling the relationship between the reservoir, the nutrient in the tank and the effluent in the bioreactor. (b) Subnet modeling nutrient uptake.

Model the nutrient uptake from the tank ( − u s X ): Once a nutrient s is in the tank (such a nutrient is modeled by place S ), it can be either transported into the cell or expelled to the effluent. The former process corresponds with the term − u s X and is modeled by Fig. 4 (b). Such transport of s from the tank to a single cell is modeled by transition t s , which matches with the exchange reaction of s in the iCHOv1 model. This exchange reaction will have a metabolic flux u s indicating the rate of consumption of metabolite s by the cell. The key element that gathers the microscopic (intracellular fluxes of iCHOv1) and macroscopic (bioreactor dynamics) variables and allows for the development of a multi-scale model is the intensity handler h s . The equation associated with h s scales the exchange flux of metabolite s by the biomass concentration X , i.e. s t = u s X (see Equation (3) ), to obtain the overall flux of s , λ [ t s t ] , consumed by the cell culture.

Model nutrients' fluxes from the bioreactor to the effluent ( − D s t ): As mentioned in the previous step, a certain amount of the metabolite modeled by place S will be carried by the effluent, preventing it from being consumed by cells. This process is depicted in Fig. 4 (a), specifically by t s o u t and the elements connected to it ( v s o u t and s s o u t ). The rate at which this process occurs is defined by the intensity handler s s o u t (see Equation (3) ).

Model the biomass in the bioreactor: The biomass evolution in the bioreactor is given by: (4) d X d t = ( μ − D ) X where μ is the specific growth rate of the cell, which coincides with the flux of the biomass reaction of the iCHOv1 model, D is the dilution rate, and X is the biomass concentration in the bioreactor.

The first part of Equation (4) , μX , is modeled with FNs in the same way as in step 3, with a slight difference. In this case, the term is positive because the growth rate produces biomass in the tank. In this part of the model, the intensity handler h x acts as the interface between the microscopic and macroscopic parts of the model ( Fig. 5 (a)). Fig. 5 (a) Subnet modeling biomass production. (b) Subnet modeling the output flux of X . Fig. 5

(a) Subnet modeling biomass production. (b) Subnet modeling the output flux of X .

The second part of the Equation (4) , − D X , models how cells are removed from the tank towards the effluent ( t x o u t ) at the rate indicated by Equation (4) (see Fig. 5 (b)).

Add the IgG antibody production: The concentration of antibodies in the bioreactor is determined by: (5) d a t d t = u a X − D a t

Equation (5) is closely related to Equation (3) . The first part consists of u a X , where u a is the flux of the antibody exchange reaction from the modified iCHOv1 model ( R 3 in Fig. 2 ). This cellular-level reaction is connected to the macroscopic variables through the intensity handler h a ( Fig. 6 (a)). Note that the term is positive, since the antibody is produced and secreted from cell to the tank. Fig. 6 (a) Subnet modeling the production of antibodies. (b) Subnet modeling the flux of antibodies from the bioreactor to the effluent. Fig. 6

(a) Subnet modeling the production of antibodies. (b) Subnet modeling the flux of antibodies from the bioreactor to the effluent.

The second term of Equation (5) , D a t , is the rate at which the antibody leaves the tank through the effluent as captured in s a o u t ( Fig. 6 (b)).

Notice that Equations (3) , (4) and (5) are all equal to zero since we assume the steady-state condition, typically reached in the continuous operational mode of a bioreactor. Thus, in order to analyze and optimize the resulting FN, a set of mathematical constraints that necessarily hold at steady-state will be derived [32] . These constraints, along with an appropriate objective function, e.g. the maximization of the flux of the antibody exchange reaction, will form a programming problem that can be used to make various predictions.

Implementation of the optimization procedure

The reactions related to the synthesis of antibodies were added to the iCHOv1 that was obtained from the BiGG Models database, ( http://bigg.ucsd.edu/models/iCHOv1 ). This modified model was converted into a Flexible Net (FN) by utilizing the open-source Python package fnyzer (Flexible Nets analYZER). FNs allowed us to model a multiscale system that integrates the metabolic network of CHO cells and the dynamics of the bioreactor. Then, the addition of flux bounds, steady state conditions and an objective function (in this case, maximizing the IgG antibodies production) resulted in a linear programming problem (LPP) whose solution represented a theoretical maximum for IgG production. The optimization software used by fnyzer to solve the resulting LPP was CPLEX [38] .

The modifications in the iCHOv1 model were performed using the 0.27.0 version of COBRApy Python package [39] . All the simulations and the analyses were conducted by fnyzer 1.3.9 [40] in combination with the CPLEX solver. The developed scripts are available at: https://github.com/jlazaroibanezz/IgG_CHO_FNs .

Simulations were run on a workstation equipped with an 11th Gen Intel(R) Core(TM) i7-1165G7 processor (2.80 GHz, 4 cores) and 15 GB of RAM. The system was running Ubuntu 22.04.3 with Linux kernel version 6.8.0. One optimization of the system takes less than 30 seconds.

Results

Model validation

To evaluate our model's predictive capacity, we compared its predicted IgG production with experimentally reported values from [15] . This comparison was based on two datasets, each containing experimental uptake rates of medium metabolites (in m m o l
g D W − 1
h − 1 ) and the growth rate of the cell culture ( h − 1 ).

The first dataset, from [41] , represents measurements taken during the late exponential phase of cell culture growth (Late Exponential phase dataset). The second dataset, from [12] , includes measurements from a different cell culture during its exponential phase (HP dataset). Both datasets correspond with a fed-batch operation mode, in which nutrients are added in a controlled manner following an initial growth phase (batch phase). For additional details on these datasets, see Annex ( Table 4 ).

The equations describing the bioreactor dynamics (Equations (3) , (4) and (5) ) specifically model the dynamics of a continuous bioreactor. We used the data of the exponential phase of the fed-batch datasets as a good approximation for our model parameters. Our continuous culture FN model aligns with FBA conditions, as both maintain constant metabolite concentrations and fluxes, ensuring a balanced metabolic network in steady-state. The growth rate in a chemostat, constrained by the dilution rate, mirrors the fixed exchange fluxes in FBA, making our optimization process a relevant approach for studying metabolism under steady-state conditions [41] , [15] . Previous studies have used genome-scale metabolic models to analyze CHO cell metabolism in fed-batch cultures, applying FBA and FVA under steady-state assumptions [12] , [41] . Similarly, we integrated experimental uptake rates from fed-batch operations into our FN model and compared predicted and experimental IgG production ( Table 2 ). Table 2 Antibody production in mmol g D W − 1 h −1 . The row “Model” reports the model predictions for the two experimental datasets assessed (Late Exponential phase and HP). The row “Experimental” contains the antibody production in [41] for the Late Exponential phase dataset, and in [12] for the HP dataset. Table 2 HP Late Exp phase Experimental 2.02 ⋅ 10 −5 2.44 ⋅ 10 −5 Model 2.04 ⋅ 10 −5 3.205 ⋅ 10 −5 Relative error 9.9 ⋅ 10 −3 0.3

Table 2 Antibody production in mmol g D W − 1 h −1 . The row “Model” reports the model predictions for the two experimental datasets assessed (Late Exponential phase and HP). The row “Experimental” contains the antibody production in [41] for the Late Exponential phase dataset, and in [12] for the HP dataset. Table 2 HP Late Exp phase Experimental 2.02 ⋅ 10 −5 2.44 ⋅ 10 −5 Model 2.04 ⋅ 10 −5 3.205 ⋅ 10 −5 Relative error 9.9 ⋅ 10 −3 0.3

Antibody production in mmol g D W − 1 h −1 . The row “Model” reports the model predictions for the two experimental datasets assessed (Late Exponential phase and HP). The row “Experimental” contains the antibody production in [41] for the Late Exponential phase dataset, and in [12] for the HP dataset.

The final FN-based model integrates the following constraints: 1) the medium formulation, 2) experimental upper bound flux constraints for the uptake rates of each amino acid in the cell culture retrieved from [41] , [12] , 3) the equations modeling the bioreactor dynamics, and 4) the mass balance and stoichiometric constraints inherent in iCHOv1. Once the model is transformed into an FN by means of the fnyzer Python package [40] , we impose the steady-state condition and set as objective function the maximization of the flux of the antibody synthesis reaction ( R 3 in Fig. 2 ). This results in a null flux for the growth rate. Conversely, when the objective function is the maximization of the growth rate, the flux of the antibody synthesis reaction becomes zero. This implies that the available medium resources are entirely allocated either to antibody synthesis or to growth, depending on the selected objective function. Constraining the growth rate with the experimental values provided in the “Experimental growth rate” row from Table 4
[12] , [41] overcomes the problem.

By imposing the steady-state condition and maximizing the IgG synthesis under the aforementioned constraints, we predicted the values reported in Table 2 (row ‘Model’). These predicted values were then compared to the experimentally measured IgG synthesis reaction fluxes (Experimental). The relative error was computed as follows: (6) R e l a t i v e E r r o r = | M o d e l − E x p e r i m e n t a l | E x p e r i m e n t a l

Our model tends to overestimate the IgG production in all cases. However, the prediction improves when the model is enriched with the experimental data from the HP dataset. Although the relative error for the Late Exponential phase dataset is higher, it is not significant, and the experimental and predicted IgG production values remain within the same order of magnitude (see Table 2 ).

IgG production optimization

After evaluating the predictive capacity of our model, we further investigated how the parameters of the macroscopic variables could influence the production of IgG in the cell culture. We selected typical intervals for the dilution rate and cell density observed in CHO cell cultures and focused on maximizing the IgG synthesis reaction flux while fixing the growth rate to the sampled dilution rates, as described in Equation (4) (equal to zero in steady-state).

The heatmap in Fig. 7 shows the flux of IgG synthesis, measured in μ M h − 1 , as a function of dilution rate ( h − 1 ) and biomass concentration (10 6 cells m l − 1 ) using data from the HP dataset. The peak IgG production is observed in the central region of the heatmap, with the highest flux value of 0.194 μ M h − 1 occurring at a dilution rate of 0.0166 h − 1 and a biomass concentration of 10.624 ⋅ 10 6 cells m l − 1 . Blank rows and columns represent regions where no feasible optimization solution exists. Fig. 7 Heatmap showing the IgG synthesis reaction flux ( μM h −1 ) as a function of the dilution rate ( h −1 ) and biomass (10 6 cells ml −1 ). The abbreviation of the limiting amino acid is indicated for each biomass and dilution rate sampled. The implemented data corresponds with the HP dataset. Fig. 7

Heatmap showing the IgG synthesis reaction flux ( μM h −1 ) as a function of the dilution rate ( h −1 ) and biomass (10 6 cells ml −1 ). The abbreviation of the limiting amino acid is indicated for each biomass and dilution rate sampled. The implemented data corresponds with the HP dataset.

In comparison, Fig. 8 shows lower overall IgG production rates when using parameters derived from the Late Exponential phase dataset. This suggests that the Late Exponential phase data impose stricter constraints, which limit the IgG reaction flux more significantly than the parameters derived from the HP dataset. Fig. 8 Heatmap with the IgG synthesis reaction flux ( μM h −1 ) as a function of the dilution rate ( h −1 ) and biomass (10 6 cells ml −1 ). The abbreviation of the limiting amino acid is indicated for each biomass and dilution rate sampled. Here we implemented the data from the Late Exponential dataset. Fig. 8

Heatmap with the IgG synthesis reaction flux ( μM h −1 ) as a function of the dilution rate ( h −1 ) and biomass (10 6 cells ml −1 ). The abbreviation of the limiting amino acid is indicated for each biomass and dilution rate sampled. Here we implemented the data from the Late Exponential dataset.

An important consideration when working with culture media is identifying the limiting metabolite. To determine which metabolite is limiting in each scenario, we implemented the following strategy: 1. Fix the IgG synthesis reaction flux to the value obtained when maximized. 2. Maximize and minimize the fluxes of the exchange reactions for each amino acid one at a time. 3. Compute the difference between the maximum and the minimum exchange fluxes for each amino acid. 4. The amino acid whose difference between maximum and minimum fluxes for its exchange reaction is zero will be the limiting metabolite.

Fix the IgG synthesis reaction flux to the value obtained when maximized.

Maximize and minimize the fluxes of the exchange reactions for each amino acid one at a time.

Compute the difference between the maximum and the minimum exchange fluxes for each amino acid.

The amino acid whose difference between maximum and minimum fluxes for its exchange reaction is zero will be the limiting metabolite.

In other words, a metabolite is said to be limiting when the maximum and minimum possible fluxes of its exchange reaction are equal. Consequently, even a subtle change in the flux of the exchange reaction of the limiting metabolite can significantly impact IgG production.

The results for the HP dataset ( Fig. 7 ) show that, at low dilution rates and cell densities, the limiting metabolite varies between tyrosine (Tyr) and phenylalanine (Phe). However, as the dilution rate and cell density increase, lysine (Lys) becomes the limiting metabolite.

For the Late Exponential phase dataset ( Fig. 8 ), threonine (Thr), tyrosine (Tyr), and valine (Val) are limiting at low dilution rates and cell densities, but lysine (Lys) becomes limiting as these parameters increase.

Assessing medium uncertainty

The composition of a commercial medium may vary between batches. This can occur due to raw material variability or differences in the manufacturing process. Such batch-to-batch differences can impact cell growth and experimental outcomes, especially for sensitive cultures. To assess such an impact, the concentration of each metabolite in the medium was assumed to have a given degree of uncertainty, while the overall concentration of the medium was kept constant.

In particular, to model uncertainty in the medium, the concentration of each nutrient in the medium was not fixed to the particular value c provided in Table 5 , but was allowed to take any value in the interval [ c ⋅ ( 1 − u ) , c ⋅ ( 1 + u ) ] where u is the degree of uncertainty. Notice that this can be easily accommodated in the FN model by replacing equalities by inequalities.

By applying the steady-state condition and optimizing the Linear Programming Problem (LPP) with an objective function, specifically the IgG synthesis reaction, the theoretical maximum IgG production was computed for values of uncertainty ranging from 0% to 25%, see Fig. 9 . Fig. 9 IgG production as a function of the percentage of uncertainty incorporated in the medium nutrients. Fig. 9

IgG production as a function of the percentage of uncertainty incorporated in the medium nutrients.

It can be seen that for low degrees of uncertainty, the production increases linearly. This is because higher uncertainty allows for more flexibility in optimizing nutrient concentrations, leading to higher yields. However, when the uncertainty reaches 10%, the increase starts to slow down, and from 15%, IgG production remains constant. This behavior is due to the uptake rate limitations of individual nutrients by the cells.

Medium optimization

Another application of the developed model is medium optimization. We determined which nutrients are present in excess and identified the minimum concentrations required to sustain a certain level of IgG production. In particular, we focused on the following optimization objectives: 1) economic medium minimization, and 2) minimization of the number of nutrients in the medium.

Economic medium minimization

For the economic medium minimization, we selected the dilution rate and biomass concentration at which the maximum IgG synthesis reaction flux is achieved, as shown in Fig. 7 , Fig. 8 . We fixed this value and minimized the weighted sum of the intensities of the incoming transitions to the tank ( t s i n in Fig. 3 ), with each nutrient weighted according to its market cost as listed in the US Pharmacopeia [42] . This is expressed as: (7) min ∑ s ∈ N u t k s ⋅ λ [ t s i n ] where Nut is the set of nutrients in the medium (see Table 5 ), k s is the weight based on the market cost of nutrient s , and λ [ t s i n ] is the flux entering the tank from the reservoir. This optimization calculates the required concentration of each nutrient entering the tank from the reservoir compartment. We conducted this analysis both for the HP and the Late Exponential phase datasets, each corresponding to the maximum flux of the IgG synthesis reaction in each dataset.

The optimization for the HP dataset produced the results shown in Fig. 10 . The CD CHO medium (see Table 5 ), typically used for CHO cell cultures, appears to contain all nutrients in excess except for methionine (Met) and tyrosine (Tyr). According to our integrated FN model, adjusting nutrient concentrations to the levels represented by the orange bars would not affect the IgG production levels while reducing the overall medium cost. Fig. 10 Bar plot representing the concentrations of each nutrient in the original CD CHO medium (blue bars) and in the optimized medium according to the FN-based model enriched with experimental data from the HP dataset (orange bars). Fig. 10

Bar plot representing the concentrations of each nutrient in the original CD CHO medium (blue bars) and in the optimized medium according to the FN-based model enriched with experimental data from the HP dataset (orange bars).

We performed the same optimization generating Fig. 11 , with the FN model including the parameters from the Late Exponential phase dataset. Fig. 11 Bar plot that represents the concentrations of each nutrient in the original CD CHO medium (blue bars) and in the optimized medium according to the FN-based model enriched with experimental data from the Late Exponential phase dataset (orange bars). Fig. 11

Bar plot that represents the concentrations of each nutrient in the original CD CHO medium (blue bars) and in the optimized medium according to the FN-based model enriched with experimental data from the Late Exponential phase dataset (orange bars).

Minimization of the number of nutrients in the medium

This section aims to minimize the number of nutrients in the medium rather than its economic cost. To achieve this goal, we first set a null flux from the reservoir to the tank, i.e. λ [ t s i n ] = 0 , to each nutrient individually, this is equivalent to removing the nutrient from the medium, while keeping the biomass and dilution rate fixed at the values that produce the maximum IgG synthesis flux. Then, antibody production is maximized subject to the previous constraints in each dataset. Nutrients in the media were classified depending on their effect in IgG production: • Indispensable nutrients : Nutrients whose removal from the medium results in a null antibody production. • Modulatory nutrients : Nutrients whose removal reduces antibody production but does not stop it completely. • Redundant nutrients : Nutrients whose removal causes no change in antibody production.

Indispensable nutrients : Nutrients whose removal from the medium results in a null antibody production.

Modulatory nutrients : Nutrients whose removal reduces antibody production but does not stop it completely.

Redundant nutrients : Nutrients whose removal causes no change in antibody production.

The simulations identified nine indispensable nutrients: histidine (His), isoleucine (Ile), leucine (Leu), lysine (Lys), methionine (Met), phenylalanine (Phe), threonine (Thr), tryptophan (Trp), and valine (Val); two modulatory nutrients: cysteine (Cys) and tyrosine (Tyr); and eight redundant nutrients: arginine (Arg), asparagine (Asn), aspartate (Asp), glucose (Glc), glutamine (Gln), glutamate (Glu), proline (Pro), and serine (Ser).

Additionally, we minimized the sum of all λ [ t s i n ] values from the reservoir to the tank after having fixed the antibody production to the maximum obtained in Section 3.2 : (8) min ∑ s ∈ N u t λ [ t s i n ]

This optimization revealed that the maximum antibody production can still be achieved when the concentrations of all the redundant nutrients are simultaneously set to zero. Consequently, the minimum number of nutrients required to sustain the maximum IgG production along with the associated biomass is 11, i.e., the number of the indispensable nutrients plus the number of the modulatory nutrients.

Optimizing the model with Equation (8) as objective function allowed us to compute the minimum sum of fluxes of each nutrient necessary to achieve the fixed IgG production. The obtained fluxes were the same as in the previous section for the Late Exponential phase dataset, see Fig. 11 . However, there were slight differences between the economic and number of nutrients minimizations when using the HP dataset, see Fig. 12 . A comparison of Fig. 10 , Fig. 12 shows that the economically minimized medium contains less cysteine and phenylalanine but more methionine and tyrosine compared to the medium with minimized number of nutrients. This difference is likely due to tyrosine's relatively lower cost compared to phenylalanine. Both amino acids are among the three aromatic amino acids, sharing similar structures and roles in cellular metabolism and IgG production. As a result, substituting a specific quantity of phenylalanine with tyrosine can help reduce overall costs. Fig. 12 Bar plot showing the concentrations of each nutrient in the original CD CHO medium (blue bars) and in the optimized medium according to the FN-based model enriched with experimental data from the HP phase dataset (orange bars). In this case the minimization of the medium was performed focusing on the number of nutrients utilized. Fig. 12

Bar plot showing the concentrations of each nutrient in the original CD CHO medium (blue bars) and in the optimized medium according to the FN-based model enriched with experimental data from the HP phase dataset (orange bars). In this case the minimization of the medium was performed focusing on the number of nutrients utilized.

Discussion

The computational model predicting IgG production was validated using experimental data from [41] and [12] . This model incorporates constraints like media formulation, amino acid uptake rates, bioreactor dynamics, and stoichiometric balance. By imposing steady-state conditions and optimizing antibody synthesis, the model revealed a trade-off between cell growth and IgG production, reflecting biological resource allocation principles. Constraining with experimental growth rates resolved unrealistic predictions. The model generally overestimated IgG production but remained within the experimental range, performing best with the HP dataset.

Higher dilution rates can be reached in the HP simulations ( Fig. 7 ) compared to the ones in the Late Exponential phase ( Fig. 8 ) before having a washout. Washout happens when the dilution rate (the rate at which fresh medium is added and culture broth is removed) exceeds the specific growth rate of the cells. If the dilution rate is too high, cells are removed from the reactor faster than they can reproduce, leading to a decrease in biomass concentration.

Our model's assessment using both datasets revealed the critical influence of dilution rate and cell density on IgG production in CHO cell cultures, identifying optimal conditions at a certain dilution rate and biomass concentration. Higher IgG production rates were observed with the HP dataset, while the Late Exponential phase dataset imposed stricter constraints, and hence, the solution space was reduced in comparison to the one generated by the HP parameter values.

Identifying limiting metabolites (tyrosine, phenylalanine, threonine, valine, and notably lysine at higher rates and densities) highlighted the importance of specific amino acids in maximizing IgG synthesis. The limiting metabolite is dependent on the biomass, dilution rate, uptake rate, growth rate and IgG synthesis reaction flux. Lysine consistently emerged as a critical limiting factor as the dilution rate and biomass increased.

These insights can guide bioprocess optimization, emphasizing parameter fine-tuning and targeted nutrient supplementation to enhance IgG yields. By understanding constraints and metabolic demands, biotechnologists can develop more efficient and scalable processes for therapeutic antibody production, ultimately improving yield and reducing costs.

The study also focuses on optimizing medium composition for CHO cell cultures to maximize IgG production efficiently while minimizing costs. Utilizing computational modeling and experimental data integration, two distinct strategies were employed: economic optimization based on nutrient costs and reduction of the number of nutrients while maintaining productivity. Across different datasets (HP and Late Exponential), only 11 essential nutrients consistently supported maximum IgG synthesis, underscoring their critical roles in cellular metabolism and protein production. When a cell line is auxotrophic for a nutrient, this cell line is not able to grow in a medium lacking this element. Thus, CHO cells are auxotrophic for the 9 amino acids that are essential in mammalian cells: His, Ile, Leu, Lys, Met, Phe, Thr, Trp, and Val [43] , [44] , [45] .

After minimizing the medium, our model predicts minimum concentrations of these amino acids, but is always different from zero. Supplementing these nutrients in the medium is necessary, and the removal of one of these from the medium would cause no growth of CHO cells at all. Our model also predicts minimum concentrations for two additional amino acids, Tyr and Cys, which are not categorized as essential, but can act as such under specific growth conditions. Furthermore, multiple CHO cell lines have been proven to be Cys auxotrophic [46] . However, several studies suggest that many CHO cell lines are also auxotrophic for proline (Pro) [47] , [48] and arginine (Arg) [49] . This limitation of the iCHOv1 genome-scale model has been reported in [44] , [50] . If a cell-line specific model were used instead, it would predict a different optimized medium to accommodate the additional auxotrophies.

This approach not only highlights metabolic flexibility in substituting costly nutrients with cheaper alternatives but also ensures production efficiency. The uptake rate of a metabolite depends on several factors, including the concentration of metabolites in the culture medium, the metabolic demands of the cells, the expression and activity of specific transporters on the cell membrane, and the overall cell density [51] , [52] . That is why including kinetic information of the membrane transporters would be useful to improve our model.

By aligning the composition of the medium with economic considerations and bioprocessing demands, the study contributes to enhancing the sustainability and cost-effectiveness of biopharmaceutical production. Here, integrated approaches combining experimental data with FNs have proven useful for refining bioprocess strategies, offering valuable advancements in industrial biotechnology and biopharmaceutical development.

Conclusions

In this research, Flexible Nets (FNs) were exploited to integrate multiple types of biological constraints in a single model. FNs offer an adaptable modeling approach that integrates intracellular fluxes with macroscopic bioreactor dynamics, such as biomass growth and metabolite concentration. They can handle dynamic and nonlinear systems, include uncertainties, and allow for regulatory modeling based on metabolite levels. Standard methods like MFA, FBA, or dynamic FBA cannot fully capture the complexity of the system, as they focus only on intracellular fluxes. Our FN model integrates both intracellular fluxes (from FBA, MFA, etc.) and bioreactor dynamics into a single framework, creating a hybrid multi-scale approach. While it shares similarities with MFA and FBA by incorporating experimental uptake rates and solving an LPP, it uniquely includes macroscopic bioreactor variables, which standard methods cannot simulate.

The developed computational model for predicting IgG production, which was validated using experimental data from two studies, effectively incorporates constraints like media formulation, amino acid uptake rates and continuous culture dynamics as well as organisms' metabolic constraints. The model revealed a trade-off between cell growth and IgG production, which was mitigated by adjusting experimental growth rates to avoid unrealistic predictions. Although the model tends to overestimate IgG production, it remains accurate within the experimental order of magnitude, particularly with the HP dataset. This underscores the importance of including specific data to refine models.

The study highlighted specific metabolic demands, with lysine and other essential amino acids identified as limiting factors, and demonstrated that IgG production can be optimized by changing macroscopic parameters from the bioreactor such as the dilution rate and the biomass. Finally, by optimizing growth medium composition and nutrient costs, the model provides insights into enhancing IgG yields and production efficiency, demonstrating the crucial role of integrating experimental data with computational models in bioprocess optimization.

Funding

This work was supported by the Spanish Ministry of Science and Innovation through the projects DAMOCLES-PID2020-113969RB-I00/AEI/10.13039/501100011033 and TED2021-130449B-I00 , by the Aragonese Government under Programa de Proyectos Estratégicos de Grupos de Investigación (DisCo research group, ref. T21-23R ), and by a PhD fellowship from Diputación General de Aragón (DGA) to J.L.

CRediT authorship contribution statement

Jorge Lázaro: Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Teresa Joven: Writing – original draft, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Diana Széliová: Writing – review & editing, Validation, Supervision. Jürgen Zanghellini: Writing – review & editing, Validation, Supervision, Resources, Conceptualization. Jorge Júlvez: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Resources, Project administration, Investigation, Funding acquisition, Conceptualization.

Declaration of Competing Interest

All authors declare that they have no conflicts of interest relevant to the content of this manuscript. No financial, professional, or personal relationships influenced the research presented in this paper.