Metabolic Regulatory Network Kinetic Modeling with Multiple Isotopic Tracers for iPSCs

The rapidly expanding market for regenerative medicines and cell therapies highlights the need to advance the understanding of cellular metabolisms and improve the prediction of cultivation production process for human induced pluripotent stem cells (iPSCs). In this paper, a metabolic kinetic model was developed to characterize underlying mechanisms of iPSC culture process, which can predict cell response to environmental perturbation and support process control. This model focuses on the central carbon metabolic network, including glycolysis, pentose phosphate pathway (PPP), tricarboxylic acid (TCA) cycle, and amino acid metabolism, which plays a crucial role to support iPSC proliferation. Heterogeneous measures of extracellular metabolites and multiple isotopic tracers collected under multiple conditions were used to learn metabolic regulatory mechanisms. Systematic cross-validation confirmed the model's performance in terms of providing reliable predictions on cellular metabolism and culture process dynamics under various culture conditions. Thus, the developed mechanistic kinetic model can support process control strategies to strategically select optimal cell culture conditions at different times, ensure cell product functionality, and facilitate large-scale manufacturing of regenerative medicines and cell therapies.


INTRODUCTION
Cell therapeutics and regenerative medicines have the potential to treat and prevent diseases, such as cancers, and cardiovascular and hematologic diseases 1,2 .The cell therapy market is experiencing unparalleled growth and the projected worth of the relevant market is expected to be over $8 billion in 2025 3 .To meet the rising demand, cost reductions in manufacturing and increases in quality for human induced pluripotent stem cells (iPSCs) on a large scale are crucial for the success of cell therapies 4 .Unlike traditional biopharmaceuticals, the cells are the product with functional identity depending on the regulatory metabolic dynamic behaviors.
The efficacy of the cells is very sensitive to culture conditions.Variability in the culture conditions can lead to reduced yields and heterogeneous cell populations.Heterogeneous cell populations have been observed to increase the potential for tumor or teratoma formulation in the patient 5 .Traditional cell culture process control strategies often are ad hoc and rely solely on experimental approaches or PID controllers ignoring long-term effects.The optimal design and control of mammalian cell culture processes, especially iPSCs, can be laborious due to the limited understanding of cellular metabolisms and end-to-end cultivation process dynamics 6,7 .arXiv:2305.00165v2 [q-bio.MN] 26 Oct 2023  In order to determine the best culture conditions, accounting for cell life cycle, relationships between culture conditions and the output trajectory (i.e., cell functional identity and yield) are needed.Mechanistic dynamics model has the capability to predict cell outcomes (e.g., metabolic flux rates) and can be used to assess product quality and strategically select the optimal culture conditions at different times.
The proposed mechanistic model describes an iPSC metabolic network regulatory mechanisms.It is associated with the existing studies, including (1) metabolic flux analysis (MFA) 8,9 , (2) 13 C-MFA 10,11,12,13,14,15,16,17 , and (3) kinetic models 6,9,18,19,20 .In MFA, metabolic flux rates (such as substrate uptake rate and metabolite secretion rate) are estimated based on experimental measurements subject to stoichiometric constraints.Under the standard assumption of (pseudo) steady state for intracellular metabolites, the sum of all fluxes producing a metabolite is equal to the sum of all fluxes consuming that metabolite.This technique is frequently used to compare the metabolism of different cell lines to assess the activity of individual pathways under different cultivation conditions 8 .However, intracellular metabolites can change during cell proliferation 21 .
Stable isotope studies, integrating with MFA, have provided more detailed information on intracellular state and metabolic pathways, namely 13 C-MFA can precisely estimate metabolic reaction rates.Typically these flux rates are estimated from measured intracellular mass isotopomer distribution (MID) patterns and external metabolites concentrations.Common 13 C and 2 H isotopic tracers were utilized for investigating mammalian cell metabolism with a specific emphasis on the interpretation of isotopic labeling patterns 16 .Many existing studies are built on the assumption of isotopic and metabolic steady state 10,11,12,13,14,15,16 .Several modeling tools have been developed for situations involving isotopic and/or metabolic non-steady states 17,10,11 .
Another approach to model cell dynamics includes kineticmetabolic models 6 .The vast majority of the kinetic models use Monod and Michaelis-Menten expression formalisms modeling the metabolic regulation mechanisms 18,19,20,6,9 .For Chinese hamster ovary (CHO) cell cultures, kinetic-metabolic models include glycolysis, pentose phosphate pathway (PPP), tricarboxylic acid (TCA) cycle, respiratory chain, as well as the regulatory functions from energy shuttles (ATP/ADP) and cofactors (NADH, NAD + , NADPH, NADP + ); this high level of details is required to predict hypoxic perturbation 19,20 .For example, the study 9 proposed a Monod model characterizing the cell flux rate response to environmental perturbation.That study used limiting substrate kinetics to calculate the extracellular metabolite consumption/production rates and further infer metabolic flux rates built on a static MFA assumption, i.e., assuming the intracellular metabolites are at the pseudo-steady state.Even though several in silico metabolomic platforms have been developed for CHO cells 18,19,20,6,9 , a metabolic regulatory network simulator for iPSC has not been developed.
As lactate accumulation emerges as an inevitable outcome of iPSC metabolism, it becomes essential to explore the consequences of accumulated extracellular lactate on intracellular iPSC metabolism.Equally important is to comprehend how glucose concentration acts as a regulatory factor, influencing lactate production.Understanding the implications of these concentration variations is not only critical for the survival and proliferation of iPSCs 22 but also essential for maintaining the quality of cell products.Hence, the main objective of this study is to develop and validate a mechanistic model of an iPSC metabolic regulatory network that can precisely predict cell responses to environmental changes, with a particular focus on variations in extracellular glucose and lactate concentrations.
The proposed model primarily focuses on central metabolism, given its pivotal role in stem cell proliferation, biosynthesis, and overall functionality.To estimate model parameters, this model utilized experimentally obtained data from well-designed monolayer K3 iPSC cultures, as described in the study 4 .During the model development, the level of detail required for key pathways was explored, namely the PPP and branched amino acids uptake rates into the TCA cycle.In this paper, the proposed kinetic model structure with critical regulatory mechanisms will be discussed and validated through comparisons with the experimental data, with the aim of ensuring the model's ability to accurately predict cell responses to environmental changes.
The proposed mechanistic model, comprising the metabolic reaction network structure and flux regulatory mechanism, can be easily customized to suit various mammalian cell culture systems, e.g., embryonic stem cells (ESCs), CHO, etc.This adaptability arises from the underlying principles and general characteristics of metabolic regulation that are common across various mammalian cells.Through the collection of extracellular metabolite concentrations and intracellular MID measurements, this mechanistic model can be adeptly employed to analyze and forecast the behavior of other mammalian cell cultures with consistent accuracy and reliability.

DATA DESCRIPTION
Briefly, the data used to estimate the model parameters were derived from an experimental study of K3 iPSCs, where monolayer cultures were conducted in 6-well plates and 60 mm cell culture dishes 4 .The data encompassed extracellular metabolite measurements and isotopic labeling details for four distinct culture conditions, utilizing three tracers: glucose, glutamine, and lactate .The four different initial culture media concentrations were: (1) high glucose and low lactate (HGLL); standard initial conditions, i.e., 18.3 mM glucose and 0 mM lactate; (2) high glucose and high lactate (HGHL), i.e., 18.3 mM glucose and 20 mM lactate; (3) low glucose and low lactate (LGLL), i.e., 5.6 mM glucose and 0 mM lactate; and (4) low glucose and high lactate (LGHL), i.e., 5.6 mM glucose and 20 mM lactate.Table 1 presents the media concentrations used for each condition, encompassing glucose, lactate, and supplemented sodium chloride concentrations.Furthermore, 2.75 mM glutamine was included in each condition.These concentrations were chosen to ensure comparable growth rates across all conditions while maintaining pluripotency 4 .
Cell density, glucose, lactate, pyruvate, and extracellular amino acid concentrations were obtained for each condition at 0-, 12-, 24-, 36-, and 48-h.Intracellular amino acid MIDs were obtained from parallel labeling experiments that used [1,2- 13 C 2 ] glucose, [U-13 C 5 ] L-glutamine, and [U- 13 C 3 ] sodium L-lactate (when lactate was added to the culture media).The time-course measurements, as illustrated in Figure 1 , included the transition to metabolic and isotopic steady states 4 .These measurements were used to develop the metabolic regulatory network kinetic model for iPSC cultures.

TABLE 1
The media formulations for K3 iPSC cultures include initial concentrations of glucose, sodium L-lactate, and supplemented sodium chloride (NaCl).NaCl was added to balance osmolarity in cultures lacking sodium lactate.Additionally, the media contains 2.75 mM L-glutamine.characterizing cell response to environmental variations, is fitted by using time-course extracellular concentration measurements and MIDs.Finally, in Section 3.4, the assessment of fit and validation of the dynamic model for the iPSC cultures is discussed, where the objective function and goodness-of-fit will be presented.

Metabolic Network Kinetic Modeling
In this section, the general regulatory metabolic network kinetic model and time-course isotopic labeling simulation are presented.The kinetic model captures the dynamic changes in cell density, and extracellular and intracellular metabolite concentrations, where the changing rates of metabolites depend on the concentration of substrates and inhibitors.To incorporate the mass isotopic data, a dynamic isotopic labeling system was constructed, and a time-course isotopic labeling simulation was developed.Overall, the metabolic reaction network kinetic model provides a comprehensive understanding of the underlying mechanisms of the iPSC cultures.
(1) Metabolic Reaction Network Kinetic Model.The data from the study for K3 iPSCs were collected during the exponential growth phase 4 .Through experimental design and statistical analysis, our previous study 4 discovered uniform specific cell growth rates () across the four batch culture conditions (HGLL, HGHL, LGLL, and LGHL), maintaining pluripotency.This suggests that glucose and lactate concentrations within these ranges did not notably influence K3 iPSC growth rates and did not have a detrimental impact on human K3 iPSC pluripotency.Thus the classical cell density formalism was used to relate growth and the cell density, denoted by () (cells/cm 2 ), at time  as, ()  = (). (1) The extracellular metabolite concentrations at time  are denoted by a vector with dimension , i.e.,   () = (  1 (),  2 (), … ,   () ) ⊤ .Similarly, the intracellular metabolite concentrations of interest are denoted by a vector with dimension  at time  as   () = (  1 (),  2 (), … ,   () ) ⊤ .Thus, at any time , the state of extracellular and intracellular metabolite concentrations is denoted by To model the dynamic evolution of cell response to environmental fluctuations during the iPSC cultures, at any time , the specific reaction flux rates, represented by a vector with dimension , depend on both extracellular and intracellular metabolite concentrations, i.e.,

𝑣 𝑣 𝑣[𝑢 𝑢 𝑢(𝑡)]
Let  denote a ( + ) ×  stoichiometry matrix characterizing the structure of metabolic reaction network.Therefore, the dynamic evolution of extracellular and intracellular metabolite concentrations is modeled through a mass balance on the system of equations,

𝑑𝑢 𝑢 𝑢(𝑡) 𝑑𝑡 = 𝑁𝑣 𝑣 𝑣[𝑢 𝑢 𝑢(𝑡)]𝑋(𝑡).
(2) To capture the cell response to environmental perturbation, a Michaelis-Menten (MM) formalism based regulation model was used to characterize the relationship of metabolic flux rates depending on the concentrations of associated substrates and inhibitors.This allows leveraging the existing biology knowledge and facilitating the learning of regulation mechanisms of iPSC metabolic reaction network from the experimental data.
Specifically, the -th flux rate at the time  is modeled as, for  = 1, 2, … , , where the set Ω   represents the collection of activators (such as nutrients and substrates) influencing the flux rate   and the set Ω   represents the collection of inhibitors dampening   .The parameters  ,⋅ ,  ,⋅ , and  ,⋅ represent the affinity constant, the inhibition constant, and the maximum specific flux rate respectively.For iPSC cultures, the experimental data referenced were used to identify the critical activators and inhibitors, as well as learning the MM regulation model coefficients.
(2) Time-Course Isotopic Labeling Simulation.Different 13 C-labeling patterns are generated by different flux distributions   [  ()].Thus, incorporating the dynamic system of isotopic labeling patterns (MID), which is a vector containing the fractional abundance of each mass state of metabolites, into the metabolic network kinetic model can aid in understanding intercellular metabolic network mechanisms.However, modeling each individual atom as one system state variable is computationally expensive.To address this issue, the elementary metabolite unit (EMU) framework was proposed, and it is based on a highly efficient decomposition method that can identify the minimum amount of information needed to simulate isotopic labeling 24 .Basically, the EMUs are created by using a decomposition algorithm and form the new basis for generating system equations that describe the relationship between fluxes and isotope measurements; see more detailed information in studies 25,26,10,27,24,11,23 .
Based on the study 24,28 , the reduced system can be obtained after decoupling based on EMUs with size  = 1, 2, … ,  and connectivity.The time-dependent network was first identified: the decoupled EMUs with size  network,   () = {  ,   ,   ()} with  = 1, 2, … , .The vector ) representing the dependence between each vertex.The corresponding weight matrix   () varies with time .The non-negative (, )th element  (,)  () indicates the flux rate of the reaction producing -th EMU by consuming -th EMU at time .
Thus, at any time , the dynamic isotopic labeling system can be defined as: where the rows of the state matrix   () correspond to the MIDs of EMUs within and

iPSC Metabolic Reaction Network
The developed iPSC metabolic network model leverages central carbon metabolism from several previously published metabolic frameworks 5,3,20,8,19,4 and further learns from data.The iPSC metabolic network included glycolysis, TCA cycle, anaplerosis, PPP, and amino acid metabolism.For simplification, the reactions of PPP were collapsed into two reactions: Oxidative phase/branch and Non-oxidative phase/branch (i.e., No. 9 & 10 reaction in Table A1 ).During the exponential phase, greater than 90% of pentose-phosphate carbons is returned to glycolysis for the mammalian cell 21,29 .Thus, the synthesis of nucleotides and nucleic acids is considered as an insignificant contributor to the model.The metabolic network is shown in Fig. 2 and the metabolic stoichiometry is listed in Table A1 .The descriptions of metabolites and the enzymes conform to the Enzyme Commission Number (EC-No.)and are provided in Tables A3 and A4 , respectively.While this model primarily focuses on central carbon metabolism, it's worth noting that the additional model's structure, including metabolic reaction pathways (e.g., one-carbon metabolism, fatty acid oxidation, and nucleotide biosynthesis) and the associated flux regulatory mechanisms, could be easily incorporated.

Biokinetic Model of Flux Regulatory Mechanism
Due to the lack of energetic state and redox level measurements, each metabolic flux rate is modeled as dependent on the substrates and inhibitors concentrations using Michaelis-Menten model formalism.Similar to the existing studies 5,3,20 , for simplification purposes and also due to a lack of available data in the literature, a single affinity constant value is used for each metabolite.Additionally, reaction reversibility was considered during the iPSC culture model development for some key reactions.The final model is fully described in Table A2 .Several regulatory mechanisms (i.e., No. R1 to R7) are highlighted in Fig. 2 , which were incorporated into the metabolic flux kinetics model to improve the model's prediction capability and better characterize the responses to environmental variations.These key final reactions are described below, where these new dependencies are highlighted in under brackets … … ⏟⏟⏟ with the corresponding regulatory mechanism No., while the original nomenclature is outside the brackets.
• Lactate accumulation has previously been reported to reduce glycolytic activity by inhibiting hexokinase (HK) and phosphofructokinase (PFK) activity in mammalian cells, where lactate acts as a signaling molecule to downregulate PFK activity 30,31,32 .After evaluation of the  A3 .The enzyme description, including the Enzyme Commission Numbers (EC-No.)for each reaction, is listed in Table A4 .
experimental data, the model for HK was updated to include this inhibitory effect of lactate on it: • Since lactate inhibits glutaminase activity -the enzyme responsible for converting glutamine (GLN) to glutamate (GLU) 33,34 -the forward ( ) flux rate for the reaction, i.e., Gln ↔ Glu + NH 4 , was updated, • Several regulatory functions, adapted from the studies 19,18,20 , were evaluated using the experimental data to characterize activations and inhibitions.The regulatory mechanisms involved in glycolysis are described as: a) hexokinase inhibition by its product G6P, see R3 in eq (9); b) activation of pyruvate kinase by F6P, see R4 in eq (10), as well as c) the inhibition of lactate dehydrogenase reverse () reaction, see R5 in eq (11), was changed to include the G6P inhibition: • Both lactate and pyruvate transport across the plasma membrane are facilitated by proton-linked monocarboxylate transporters (MCTs) 29,35 .More favorable lactate transport kinetics may decrease pyruvate consumption under high lactate culture conditions 4,36 .Thus, the model of PyrT was updated to: • Since the transportation of extracellular-glutamine via cell membrane can be inhibited by intracellularglutamine, the model for GlnT was updated to:

Model Fit and Goodness-of-fit Assessment
The for  = 0, 12, 24, 36, 48 hours (h).To ensure the model's accuracy and minimize the impact of measurement errors, values that fell below the detection limit of the Cedex Bioanalyzer are treated as missing data.The intracellular isotopic labeling measurements at time  ′ are denoted as MID   ′ for  ′ = 24 and 48 hours.The MID measurements obtained from [U- 13 C 3 ] lactate for the high glucose high lactate and low glucose high lactate cultures exhibited considerable measurement errors, as reported in the study 4 .Consequently, these measurements were intentionally excluded from the training dataset used for developing the proposed model.
For model fitting, the objective is to minimize the weighted sum of squared residuals (SSR) between the available experimental data, which includes both extracellular metabolite concentrations (   ) and intracellular isotopic labeling (MID   ′ ), and the corresponding model-predicted values (   and MID   ′ ), where the weights are the inverse of the variances of the experimental measurements on extracellular metabolites and MIDs, i.e., (   ,    ′ ).The fitting criteria in eq (14) allows the model to systematically and effectively fit the metabolic regulation network at different times.
The initial conditions for extracellular metabolites were obtained from culture data.As the experimental data only includes the relative abundance and the MIDs, estimates for the initial intracellular metabolite concentrations were sourced from the BRENDA Enzyme Database 37 and related references for mammalian cells cultured under similar conditions.Similarly, the initial kinetic parameter values were taken from the BRENDA Enzyme Database 37 and relevant references for similar metabolic networks and pathways.
The goodness-of-fit for the developed iPSC metabolic kinetic model was assessed based on predictions using a chisquare test.Basically, the null hypothesis is that the fitted model can faithfully represent the iPSC culture metabolic mechanism.This test assumes that the minimized varianceweighted SSR follows a chi-square distribution with  degrees of freedom, where the degree of freedom  equals the number of observations minus the number of fitted parameters.In this study, the chi-square statistical test has -value much greater than 0.05 indicating the model predictions are not significantly different than the measured values.

RESULTS AND DISCUSSION
The proposed mechanistic model was developed based on a small set of experimental data (see Section 2) and the model predictive performance was then validated in two manners.First, to mimic the dynamic data collection and assess the rolling forecasts required for process control, at any time , historical data were used to fit the model up to  and then the fitted model was used to predict the remainder of the culture behavior.In Section 4.1, the model's ability to capture the dynamic evolution of iPSC metabolic characteristics was evaluated.The model was trained on experimental data collected over different time intervals (0-h to 12-h, 0-h to 24-h, and 0-h to 36-h), and used to predict the iPSC culture for up to 48-h.
Second, to assess the model's ability to generalize and extrapolate to new culture conditions, three datasets, randomly selected from (1) high glucose and low lactate (HGLL); (2) high glucose and high lactate (HGHL); (3) low glucose and low lactate (LGLL); and (4) low glucose and high lactate (HGHL), were used to predict the fourth dataset with different initial conditions in Section 4.2.Basically, three datasets were used to train the model, and then the model was used to predict the fourth dataset behaviors.This training/test dataset selection and evaluation strategy is built on the philosophy of cross validation (CV) that is often used in the literature on model selection 38,39 .Each culture condition was examined in succession using the other three datasets.Therefore, by validating the model's prediction performance under different experimental conditions, its robustness and usefulness for predicting iPSC characteristics in a variety of settings is verified.

Prediction of iPSC Culture Process
During the model validation process, the proposed model simulator was trained using experimental data collected over various time intervals (0-h to 12-h, 0-h to 24-h, and 0-h to 36-h).
The model's ability to predict the dynamic trajectories of key metabolites, including glucose, lactate, glutamate, glutamine, and pyruvate, as well as cell density, was evaluated by comparing the predicted values to experimental observations up to 48 hours.Fig. 3 provides a representative result of the cell culture process prediction for the case starting with HGLL.The results for the remaining settings, including HGHL, LGLL, and LGHL, are presented in Appendix Fig. B1 to B3 .
The Chi-square test was employed to assess the goodnessof-fit of the proposed mechanistic model trained using the historical data collected in different time intervals.The statistical test for the model trained on 0-h to 12-h of experimental data yielded a value of 33.6 with 47 degrees of freedom (value = 0.93).Similarly, for the model trained on 0-h to 23-h and 0-h to 36-h of experimental data, the test statistics was 82.6 with 107 degrees of freedom (-value = 0.96) and 138.1 with 167 degrees of freedom (-value = 0.95), respectively.All of the p-values exceed 0.05, suggesting that the proposed metabolic kinetics model fits well with the experimental dataset.The results in Fig. 3 demonstrate that the fitted model can provide improved predictions for most of the metabolite trajectories as more data are collected over time.
Overall, the predicted profiles of the iPSC cultures can closely track the dynamic patterns of the measured profiles and capture the cell culture dynamics.The experimental iPSC data were obtained from plates and dishes with limited space.Even though the overall growth rates were the same across the four conditions by experimental design, the growth rate decreased with time as nutrients were consumed and metabolic wastes accumulated.It is important to note that the developed mechanistic metabolic kinetic model is a simplified representation of the iPSC carbon central network, there could likely be some key metabolic regulations and biochemical reactions that were not incorporated.The less accurate predictive performances for lactate (under low lactate conditions) and glutamate may be due to the omission of key reactions or regulatory functions.For example, as suggested by the calibrated model in previous studies 19,18,20 , the membrane transportation of glutamate for mammalian cells is significantly impacted by the energetic state (ATP and ADP), and the forward and reverse conversion of lactate to pyruvate is influenced by the redox level (NAD and NADH).
Since the redox level and energetic state-related measurements were not collected for the iPSC cultures, these factors were not accounted for in the current model.The study 9 found for CHO cells, a redox parameter assisted with lactate predictions; however, the redox parameter implemented by this study was unmeasurable and thus would be difficult to incorporate into a control model.Furthermore, the lower prediction performance for glutamate may be attributed to its involvement in multiple reactions (see Fig. 2 ), leading to error accumulation

Prediction across Different Initial Conditions
To assess the predictive capability of the proposed mechanistic model across various initial culture conditions and enhance the understanding of iPSCs response to environmental changes, specifically glucose and lactate, a cross-validation approach was employed.This approach first trained the model with three of the four datasets (i.e., HGLL, HGHL, LGLL, and HGHL).The remaining dataset behavior was then predicted.This approach provided valuable insight for optimizing Design of Experiments (DoE) strategies in the realm of in silico experimentation.
The model cross-validation prediction for the HGLL culture was evaluated using the model trained on HGHL, LGLL, and LGHL datasets, as illustrated in Fig. 4 .The other three acrossvalidation predictions are provided in Appendix Fig. C4 to C6 .Cell density and key metabolites (i.e., glucose, lactate, glutamate, glutamine, and pyruvate) were predicted and are shown with the average of the measured values with standard deviations.To evaluate the goodness-of-fit, the chi-square test SSR statistics is 783 with 802 degrees of freedom with -value = 0.67 much greater than 0.05.This p-value indicates that the metabolic kinetic model can faithfully represent the iPSC culture regulatory mechanisms under different levels of glucose and lactate concentrations.
The predictions of intracellular MID from the [1,2- 13 C 2 ] glucose and [U-13 C 5 ] glutamine tracer at 48-h under the control (HGLL) culture condition are shown in Fig. 5 and Fig. 6 , respectively.Even though there is some prediction error, the simulation model can correctly predict the dynamics and interdependencies of multivariate iPSC culture process metabolism.The developed metabolic kinetic model can faithfully predict the cell response to environmental perturbation, which could guide the strategic feeding strategy for the integrated iPSC culture process for future research.In addition, by incorporating even the simplified PPP reactions, the M+1 isotopic labeling pattern for pyruvate and related metabolites from [1,2- 13 C 2 ] glucose tracer can be well-predicted.
The prediction of metabolic flux maps for K3 iPSC under control culture at 24-h and 48-h are shown in The prediction of metabolic flux dynamics for K3 iPSC culture reveals several significant observations.First, the glycolytic efficiency, which measures how effectively K3 iPSCs  utilize glucose as an energy source, remained consistently high across all conditions, exceeding 1.6 and approaching the theoretical maximum of 2.0 moles of lactate produced per mole of glucose consumed.This suggests that K3 iPSCs are efficiently converting glucose into lactate with minimal byproduct production.Second, the [U- 13 C 3 ] lactate tracer showed that K3 iPSCs consumed and metabolized lactate in the high lactate cultures.However, in all four culture conditions, there was a net production of lactate, indicating that K3 iPSCs produced and also consumed lactate in the high lactate culture conditions.Third, in the control experimental setting, as the iPSC cultures approached the late exponential phase, there was a gradual increase in the flux rates of the TCA cycle, consistent with findings reported in a previous literature study 21 .Finally, the mechanistic model successfully predicted the reduced glutamine consumption fluxes under high lactate culture conditions, along with a decrease in the conversion of glutamine to AKG.Consequently, the model also predicted slightly lower fluxes through key TCA cycle reactions such as -ketoglutarate dehydrogenase (AKGDH), succinate dehydrogenase (SDH), malate dehydrogenase (MDH), fumarase (FUM), and citrate synthase (CS).
Table 2 lists the predictions of biomass-specific uptake and production rates for critical extracellular metabolites at 12-h and 36-h.The table also includes the average flux between 12h and 36-h, as determined by the Extracellular Time-Course Analysis (ETA) software 40 , for comparison.The general mass balance equation, constructed by ETA, which describes the dynamics in the concentration of the -th extracellular metabolite under batch growth conditions, is given by where   () represents the concentration of the -th metabolite at time ,   is the specific rate of metabolite consumption or production, and () stands for cell density at time .This equation assumes no degradation of metabolites.Notice that ETA results are based on the assumption of metabolic steadystate, meaning that the consumption or production rate   of any -th metabolite remains constant and it is not influenced by the changes in extracellular or intracellular concentrations.
This assumption explains why most flux estimations using the ETA approach for 12-h to 36-h fall within the range of the model predictions at 12-h and 36-h.Under conditions of limited substrate availability, there are significant changes in consumption rates, as exemplified by glucose in LGLL and LGHL cultures, as well as pyruvate in all four conditions.Unlike averaging approaches (e.g., ETA model), the proposed metabolic kinetic model effectively captures these essential dynamic behaviors.This capability holds promise for facilitating future research on the end-to-end cell culture process optimal control to improve high-quality iPSC production.Additionally, as stirred suspension cultures become a common choice for large-scale iPSC manufacturing 41,42 , this study can also contribute to the understanding of cell behavior within aggregates, as the well-established pattern involves low glucose and high lactate concentrations at the aggregate center 43 .
By modeling the regulatory mechanisms outlined in Section 3.3, the results presented in Table 2 are consistent with observations reported in existing studies 19,18,20,30,31,32,33,34,29,35 .The ability of the proposed model to accurately predict cell dynamic behavior and response under varying environmental conditions reaffirms its reliability and underscores its potential to enhance the optimization and control of iPSC culture processes.Specifically, the model indicates that the consumption rate of glucose increases with higher extracellular glucose concentrations (e.g., HGLL vs. LGLL and HGHL vs. LGHL), while it decreases with higher extracellular lactate concentrations (e.g., HGLL vs. HGHL and LGLL vs. LGHL).These trends can be attributed to regulatory mechanism R1 as described in equation (7).Regarding lactate, the model reveals that the net production rate is lower when a certain amount of lactate exists in the cultural environment (e.g., HGLL vs. HGHL and LGLL vs. LGHL).In the case of pyruvate, under high lactate culture conditions, more favorable lactate transport kinetics of MCTs result in decreased pyruvate consumption (e.g., HGLL vs. HGHL and LGLL vs. LGHL).These observations are associated with regulatory mechanism R6 as explained in equation (12).Furthermore, the presence of lactate dampens the conversion of glutamine to glutamate.Consequently, under high lactate culture conditions, the consumption rate of glutamine and the production rate of glutamate decrease (e.g., HGLL vs. HGHL and LGLL vs. LGHL).These behaviors are linked to regulatory mechanisms R2 in equation ( 8) and R7 in equation ( 13).

CONCLUSIONS
In this paper, a metabolic kinetic model is described that characterizes the time-varying dynamics and regulatory mechanisms of the iPSC cultures.This model detailed central carbon metabolism, including glycolysis, TCA cycle, PPP, anaplerosis, and key amino acid metabolism.The iPSC metabolic regulatory network was calibrated using extracellular metabolite concentrations and intracellular isotopic data for multiple tracers (i.e., [1,2- 13 C 2 ] glucose and [U- 13 C 5 ] glutamine tracers).These time-course measurements were collected at multiple time points for the different culture conditions.The validation results demonstrate that the developed metabolic kinetic model can provide a reliable prediction on iPSC metabolic response to changes in extracellular glucose and lactate levels.
This model holds substantial potential for advancing the understanding of complex intracellular regulatory mechanisms within iPSC cultures.Furthermore, it empowers continuous monitoring and precise control of cell cultures, particularly advantageous in the realm of large-scale manufacturing.In summary, the proposed mechanistic model holds the potential to facilitate the optimization and control of iPSC cultures, potentially even at large scales, ensuring the survival, productivity, and quality of iPSC-derived cell products.LGLL (nmol/10

A APPENDIX: TABLE
The following tables summarize detailed information about the constructed iPSC metabolic network and developed kinetic model.
• The iPSC metabolic network contains the major reactions for glycolysis, the TCA cycle, anaplerosis, PPP, and amino acid metabolism; see the specific reactions in Table A1 .The reactions of PPP are collapsed into two: Oxidative phase/branch and Non-oxidative phase/branch (i.e., No. 9 & 10 reactions).
• The iPSC metabolic flux rate regulation biokinetic model is summarized in Table A2 .The Michaelis-Menten model is used to characterize how key activators and inhibitors influence regulatory mechanisms in each reaction.
• The descriptions of the metabolites, considered in the developed metabolic kinetic model, are listed in Table A3 .
• For the developed metabolic kinetic model, the descriptions of the enzymes with Enzyme Commission Number (EC-No.)are provided in Table A4 .They are associated with metabolic flux reaction rates.

FIGURE 2
FIGURE 2 An illustration of the iPSC regulatory metabolic network (Created with BioRender.com).Glycolysis, PPP, TCA and anaplerosis and amino acid utilization are shown.Additionally, the regulatory mechanisms included in the dynamic model are shown.The metabolites descriptions are listed in TableA3.The enzyme description, including the Enzyme Commission Numbers (EC-No.)for each reaction, is listed in TableA4.

FIGURE 3
FIGURE 3 Dynamic model trained with historical data collected under different time intervals -Cell characteristic predictions for high glucose and low lactate cultures.(A) Cell density, (B) Glucose, (C) Lactate, (D) Glutamine, (E) Glutamate, and (F) Pyruvate.Times 0-h to 12-h (green dotted line); Times 0-h to 24-h (brown dashed line); and Times 0-h to 36-h (purple solid line).The detection limit of the Cedex Bioanalyzer is shown as the red dash-dot line.

Fig 7 .
The prediction performance of metabolic concentration trajectory and flux maps for K3 iPSC under HGHL, LGLL, and LGHL are provided in Appendix Fig.D7to Fig.D9.Since no restrictive steady-state assumption was required for the developed metabolic kinetic model, the flow-in flux is not required to be equal to the flow-out flux for each metabolite.

FIGURE 4
FIGURE 4 Dynamic model trained with three other datasets with varied initial glucose and lactate concentrations for cell characteristic predictions in high glucose and low lactate cultures.(A) Cell density, (B) Glucose, (C) Lactate, (D) Glutamine, (E) Glutamate, and (F) Pyruvate.The detection limit of the Cedex Bioanalyzer is shown as the red dash-dot line.

FIGURE 5
FIGURE 5 MID measurements and predictions for the [1,2-13 C 2 ] glucose tracer at 48-h for high glucose low lactate cultures using the dynamic model trained on the other three data sets compared to literature measurements.MIDs shown have been corrected for natural abundance.

FIGURE 6
FIGURE 6 MID measurements and predictions for the [U-13 C 5 ] glutamine tracer at 48-h for high glucose low lactate cultures using the dynamic model trained on the other three data sets compared to literature measurements.MIDs shown have been corrected for natural abundance.

FIGURE 7
FIGURE 7 Metabolic flux maps for K3 iPSC for the high glucose and low lactate cultures at 24-h and 48-h.Predicted fluxes are given in nmol/10 6 cells⋅h.The line thicknesses represent the relative fluxes.

Figures
Figures B1 to B3 depict the model's predictions for iPSC cultures in three scenarios: high glucose and high lactate (HGHL), low glucose and low lactate (LGLL), and low glucose and high lactate cultures (LGHL).Predictions are based on training the model using data collected over various time intervals.

Figures
Figures C4 to C6 depict prediction results across various initial conditions, including high glucose and high lactate (HGHL), low glucose and low lactate (LGLL), and low glucose and high lactate (LGHL).The model's fitting is based on data from the remaining three cases.

FIGURE D8
FIGURE D8 Metabolic flux maps for K3 iPSC for the low glucose and low lactate cultures at 24-h and 48-h.Predicted fluxes are given in nmol/10 6 cells⋅h.The line thicknesses represent the relative fluxes.

Component HGLL HGHL LGLL LGHL Glucose
23hematic for cell metabolic and isotopic states when time-course measurements of extracellular concentrations and intracellular MIDs were taken during the experimental iPSC cultures (Adapted from Jazmin and Young23and created with BioRender.com).Metabolic steady state is reached prior to isotopic steady state 4 .
The concentration matrix   () is a diagonal matrix whose elements are pool sizes corresponding to EMUs represented in  ()  .The construction of   () and   () are based on the decoupled EMU reaction network   ().The system matrix   () with size | ()  | × | ()  | and matrix   () with size | ()  | × | ()  | describe the metabolic network with elements defined as follows: () at time .The input matrix   () is analogous but with rows of input/carbon sources EMUs within  ()  at time .
proposed iPSC metabolic network kinetic model focuses on characterizing central carbon metabolism, including significant reactions from glycolysis, the TCA cycle, anaplerosis, PPP, and amino acid metabolism, along with regulatory mechanisms, across distinct cultural environments characterized by varying levels of glucose and lactate concentrations.The final model developed is shown in TableA1.This model includes 30 metabolic reactions with 32 variables (metabolites' concentration) in the reaction equations.The iPSC metabolic flux rate regulation biokinetic model, presented in TableA2, incorporates key activators and inhibitors for each reaction, describing the regulatory effects through the Michaelis-Menten model.The proposed iPSC culture kinetic model estimates model parameters using both extracellular and MID data collected over time under the four different culture conditions.The extracellular metabolite concentration measurements at time  are denoted as

TABLE A1
Reactions of the metabolic network

TABLE A3
Description of the Metabolite

TABLE A4
Description of the Enzyme B APPENDIX: CELL CHARACTERISTIC PREDICTIONS