Using Kinetic Modelling to Infer Adaptations in Saccharomyces cerevisiae Carbohydrate Storage Metabolism to Dynamic Substrate Conditions

Microbial metabolism is strongly dependent on the environmental conditions. While these can be well controlled under laboratory conditions, large-scale bioreactors are characterized by inhomogeneities and consequently dynamic conditions for the organisms. How Saccharomyces cerevisiae response to frequent perturbations in industrial bioreactors is still not understood mechanistically. To study the adjustments to prolonged dynamic conditions, we used published repeated substrate perturbation regime experimental data, extended it with proteomic measurements and used both for modelling approaches. Multiple types of data were combined; including quantitative metabolome, 13C enrichment and flux quantification data. Kinetic metabolic modelling was applied to study the relevant intracellular metabolic response dynamics. An existing model of yeast central carbon metabolism was extended, and different subsets of enzymatic kinetic constants were estimated. A novel parameter estimation pipeline based on combinatorial enzyme selection supplemented by regularization was developed to identify and predict the minimum enzyme and parameter adjustments from steady-state to dynamic substrate conditions. This approach predicted proteomic changes in hexose transport and phosphorylation reactions, which were additionally confirmed by proteome measurements. Nevertheless, the modelling also hints at a yet unknown kinetic or regulation phenomenon. Some intracellular fluxes could not be reproduced by mechanistic rate laws, including hexose transport and intracellular trehalase activity during substrate perturbation cycles.


Introduction
Saccharomyces cerevisiae, also commonly known as baker's yeast, has been used by mankind for thousands of years for the production of relevant beverages, foods and chemicals. However, despite its extensive use in industry [1][2][3], scaling new S. cerevisiae production processes to industrial scale poses several interesting and fundamental challenges. The source of most challenges is spatial inhomogeneities due to mixing limitations in large-scale bioreactors leading to gradients throughout the reactor. A cell dispersed in the reactor is therefore exposed to rapid changes in its extracellular environment, which in turn will impact intracellular metabolic regulation [4,5]. Similarly, the natural habitat will commonly have oscillations with respective to perturbations of environmental conditions such as temperature, pH and substrate availability.
Although dynamic conditions are encountered for industrial applications as well as environmental habitats, many physiological studies on yeast are performed under (pseudo-) steady-state conditions. Clearly, with the vast available reference data, reliable measurements and reproducibility, steady-state (SS) experiments are very useful in the quantification of intracellular fluxes. However, for the identification of in vivo kinetic parameters, dynamic metabolic experiments are required [6]. To bridge this gap, dynamic perturbation experiments can be performed, and many studies have focused on elucidating the metabolic response from single pulse (SP) experiments [6][7][8][9][10][11][12].
While this SP approach is very useful for the identification of kinetic parameters of networks adapted to the pre-perturbation limited steady-state, it cannot describe adaptations that may occur upon the subsequent perturbations observed under industrial conditions [5]. To emulate such an environment, a system of periodic substrate perturbations was developed. This regime produces repetitive substrate concentration gradients in time, which allow for accurate and reproducible sampling of the intracellular metabolism ( Figure 1) [13,14]. Suarez-Mendez et al. (2014) used this repeated substrate perturbation setup to monitor the in vivo metabolic activity during cycles of 400 s. At this timescale, it is assumed that the metabolic response within one cycle is mainly governed by metabolic interactions, as enzyme concentrations will remain basically constant during these 400 s [15]. In a cycle, feed was provided block-wise, i.e., 20 s feeding followed by 380 s of no feed ( Figure 1).  After a chemostat phase (reference steady-state) of 50 h, a block-wise feed is applied in a 400 s cycle at the same average substrate supply and dilution rate for another 50 h (adapted from [13]). On the left, a schematic overview of the feed rate during chemostat and repeated substrate perturbation regimes is shown. On the right, the resulting extracellular substrate concentration profile in the fermentation broth is shown.
Under such dynamic substrate conditions, S. cerevisiae cultures show different metabolic phenotypes compared to SP or SS cultures [13]. After a pulse in the repeated substrate perturbation regime, an increase rather than a decrease in ATP, no ethanol production and no accumulation of glycolytic metabolites was observed. These differences in metabolic response suggest a proteomic adaptation induced by the prior dynamic growth conditions [16][17][18]. Especially, translational regulation can lead to condition-specific proteome compositions [19]. In fact, distinct proteome compositions have been observed under changing glucose availability conditions [20], both for sugar transporters and intracellular enzymes [21,22], and distinct isoenzymes have different kinetic properties that can include glucose sensitivity as well [23]. However, the mechanisms behind this adaptation are not well understood.
For the identification of kinetic parameters and putative regulation mechanisms, quantitative data from different approaches will be required. Especially, to generate comprehensive models, in addition to intracellular concentrations, carbon tracing is required to identify bidirectional, cyclic or parallel reactions [24,25]. Work with 13C-labelling indicates that storage metabolism is a major metabolic sink upon changes in the glycolytic flux, with on average 15% of the carbon flux being diverted through the glycogen and trehalose cycles [26]. This diversion of flux is in accordance with earlier studies into the importance of the trehalose cycle under SP conditions [9].
Understanding metabolism, especially under dynamic conditions, requires integration of stoichiometry and enzyme kinetics. Such kinetic metabolic models have been developed for S. cerevisiae using either in vitro and in vivo parameters [22,27], additional allosteric regulation [28] or subpopulation dynamics [9]. Nonetheless, the conditioning of the cells was mostly at steady-state [29], leading to putative mismatches when applied to large-scale, dynamic cultivation conditions. Under dynamic conditions, further pathways have been described to play a regulatory role-glycogen and trehalose metabolism. For both, the reaction stoichiometry is known, and in vitro parameters have been broadly studied [30][31][32][33], but no in vivo based parameter values have been derived.
Especially for cyclic pathways, such as the trehalose cycle, quantifying in vivo parameters can be challenging as both in-and outfluxes influence the concentration change and no in-or outflux is directly observable. This correlation, plus the fact that the networks are getting larger, leads to a danger posed by local minima and ill-conditioning and consequently sloppy parameter estimates [34]. To overcome this challenge and identify a minimal set of necessary changes in kinetic constants, the divide-and-conquer approach has been developed. Here, a decomposition of the global estimation problem into independent subproblems [35] is used. Furthermore, to consider the already known parameter values for the enzymes under study [22,29], L1 or Tikhonov regularization can favour a given parameter set as long as experimental data are properly reproduced [36][37][38].
Here, we specifically studied the impact of proteome adaptation to substrate perturbations on the changed metabolic response. To this end, we expanded upon existing state-of-the-art kinetic models, combining both metabolome and fluxome, as well as 13 C enrichment data, to evaluate which proteome changes are most relevant to explain the experimentally observed change in metabolic response.

Strain and Growth Conditions
The haploid yeast Saccharomyces cerevisiae CEN PK 113-7D strain was grown at 30 o C, pH5, first in a batch phase and then in chemostat at dilution rate 0.1 h −1 [13]. The repetitive substrate perturbation regime began after five residence times and consisted of 20/380 scycles in which a feed was added in the first 20 s. The concentration of this feed was 20 times higher than the one of the chemostat phase to ensure that the culture would overall receive the same amount of glucose. Data were collected after 20 cycles. For further reference, see [13].

Experimental Datasets Used in This Work
The experimental datasets used in this work consisted of metabolite concentration measurements [13] and the respective calculated reaction rates [26]. Samples were collected more frequently during the first part of the cycle, since substrate concentration, and thus network dynamics changed more rapidly during this part [13]. Suarez-Mendez and colleagues [13] took samples for the extracellular measurements using a cooled syringe to enable quenching and fast filtration (see details in Mashego et al. [39]). The filtered sample was then analysed using HPLC and enzymatics analysis [13]. For the measurement of intracellular metabolite concentrations, samples were rapidly withdrawn using a selfdeveloped sampling device. The broth was first quenched in cold methanol; cells were washed and metabolites were extracted by boiling in ethanol (see details: Wu et al. [40] and Douma et al. [41]). The extract was then analysed using GC-MS and LC-MS [42].
Extracellularly, concentrations were measured for carbohydrates glucose and trehalose. Intracellularly, concentrations were measured for carbohydrates involved in glycolysis, the trehalose cycle, PPP, the glycerol branch and the TCA cycle and for adenosine nucleotides. Dynamic fluxes were estimated as piece-wise linear functions [43] using a consensus stoichiometric model for yeast [44]. Fluxes were estimated for glycolysis, carbohydrate storage metabolism (trehalose and glycogen cycles), PPP and the TCA cycle [26].

Model Description
A kinetic model of yeast central carbon metabolism was adapted in this work [29]. The original model contained the reactions that compose glycolysis, glycerol branch, a simplified trehalose cycle. Reactions of the PPP, the TCA cycle and uptake of glycolytic metabolites for biomass production were lumped as sink reactions (similar to [45]). The following modifications were made to represent the complexity of carbon storage metabolism seen in the data and adapt the sink reactions of the TCA to the repeated substrate perturbation setup: 1.
New reactions were added to represent a complete trehalose cycle and glycogen synthesis and degradation: • The A-glucoside transporter (AGT1) mobilizes trehalose between the extracellular space and cytosol [32]. Its reaction rate was modelled using reversible uni-uni MM kinetics. Since the experimental data pointed at a decay in its activity during the cycle but it did not contain any information on possible inhibitors, an inhibitory effect of T6P was added as a proxy of an increasing flux through the trehalose cycle. • A vacuolar transport of trehalose was added to mobilize trehalose between cytosol and vacuole-like compartments. Even though trehalose can be compartmentalized in vesicles in the cytosol, the kinetics of the process are not known.
Here it was assumed that reversible MM kinetics determine this process, as with AGT1. • Acid trehalase (ATH1, EC 3.2.1.28) degrades trehalose to glucose. It acts in more acid environments that the cytosol, such as the vacuole or the intracellular space [32], even though its location is still under debate. This reaction was modelled using irreversible MM kinetics. Similar to AGT1, inhibition by T6P was added. • UDP-Glucose phosphorylase (UDPG, EC 2.7.7.9) carries out the reaction from G1P to UDP-glucose, which is later used as substrate for glycogen synthesis. This reaction was adapted from [46] and modelled using an ordered bi-bi mechanism. • Glycogen synthesis was not modelled by enzymatic kinetics but interpolated from the experimental data in this study, with an added UDP-glucose saturation factor. • Glycogen degradation was also interpolated from the experimental data in this study, with an added UDP-glucose saturation factor.

2.
The sink reactions were optimized for chemostat growth [47] in the previous model. At a dilution rate of 0.1 h −1 , the fluxes observed were higher than the ones seen under the repeated substrate perturbation regime. As a result, the flux simulated in repeated substrate perturbation towards the TCA cycle via the sink of pyruvate was overestimated, resulting in a lesser flux towards the fermentative direction and more CO 2 being produced than measured. A factor was added to the reaction accounting for the pyruvate sink to reduce its flux and fit the CO 2 produced in the experiment.
In the following section, further details on model implementations are discussed.

System of Ordinary Differential Equations
The model consists of a series of ordinary differential equation representing the mass balances for each metabolite in the model. The model contained three compartments: cytosol, vacuole and extracellular space. The metabolites that are part of glycolysis, trehalose and glycogen cycles, glycerol branch and cofactors metabolism were located in the cytosol. Trehalose could be compartmentalized in the vacuole or secreted to the extracellular space, and glucose was modelled extracellularly as well. Meanwhile, the amount of carbon structure inside the cytosol depended on the inflow of glucose and outflow of the system. Moiety conservations were used for cofactors as in [22]. The sum of adenosine and nicotinamide adenine nucleotides (ATP + ADP + AMP and NAD + NADH, respectively) was kept constant in the cell. The model mass balances can be seen in detail in Supplementary File S1.

Reaction Rate Equations
The reaction rate equations used in this model followed Michaelis-Menten kinetics in most cases, but there were exceptions: PFK kinetics are affected by multiple regulators and the alternation between tense and relaxed state [27], PYK and PDC follow Hill-type kinetics [28] and glucose transport occurs by facilitated diffusion (equilibrium constant equals 1). Additionally, multiple allosteric regulations occur in the network, both activation and (competitive) inhibition. Reactions were made reversible, except hydrolysis reactions, due to their remarkably negative Gibbs energy (reference) and the sink reactions in the model. Reaction rates were expressed in (mM s -1 ). This unit refers to the intracellular volume. It could be correlated with the biomass dry weight via the biomass volume fraction (defined as 0.002 L per gram biomass dry weight). The kinetic rate expressions can be seen in detail in Supplementary File S1.

Simulation Setup
The simulations were performed in three steps aimed at resembling the experimental process that cells underwent in the experiments in [13]. The first step consisted of a chemostat, which also served to confirm that the system remained in a physiological realistic steady state. The second step consisted of the substrate perturbation cycles. A total of 20 repetitive cycles were simulated in which glucose was fed for the first 20 s of the cycle without any outgoing flux. For the rest of the cycle, no glucose was fed, and the outgoing flux lasted until the same amount of volume increased in the first 20 s had been emptied, by approximately 260 s of the cycle. Afterwards, both incoming and outgoing fluxes were kept at zero. After running for 20 cycles, the resulting simulation was compared to the experimental metabolite concentrations in reaction rates obtained in [13]. The third step concerned the simulation of the enrichment profiles reported in [26]. This extra simulation was run after the second step, and 99% of the inflow of glucose was 13 C labelled. All the simulations were carried out with the abovementioned mechanistic model. Matlab version 9.3.0.713579, R2017b and the ode15s solver were used. A graphical view of how simulations were constructed can be found in Supplementary File S1.

Implementation of 13C-Labeling Data Simulations
As described by [25], reactions of a metabolic network can be correctly represented with mass isotopomeric models if there are no cleavage reactions present because in the latter the position of the labelled carbon(s) is decisive to define the isotopomers of the following metabolites. In the interest of providing an accurate model without overcomplicating it, the kinetic model is thus expanded to labelled carbon enrichment instead of simulations of isotopic transients, as the simulation of the isotopic transients requires determination of the distribution of the possible C-labelled atoms for every metabolite of the network and additionally accounting for the bidirectionality of reactions in the isotopomer balance equations, which is no longer as trivial as building an enrichment model [48].
For each carbon-based metabolite in the model, a mass balance was added to account for its respective 13 C-labelled fraction. In this secondary mass balance, the input and output reaction rates were the same as for the total metabolite concentration mass balance, but it was multiplied by the respective fraction of labelled metabolite. Whereas the metabolite total concentrations depend solely on the enzymatic rates, the metabolite labelled concentrations also depend on the fractions of other labelled metabolites (labelled concentration/total concentration) that the reactions use as substrate. Moreover, now the enzymatic fluxes need to be adjusted for reversible enzymatic reactions.
For mass balances equations of total metabolite concentrations, reversible enzymatic fluxes are defined with a positive value in one direction and negative in the reverse direction. However, to implement the mass balances equations of labelled metabolite concentrations, these are multiplied by metabolite's labelled fractions, so they need to be always positive. Thus, enzymatic fluxes that change directions during the simulation are implemented as a forward and backwards flux. These are only used in the equations for labelled metabolite concentrations and are defined in the model with conditional statements. An example of the mass balance of labelled and unlabelled acetate can be seen below. Reaction reversibility is already accounted for inside the calculation of the reaction rate: where L refers to the labelled metabolite fraction. Enrichment simulations were performed after the 20 repetitive cycles were simulated. The experimental data consisted of percentage level of 13 C enrichment over the cycle time, obtained from [26].

Parameter Values Used in This Work
The initial parameter values were obtained from the original model [29]. These parameter values had been fitted to experimental data (metabolomic and fluxomic) from data at different steady states [47] and SRE [9]. A subset of the parameter values including enzymes HXT and GLK, the ones involved in the trehalose cycle and ATPase kinetics were estimated in this work. The values of the kinetic rate expressions used in this work can be seen in Supplementary File S1.

Estimation of In Vivo Parameters
Some reactions in the model underwent changes during the substrate perturbation cycles. For instance, the proportions of the isoenzymes HXK/GLK changed. To account for its effect on kinetic parameters, such as K M or K cat , reaction parameters were estimated for HXT and HXK/GLK. Additionally, trehalose cycles parameters were estimated, since the final structure of the cycle was different than the initial and the ATPase reaction rate constant needed to change as well. For initial parameter guesses, the initial parameter values from [29] were used. The nonlinear least-squares solver lsqnonlin from the Optimization Toolbox, using an interior reflective Newton method [49], was used to estimate the parameters by minimizing the error between measured and simulated data during the transient experiment.

Design of the Cost Functions: Combination of Enzymes and Weighting Factors
Experimental quantification of isoenzymes pointed at the couple HXK/GLK experiencing the biggest deviation prior to and after the substrate perturbation cycles, but minor changes were also confirmed for the other enzymes in glycolysis. Still, estimation of all the kinetic constants in glycolysis simultaneously was undesirable due to the appearance of parameter dependencies [50] which could lead to unphysiological parameter values. Therefore, to describe the changes in the experimental data with only the essential number of enzymes changing, the parameters were fitted to the data on multiple occasions. In each of those, a different selection of enzymes and cost function weighting factors was used: • Selection of enzymes: Multiple enzyme combinations were tested. These combinations contained the trehalose cycle and added different enzymes from glycolysis each time.
The combination selected was the one that described experimental data properly while making physiological sense (such as including the changes in HXK/GLK) and having the smallest number of enzymes possible. Simultaneously, random combinations of enzymes were also tested to confirm results and give robustness to the method.
• Combination of weighting factors: It was not clear at first if it would be possible to describe all the experimental data simultaneously. For this purpose, each of the abovementioned enzyme combinations was run repeated times, each of them with a different set of weighting factors. The errors for every metabolite were normalized so that they would contribute with the same weight to the cost function. Additional weighting factors changed these weights in three orders of magnitude at most.

Design of the Cost Functions: Regularization
Parameter dependencies could still appear for a selection of enzymes or within a single enzyme, even though less generalized than if all enzymes had been optimized together. To avoid this problem, parameters were estimated again after the previous round of data fit which was only based on selecting enzymes and cost function weights. This time, L1-type regularization was implemented to force the parameter estimates closer to the initial parameter set, as long as experimental data could be properly fit, helping to identify the important parameter changes for a specific dataset [36,37]. In this way, only the necessary parameters needed to change to fit the repeated substrate perturbation data. The regularization factor λ was applied to the cost function in the following manner: error estimation = error data + error parameters (3) error data = data simulated − data experimental (4) error parameters = −λ parameters reference − parameters estimated (5)

Cells Grown under Continuous and Dynamic Substrate Conditions Demonstrate Different Enzymatic Levels and Metabolic Responses-Experimental Observations
As mentioned before, cells exposed to a block-wise feeding (repeated substrate perturbation regime) showed a remarkably different response compared to glucose-limited cells from chemostat conditions. The proteome during both conditions was measured and subsequently analysed on their composition [51]. Major changes were observed within the glycolytic and transporter enzymes, specifically in the expression levels of hexose transporters (Hxt2p, Hxt3p, Hxt5p, Hxt13p), hexokinase (Hxk1p, Hxk2p) and glyceraldehyde dehydrogenase (Tdh1p, Tdh2p, Tdh3p) ( Figure 2). The expression of hexose transporters has been shown to correlate with the (maximal) substrate uptake rate [52]. The observed decrease in protein concentration can be interpreted as an adaptation to limit rapid influx of glucose upon glucose pulse. In contrast, glucokinase (GLK) is slightly upregulated. HXK is highly regulated through inhibition by T6P; however, GLK is not inhibited by trehalose-6-P (T6P) up to a level of 5 mM [53]. As such, the lower concentration of HXK in combination with the upregulation of GLK will likely result in an adaptation in the regulation of the glycolytic flux. The downregulation of upper glycolysis (HXK), in combination with the upregulation of lower glycolysis (TDH), may additionally allow for improvement of flux capacity through glycolysis upon glucose influx [9]. Ref. [54] observed that long-term adaptation of the proteome composition had a major influence on the adjustment of the metabolic response of the cell. To assess whether the observed metabolic response can indeed be explained by the measured proteome changes, combinatorial enzyme selection and regularization was used to identify key parameter adaptations. Next to this, the effect of individual iso-enzymes was considered and evaluated as a factor influencing the observed metabolic response. Furthermore, the kinetics and implementation of the storage metabolism were evaluated.

Carbon Storage Physiology Differs between Continuous and Dynamic Substrate Conditions
A kinetic model of yeast glycolysis was previously developed to fit various SS and SP datasets [29]. For the dynamic substrate conditions, the model was extended with regard to trehalose metabolism in different compartments and glycogen synthesis and degradation ( Figure 3).
The model simulations could reproduce most of the experimental data properly after several kinetic constants from HXT and HXK were estimated (this is explained in detail in the next section), but with a few exceptions ( Figure 4A). For instance, glucose 6-phosphate (G6P) and fructose 6-phosphate (F6P) simulated concentrations were smaller, which was already documented in the original model and attributed to underdetermined phosphofructokinase (PFK) reaction kinetics [29]. Other metabolites, such as fructose bisphosphate (FBP), glucose 1-phosphate (G1P) or trehalose 6-phosphate (T6P), also deviated. This could be explained by affinity constants undergoing changes during substrate perturbation cycles, here unaccounted for. Furthermore, reaction rates, estimated from 13 C enrichment data, were in close agreement ( Figure 4B), except for the maximum simulated rate for HXT with a simulated maximum lower than observed experimentally.   Figure 3. Kinetic metabolic model with a detailed description of trehalose and glycogen metabolism. Sink reactions account for fluxes towards the TCA, PPP and biomass synthesis. This model was adapted from [29]. The diagram style was adapted from [9]. For a detailed view on model mass balances and reaction kinetics, see Supplementary File S1.  In the substrate perturbation simulations, the increase in residual glucose during feeding resulted in transient changes in glycolytic metabolites, which returned to the initial value by the end of the cycle ( Figure 4C). Upper glycolysis metabolites, except glucose, reached their maximum concentration within 50 s due to recirculation via the storage metabolism. Glycerol branch and storage kinetics followed a similar pattern, but with a delay. Due to the slower reaction rate for enzyme GAPDH [55], the entry in lower glycolysis was delayed, and BPG reached its maximum at about 130 s. Nonetheless, the increase in FBP activated pyruvate kinase (PYK), reducing concentrations of 3-phosphoglycerate (P3G), 2phosphoglycerate (P2G) and phosphoenolpyruvate (PEP). As FBP decreased, its activation dissipated, and these lower glycolysis metabolites reached maximum concentrations in about 230 to 250 s of the cycle. This trend is in agreement with the known allosteric regulation of FBP on PTK [12].
Trehalose metabolic dynamics were different between repeated substrate perturbations and SP. During the repeated substrate perturbations, the maximum flux towards production of trehalose was less than 10% of the HXK reaction rate, in comparison to the 30% observed in SP [9], and part was secreted to the extracellular space, what is commonly regarded as stress protection [56,57]. Nonetheless, glycogen took up a greater portion during the repeated substrate perturbation regime, implying that carbon storage is predominant over the stress response by the trehalose cycle and suggesting that the cell is indeed adapted to the repeated substrate perturbation setup.
Small changes in protein expression occur between cells in a population. One way to examine the possible resulting phenotypes of the network upon perturbation is by means of ensemble modelling [58,59]. To test the robustness of the model, 10000 simulations were performed with random parameter values deviating within a range of 10% of the model parameter set. The concentration and reaction rate profiles were very consistent (see Supplementary Figures S1 and S2, respectively), especially for reaction rates, where the relative deviation between fluxes was very small. This suggested that the model dynamics are consistent within the parameter range tested.

Glucose Transport and Phosphorylation Identified as Key Adaptations from Combinatorial Parameter Estimation
Glycolytic enzyme expression changed from chemostat to dynamic substrate conditions, most notably for HXT and HXK ( Figure 2). As each iso-enzyme has specific kinetic properties, the catalytic (K cat ) and Michaelis-Menten (K M ) constants also differ [22]. To identify changes in kinetic parameters, parameter estimation based on the substrate perturbation datasets was performed and interpreted in light of the measured proteome changes. However, estimating all kinetic parameters simultaneously can lead to multiple local minima and ill-conditioning [34]. To bypass this problem and identify which are the key parameters that change between the cells adapted to continuous and dynamic substrate conditions, respectively, we adapted the scale and setup of the parameter estimation problem. Two stages were applied ( Figure 5): (1) Parameters were estimated for multiple combination of enzymes in parallel assays to isolate which enzymes were key to reproduce the data properly. (2) Regularization was implemented, i.e., parameters were re-estimated for the best selection of enzymes found and with a penalty for deviation from the reference parameter set.
In Step 1 ( Figure 5A), good fits were achieved when HXT, HXK/GLK and the trehalose cycle were included in the parameter estimation, suggesting that these are the relevant enzymes undergoing changes. In Step 2, an optimal fit between model error and parameter deviation was achieved by adding a regularization factor, shown in ( Figure 5B) for GLK kinetics. This overcame dependencies between kinetic constants of the reaction and pointed to K m,GLC , K i,T6P and V max being the key parameter alterations with respect to the reference parameter set from [29] (Figure 5C). Compared with other toolboxes available to perform parameter estimation in complex kinetic metabolic models [22,34,36,50,60], this pipeline incorporates regularization for known in vitro parameter values into a combinatorial enzyme selection approach, with the added value that intracellular flux data is used. Flux data have been available for only a decade, and few works have used it for yeast kinetic model development and validation [9,29].
As a result of this pipeline, some kinetic constant changes were suggested for the glucose transport and phosphorylation reactions (Table 1). For HXT, the maximal reaction rate (V max ) decreased from 8.13 to 1.7 mM s −1 , which is actually close to the value in other published models [27,28] and consistent with the experimental HXT concentration decrease (Figure 2). Changes in HXT isoenzyme proportions could also lead to changes in affinity; in vitro K M measurements have shown a high variability for the HXT1-7 subunits [21,61]. Nevertheless, the shift in isoenzymes here did not lead to drastic changes in affinity; the K M only changed from 1.01 to 0.90 mM. Parameters are regularized such that the data is still well reproduced (see Supplementary Figures S3 and S4). (C) change in key GLK parameters identified upon regularization; the deviation between the estimated parameter and the initial value taken from (bioRxiv2022) is shown in the y-axis (in logarithmic scale). Black and empty circles show the estimates prior and post regularization when parameter dependencies are minimized. The x-axis shows specific parameters. For the hexose phosphorylation reaction, several kinetic constants changed. Km glucose decreased from 0.35 to 0.11 mM, in line with the experimental increase in GLK/HXK ratio (Figure 2), given that the affinity constant for glucose was found lower in GLK than HXK both in vitro-and in vivo-like conditions [22,62]. Furthermore, the T6P inhibition constant (K i ) increased from 0.0073 to 0.0183, and K i,T6P was also found to change between the two isoenzymes [53]. Next to these parameters, V max increased from 6.25 to 15.75 mM s −1 , which is higher than values reported before in yeast glycolytic models [27,28]. Since the k cat for GLK is much lower than for HXK [22], such an increase in V max was not expected.
Outside glycolysis, changes were required in the trehalose cycle and ATP maintenance reaction. Note that the original model did not include glycogen metabolism, nor compartmentation of the trehalose cycle reactions. The missing reactions were added in this work, and the trehalose cycle parameters were re-estimated to account for the effect of the previously lumped reactions. Finally, the ATPase maximum reaction rate decreased to fit the adenosine nucleotide concentrations (ATP + ADP + AMP). This might be related to the fact that the initial model simulated the response to a GP of 20 g L −1 of glucose [29], which is seen by the cell as a stress condition due to the rapid increase in extracellular glucose concentration [9]. The final parameter set can be found in the Supplementary File S1.

Glucose Sensing Influences Hexose Transporter Kinetics during Substrate Perturbation Cycles
Glucose uptake has been widely modelled as an equilibrium-driven passive transport reaction [22,27], where its kinetics are determined by isoenzyme-specific V max and K M parameters [21,62,63]. Here, we have found that these kinetics alone cannot explain the experimental data, for which a glucose sensing mechanism [23] needs to be active.
By sampling the parameter space and using passive transport reaction kinetics, we found that no combination of parameters could fit the data ( Figure 6). Especially, none of the generated models could reproduce a net uptake reaction of almost zero at the end of the cycle (400 s) and reach the value of 0.72 mM s −1 at 20 s when the uptake rate reaches its maximum ( Figure 6A). Since throughout the entire substrate perturbation cycle the residual and maximum glucose concentration are 0.1 and 0.45 g L −1 , respectively, adjusting parameters to lower the effect of the transmembrane glucose gradient for one also reduces for the other. Interestingly, the only way to reproduce the experimental uptake ( Figure 6B) was by including a minimum glucose concentration term in HXT kinetics (Csmin in Equation (7)). This term acted as a threshold value when glucose import occurs.
where GLC ec2 = GLC ec − GLC ec,min (7) if GLC ec,min < GLC ec (8) Suarez-Mendez et al., (2014) already noticed a similar phenomenon when modelling glucose uptake dynamics. Here, we assumed that the threshold was an effect of glucose sensing under certain conditions [23]. Glucose sensing acts independently of glucose uptake [64] and is known to activate a cascade of reactions and ultimately lead to altered gene expression in yeast [65][66][67]. This could imply that, in line with [68], glucose sensing is observable under the repeated substrate perturbation condition but not in the SP experiments in which residual glucose concentrations are remarkably higher [9].

13C-Labelled Metabolite Mass Balances Validate the Model but Suggest Caveats in Carbohydrate Storage Metabolism
Metabolic and flux profiles agreed between simulations and experimental data ( Figure 4). Even though 13C isotope labelling was used in the flux estimation [26], these data have not yet been implemented in kinetic models. Here, we aimed at validating the model by implementing individual mass balances for each labelled metabolite (carbon structures) in the network. We found a considerable degree of agreement when simulating enrichment profiles. For the first 100 s of the cycle, the percentage of enriched metabolite rose to about 80% for most metabolites ( Figure 7A), It then decayed as recirculation of unlabelled trehalose and glycogen became more prominent ( Figure 7B-D).   Conversely, the accuracy of model simulations was also limited. Enrichment of T6P decreased slower and glycolytic metabolites faster than expected during the late cycle, which could indicate that there is a surplus of glycogen recirculation which was mostly unlabelled. This might be explained by the current glycogen metabolism kinetics, which were simplistic here. Small deviations from the experimental value can have a great impact on the late cycle stage, given that fluxes in the network are generally low. We initially aimed at representing glycogen synthesis and degradation as mass action or Michaelis-Menten kinetics, but unfortunately, this did not resemble the experimental reaction rates due to the high and relatively constant glycogen concentrations. Therefore, simplified phenomenological expressions were used (see details in Supplementary File S1). Besides glycogen metabolism, enrichment of lower glycolysis metabolites was faster than expected, which could be attributed to the changes observed in other isoenzymes such as TDH (Figure 2).

Missing Regulation in Trehalose Metabolism
Another location with uncertainty is the trehalase reaction, which is carried out by an acid and neutral enzyme (ATH1 and NTH1, respectively) [69,70] and whose in vivo fluxes were quantified in the repeated substrate perturbation condition in [26]. In the model simulations, NTH1 trehalase activity was reproduced but only if cytosolic trehalose concentration was artificially low (Figure 8A-C) and redirected to the other compartments. Nonetheless, trehalose is expected to locate more in the cytosol than the vacuole [71]. This occurred as a result of NTH1 reaction being modelled as simple Michaelis-Menten kinetics ( Figure 8A) [46]. To fit these kinetics, cytosolic concentrations were kept very low with a comparatively high increase during the cycle. To further doubt the current model understanding, the estimated Km.TRE decreased from 2.11 to 0.13 mM, but it was experimentally quantified to be 3-8 mM [72]. We are uncertain is what this missing regulation could be. A post-translational regulation acting on NTH1 could be a possible explanation [32,72], but new data on the state of the enzyme would be required to confirm this claim.

Conclusions and Summary
In this work, we combine metabolome, fluxome and proteome data to develop a metabolic model. Based on detailed analyses, we suggest enzymatic reactions whose reaction kinetics adapt to the dynamic substrate conditions and locations where our knowledge is limited. Testing different subsets of parameters for recalibration highlighted transporters and phosphorylation reactions as crucial for the adaptation. This in silico approach is comparable to the experimental approach of metabolic reverse engineering [73] but much faster and less laborious as no experiments with combinatorial genome modifications are required. The combinatorial approach can also be applied to other industrially relevant downscaling setups which are relevant to finding out key parameter changes in a relatively simple manner and further optimize the bioprocess.
Furthermore, disagreements between experimental and simulated data suggest that the assumed mechanistic kinetics cannot sufficiently describe the intracellular flux and metabolome. Here, glucose uptake could not be explained by facilitated diffusion only (Equation (7)) but required a glucose threshold concentration (Figure 4). In addition, some reactions of the storage metabolism required non-mechanistic adjustments to reproduce the observed labelling enrichments.

Author Summary
Kinetic metabolic models are used to understand how biological systems deal with dynamic perturbations in their environment. A well-known case of their application is the microorganism Saccharomyces cerevisiae, which was domesticated by mankind thousands of years ago, and is used to produce a wide range of products, such as bread, beverages and biofuels. When cultured in industrial-scale bioreactors, this cell factory is impacted by environmental perturbations which can challenge the bioprocess performance. The repeated substrate perturbation regime has been proposed as an experimental setup to downscale these industrial perturbations. Intracellularly, these perturbations impact central carbon metabolism, including carbon storage. Even though kinetic metabolic models have been developed to study the effect of single extracellular perturbations, they have not explored repeated substrate perturbations and their implications on carbon metabolism. We developed a model construction and parameter identification pipeline and used it to expand the existing models to represent carbon metabolism under dynamic substrate conditions. We used computer simulations to point at adaptations in yeast metabolism and locations in the model where our understanding is not entirely accurate. We found that combining multiple types of data, despite being challenging, can be very beneficial by providing a comprehensive and realistic representation of the cell.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/metabo13010088/s1, Figure S1: Simulation of metabolic concentrations is robust to parameter changes within 10% of their estimated value; Figure S2: Simulation of reaction rates is robust to parameter changes within 10% of their estimated value; Figure S3: Metabolite concentrations (mM) over time (s). Model fit; Figure S4: Reaction rates (mM s −1 ) over time (s). Model fit; Figure S5: Glucose sensing is needed to explain HXT kinetics: Individual parameter effect; File S1: Ordinary differential equations, Rate equations, Parameter set, Summary rate equations and Simulation setup.  Data Availability Statement: The data used in this work can be found in the github repository of the author (github.com/DavidLaoM/y3m2_ff, accessed on 3 January 2022). The experimental proteome dataset can be found in [51].