Model Parameterization with Quantitative Proteomics: Case Study with Trehalose Metabolism in Saccharomyces cerevisiae

Abstract

When Saccharomyces cerevisiae undergoes heat stress it stimulates several changes that are necessary for its survival, notably in carbon metabolism. Notable changes include increase in trehalose production and glycolytic flux. The increase in glycolytic flux has been postulated to be due to the regulatory effects in upper glycolysis, but this has not been confirmed. Additionally, trehalose is a useful industrial compound for its protective properties. A model of trehalose metabolism in S. cerevisiae was constructed using Convenient Modeller, a software that uses a combination of convenience kinetics and a genetic algorithm. The model was parameterized with quantitative omics under standard conditions and validated using data collected under heat stress conditions. The completed model was used to show that feedforward activation of pyruvate kinase by fructose 1,6-bisphosphate during heat stress contributes to the increase in metabolic flux. We were also able to demonstrate in silico that overexpression of enzymes involved in production and degradation of trehalose can lead to higher trehalose yield in the cell. By integrating quantitative proteomics with metabolic modelling, we were able to confirm that the flux increase in trehalose metabolic pathways during heat stress is due to regulatory effects and not purely changes in enzyme expression. The overexpression of enzymes involved in trehalose metabolism is a potential approach to be exploited for trehalose production without need for increasing temperature.

1. Introduction

Trehalose has often been associated with the eukaryotic model organism baker’s yeast during heat stress, as the microbe is observed to accumulate high concentrations of this protective molecule for survival [1,2,3]. The protective nature of trehalose is due to it being a stable and generally unreactive sugar, acting as a robust energy storage vehicle [4]. Owing to the aforementioned protective properties, trehalose has high commercial value as it is used in various industries from pharmaceuticals to food and cosmetics [4,5,6]. In the pharmaceutical industry trehalose has been used for the storage of vaccine at room temperature as it has been found to stabilise vaccines. While in cosmetics, trehalose serves as liposome stabilizer and as stable sweetener for food. It has also recently been found to offer protection against bone loss in mice [7]. Its production in the industry relies on the use of enzymes from extremophiles expressed in other microbes [6].

In addition to the increase in flux towards trehalose production, another phenomenon observed during heat stress adaptation in S. cerevisiae is the increase in glycolytic flux [2,8]. Postmus et al. [8] set out to investigate the possible factors that could contribute to the increase in flux, including gene expression, enzyme activity, protein expression, and the profile of metabolites. Interestingly, their data indicated minimal changes in gene and protein expression, but extensive changes in metabolic profile. Therefore, they attributed the changes in flux to the augmentations in the metabolic environment of the enzymes, such as the close to 10-fold increase in fructose 1,6-bisphosphate and 20-fold increase in pyruvate. Additionally, they postulated that the maintenance of the high flux is a result of the feedback activation of phosphofructokinase by fructose 2,6-bisphosphate and feedforward activation on pyruvate kinase by fructose 1,6-bisphosphate.

The development of mathematical models to aid metabolic engineering is not uncommon for biologists nowadays, with various frameworks and approaches available for building the models [9,10]. However kinetic models are often limited to small or medium sized models due to the amount of information needed to construct them, ranging from rate laws for the reactions to the kinetic parameters. In order to circumvent this, we developed Convenient Modeller [11], which uses a generalized Michaelis-Menten equation in the form of convenience kinetics [12] and parameter estimation to fill in the missing kinetic parameters. Kinetic parameters, which form a vital part of kinetic models, are constant values in the equations. They can be measured experimentally in the laboratory, and some of them are stored in databases such as BRENDA [13] and SABIO-RK [14]. Convenient Modeller provides users the option to circumvent the process of determining individual kinetic parameters by wet lab experiments.

In this study, a kinetic model consisting of trehalose metabolism in Saccharomyces cerevisiae as well as the upper portion of glycolysis was developed, using only convenience kinetics to form the rate equations. This model was fitted to steady state data of metabolites and fluxes. On top of that, the model uses quantitative proteomics data collected under two different conditions, for training and validation. Protein data from standard conditions were used to parameterize the model, while separate data collected under 37 °C were used to help simulate the biological condition of heat stress and validate the model. Not to be confused with the more extreme heat shock, this heat stress is at a more physiological temperature of approximately 37 °C [15].

To investigate if the regulatory effects in the glycolytic system play a role in increasing the flux during heat stress, a kinetic model of trehalose metabolism in S. cerevisiae was constructed. This model should also help to determine the best way to produce trehalose without the need for temperature increase. The model was built using metabolomics, fluxomics and proteomics data collected under standard conditions and subsequently validated using metabolomics and proteomics data collected under heat stress condition. The completed model was subjected to in silico regulation analysis and overexpression study. It was found that activation of pyruvate kinase does contribute slightly to the control of flux in the trehalose cycle and glycolysis, and that regulatory effects on enzymes involved in glucose entry play a significant role in affecting flux within the system. The model additionally predicted that the overexpression of enzymes directly involved in the production and degradation of trehalose would lead to an increase in its concentration in the system. 2. Materials and Methods 2.1. Model Construction

The model was built using Convenient Modeller (https://github.com/chuanfuyap/Convenient-Modeller, accessed on 15 May 2020) [11], which uses convenience kinetics (detailed below) [12] to build the rate equations. Convenient Modeller is a tool for building kinetic models for cell metabolism, requiring input from the users on the substrates and products involved in the metabolic reactions and the enzymes that catalyse the reactions, which are then translated into rate equations automatically. The convenience kinetics assumes that all the reactions are reversible and have a random binding order. A genetic algorithm (detailed below) is used to estimate the kinetic parameters given metabolomics, fluxomics and/or proteomics data. In brief, given information on substrates, products and enzymes, the tool can generate the rate equations that constitute the model; and by providing quantitative information such as metabolite concentrations, protein concentration and flux values, the kinetic parameters of the model can be estimated. Models built in Convenient Modeller use the Systems Biology Markup Language format [16].

2.2. Enzyme Kinetics

For the construction of models, Convenient Modeller makes use of convenience kinetics, a generalised reversible Michaelis-Menten equation. Specific rate laws for all reactions in the system are not essential, because their specific properties tend to be dissipated in a large system.

v(sub, prod)=Etotal.fregk+cat ∏i subi˜ni−k−cat ∏j prodj˜nj∏i(1+subi˜+…+subi˜ni)+∏j(1+prodj˜+…+prodj˜nj)−1,

where sub is the substrate concentration; prod is the product concentration;_Etotalis the enzyme concentration;_fregis a pre-factor to account for activation, using 1 + d/^kA, d is the activator concentration;^kAis the activation constant; or inhibition, using^kI/^kI+d, d is the inhibitor concentration,^kIis the inhibition constant, where activator and inhibitor are metabolites that increase or decrease the reaction rate;_k+/−catare the forward and reverse turnover rates;sub˜= sub/_ksubM;prod˜= prod/_kprodM;_ksub/prodMare the Michaelis-Menten constants for either substrate or product; n is the stoichiometric coefficient for the reaction.

2.3. Parameter Estimation

After establishing the system’s network of reactions, the kinetic parameters are needed to complete the model. All the kinetic parameters of the model were estimated using a genetic algorithm [17].

Genetic algorithm (GA, Figure 1) is a machine learning algorithm used in optimisation problems inspired by the Darwinian evolutionary principle. In brief, GA ‘replicates’ the evolutionary process with a population of potential solutions (the kinetic parameters) that go through multiple generations of mutations, selection and reproduction. The potential solutions are evaluated at every generation with a fitness function (this is GA terminology for objective function), and the population evolves over generations to produce the best solution. The fitness function of choice here is Mean Absolute Percentage Error.

Σ|Observed−SimulatedObserved|·100%n,

where observed is the training data’s value for metabolite and/or flux, simulated is there model simulated value for metabolite and/or flux and n is the total number of metabolites and fluxes in the model.

Unlike most GA, binary encoding is not used for chromosomes, instead every single kinetic parameter in the chromosome is encoded as an exponent for a base value of 10 (decadic logarithm value). The reason for this encoding is that exact values for a kinetic parameter is not a priority, rather the orders of magnitude. Parameter values are generated randomly for the first generation. To spawn the next generation, random chromosomes are chosen as parents using tournament selection, in which the population is randomly sampled several times (depending on the population size) to select for fit parents. After parents are selected, crossover is done using uniform crossover to produce offspring for the next generation. In this version of uniform crossover, a fixed ratio of parameters is exchanged instead of exchange at fixed points. On top of crossover, random chromosomes are chosen to be mutated at every generation as well. Mutation is carried out by randomly increasing or decreasing one or several parameters in the chromosome. A plague function is also written into this GA, whereby after a set amount of generations determined by the user, the population would get reduced to the initial size, purging chromosomes with low fitness score. Users can also set a plateau limit, which is the maximum number of generations for when the highest fitness score stays the same. After reaching the plateau limit, there is an adaptive function in the GA where every chromosome is mutated to try and escape this plateau, if it fails, the search comes to an end, spawning a model from the fittest chromosome. The ending criteria in a perfect scenario would be when a chromosome is able to perfectly replicate the fitting data as its output, resulting in zero error score. For the model presented in this paper, 106 parameters were estimated with the GA, which took about 130 h to run on a computer cluster with 24 cores. 2.4. Fitting Data

The fitting data for the metabolites were obtained from two separate studies, one on S. cerevisiae in normal conditions measured using mass spectrometry [18] and another measured using nuclear magnetic resonance [19]. In this work, normal or standard condition refers to culture of yeast under 30 °C with the use of complete culture medium, such as yeast extract peptone dextrose (YPD) or yeast nitrogen base which can be supplemented with amino acids. Although the work done by Puig-Castellví et al. [19] used arbitrary units for the metabolites, the values of shared metabolites were not too dissimilar from those measured by Smallbone et al., [18] so they were treated as mM. The fluxes were obtained from data from the reference strain measured by Blank and colleagues [20]. All of the data used for parameter estimation were values measured during steady state.

The protein concentrations used in the model (for both parameter estimation and validation) were obtained from a quantitative and temporal study on S. cerevisiae’s proteins undergoing heat stress [21].

The fitting data are from different sources, as a result of being generated with different strains and media conditions and culture methods. The S. cerevisiae strains used are Y23925, BY4741, CEN.PK 113-7D and BY4742 for Smallbone et al., Puig-Castellví et al., Blank et al. and Jarnuczak et al. respectively. The media condition used are YPD medium, yeast minimal medium and Yeast Nitrogen Base medium for Puig-Castellví et al., Blank et al. and Jarnuczak et al., respectively. The medium used by Smallbone et al. was not specified other than that it was without nutrient limitation. Smallbone et al. used continuous culture, while the other three studies used batch cultures. The differences mentioned above are phenotypically neutral and should have no significant effect on the trehalose metabolism modelled here. 2.5. Manipulation and Simulation of Model In order to validate the model, the enzyme values of the model fitted under standard conditions were changed to those of the new condition, which are the enzyme values under heat stress condition. The overexpression was simulated by doubling the enzyme of interest’s value from its standard condition value.

The fitted model was modified using the JSBML library [22] with Java code to make the changes needed for validation, regulation analysis and overexpression model interrogation; the changes include altering concentration values of the enzymes, changes to the equations and kinetic parameters involved, and in the case of regulation analysis removal of activator and inhibitor constants. The simulation of the models up to steady state was done using the SBML ODE Solver Library [23].

3. Results Heat stress on yeast is a phenomenon where yeast cells are introduced to an increase in temperature up to 37 to 38 °C, beyond their optimal conditions, leading to multitude of changes in their physiology including metabolically. Notably, this includes an increase in the cytoprotective metabolite trehalose. Of these changes, the increase in glycolytic flux is of interest in this study as the precise cause is yet to be determined. Additionally, trehalose is a molecule of commercial interest, so theoretical exploration is also performed to determine ways to increase its production in yeast.

A kinetic model focusing on trehalose metabolism parameterized with quantitative metabolite concentrations, flux, and protein abundance data collected before heat stress was applied to the cells, is presented here. The model was validated using a separate set of quantitative proteomics from the same study [21], but collected from cells undergoing heat stress, allowing the model to simulate heat stress. The validated model was used to study the possible causes of flux increase during heat stress as well as ways of increasing trehalose production in S. cerevisiae without the need for an increase in temperature.

3.1. Model of Trehalose Cycle & Upper Glycolysis

The completed model in this study included S. cerevisiae’s upper portion of glycolysis, and trehalose metabolism (File S1) [24] (Figure 2). Trehalose metabolism is a cycle that branches off from glucose-6-phosphate of glycolysis and re-enters glycolysis as trehalose, broken down into two glucose molecules. The main aim of this model is to study trehalose metabolism, therefore the lower glycolysis reactions from glyceraldehyde 3-phosphate to pyruvate have been simplified to a single reaction from five original reactions, which is here represented as modified pyruvate kinase reaction (PYK mod). This change has no significant effect on the trehalose metabolism as the same flux was maintained, and rather than combining the metabolites in the reactions into one which includes nicotinamide adenine dinucleotide (NAD) and additional adenosine triphosphate (ATP), they were excluded as the pyruvate kinase reaction does not interact directly with them. In total the model contains 23 metabolites (6 external metabolites, boundary condition true in SBML), 20 enzymes and reactions. Within the reaction network there are three activation regulatory effects and four inhibitory effects [25]. All these summed up to 106 kinetic parameters in the model that are fitted using omics data (Table S1).

3.2. Model Fitting and Validation

The model was fitted using machine learning via a genetic algorithm, with heterogeneous data for metabolites (Table 1) and fluxes values (Table 2) (sources discussed in Methods) to obtain the kinetic parameters. As the data source for flux does not have values for the trehalose flux, we have assumed it to be 1% of the flux towards it through the T6PS reaction, which was rounded up to be 0.2 mmol/g/h. Following the parameter estimation step, the enzyme concentrations in the model were all replaced with those measured under heat stress condition to simulate the condition. A single source for protein values was used; values measured under 30 °C were used for fitting, while those measured under heat stress at 37 °C were used for validation; the values were taken from the final timepoint of the dataset, 240 min (Table S2) [21]. Reassuringly, the model was able to successfully replicate the main metabolic responses observed in heat stress for S. cerevisiae [2,8,19,26,27], namely increase in overall fluxes in the network and high accumulation of trehalose (Table 1 and Table 2).

The root-mean-square error of the model’s simulated values and fitting data at 30 °C is 1.56, with the majority of the deviations contributed by fluxes (Table 2). The differences between simulated and experimentally measured values for fluxes at 30 °C are all below 13% except for alcohol dehydrogenase with a deviation of 22%. For the metabolites, only ATP deviated significantly (Table 1). Qualitatively, changes in the direction for the metabolites are in the right direction for the majority of metabolites considered, with the exception of trehalose-6-phosphate (T6P), glucose 6-phosphate (G6P) and fructose 6-phosphate (F6P) (Figure 3).

3.3. Regulation Analysis In order to determine if modifiers play a role in increasing the overall fluxes during heat stress, regulation analysis was performed on the validated model. This was done by removing individual or groups of enzyme modifiers in the network of reactions while maintaining everything else and simulating the model until it reached a steady state. Data for removal of glucose transport inhibition is not shown because its removal leads to an unstable model that is unable to achieve steady state.

Of all the reactions investigated, the removal of modifiers involved in the upper part of glycolysis (GLT inhibition, HXK inhibition, T6PS activation, GSY activation and GPH inhibition) caused significant changes in fluxes relative to the original model’s flux values when protein concentrations collected under 37 °C were used (Figure 4, Table S3). When there is no inhibition in upper glycolysis, the trehalose cycle’s fluxes see a drop, along with a massive decrease in the reaction where glycogen is broken down into glucose-1-phosphate, catalysed by glycogen phosphorylase and with a huge increase in flux towards the pentose phosphate pathway. When the activation of pyruvate kinase is removed, there is a slight increase in flux within the trehalose cycle, but with drops observed in the entry to the trehalose cycle and lower glycolysis.

3.4. Overexpression Simulations and Its Yields

To determine the best target for metabolic engineering in S. cerevisiae in order to achieve higher trehalose production, overexpression simulations (doubling the enzymes’ concentration) were performed with the model. Overexpression was chosen as the in silico modification of choice as this is a proven method in increasing metabolite production in metabolic models [26,27]. In order to narrow down the search of enzymes to overexpress, metabolic control analysis [28,29] was performed. It suggested that UDP-glucose phosphorylase and trehalose-6-phosphate synthase increase would result in higher trehalose accumulation. That was proved incorrect for UDP-glucose phosphorylase (Figure 5), but true for trehalose-6-phosphate synthase. However, overexpression for the trio of enzymes involved in trehalose production and degradation of trehalose delivered the best results in trehalose accumulation with 10000% increase in yield (Figure 5). Yield here is calculated by dividing the trehalose value at steady state by the extracellular glucose value.

4. Discussion

The chemically unreactive and stable disaccharide trehalose is found in variety of organisms including yeast. Its chemical properties provide protective effects to microorganisms allowing it to survive stressful situations. The same survival feature is exploited in the biotechnology industry allowing it to serve as protectant for storage of food and in pharmaceutical industries. This has brought about a lot of interest in the molecule, which resulted in the development of kinetic models studying trehalose metabolism in S. cerevisiae [25,30]. Voit’s work was used to determine metabolic regulation involved in the trehalose metabolism [25], while Fonseca and colleagues developed a model using time series of metabolic profiles to predict changes in protein levels [30]. The two models mentioned above, and the model presented in this work, all contain the complete trehalose cycle including the regulatory effects. The difference starts from the glucose-6-phosphate, where our model conserves much of the lower glycolysis maintaining three exit fluxes; the model in Voit’s work removes the entire lower glycolysis and treats it as exit flux, while Fonseca’s model has two exit fluxes with ethanol, acetate and glycerol production in one reaction and a leakage flux for carbon used in TCA cycle and respiration. Aside from the network, from the technical standpoint the other two models use a different modelling framework which uses simplified rate action kinetics; Voit employed S-system equations and Fonseca used generalized mass action equations. While these kinetic representations can be used to describe the reactions, they are less accurate biochemically as they do not account for enzyme saturation, a biochemically relevant phenomenon, which is accounted for with convenience kinetics used in this study.

Our model has used quantitative protein data to inform kinetic parameters as well as to emulate the effects of heat stress in the model. The estimation of unknown parameters was done using a genetic algorithm as this is commonly used for metabolic modelling [31,32,33]. The use of proteomics over transcriptomics is a better approach for heat stress simulation as it was found that transcription changes are mostly transient for stress adaptations in S. cerevisiae [34] but proteins remain longer in steady state to adapt to environmental stress [21]. The model built here was done by estimating all the parameters in the model, however this approach can become computationally expensive as the model scales up. Other researchers in the systems biology modelling community have developed methods that select only subsets of parameters to be estimated [35,36,37]. These methods involve the use of sensitivity analysis to identify the parameters that have the biggest influence on the system, which would form the subset used in the estimation procedure. This idea can be adapted to the Convenient Modeller framework in order to reduce the computation time needed, where the parameters not estimated can be given a default value of one. An alternative approach to model building is to use experimentally measured parameters directly in the model, or to make use of such values as estimation boundary during the estimation step. This can also be achieved in Convenient Modeller, however the work here did not aim to reproduce the parameters measured experimentally, which is usually done in isolation of other reactions. Instead, our model was generated as a holistic system where all the parameters interplay with one and other. This also follows the concept of “sloppy model”, discussed by Brown et al. [38,39] showing that a model remains useful despite some parameters varying from the “true” values.

The combination of convenience kinetics with heterogeneous sources of fitting data have been proven to be an effective modeling method here, as it produced a model with a good fit to the training data. Steady state data was used as there was insufficient time-series data for all the metabolites and fluxes in the trehalose metabolic pathway; should such dataset be produced, it can be used in conjunction with the temporal proteomics data to further improve the model. Despite using only steady state data the model was also capable of reproducing the metabolic phenomenon commonly observed during heat stress adaptation (Table 1 and Table 2). When S. cerevisiae is challenged with high temperature, it is commonly observed that fluxes would increase [2,8] and trehalose would accumulate [3,19,27]. When protein amounts collected during heat stress were input, the model was able to reproduce the increase in flux generated during heat stress, albeit to a lesser degree than those experimentally measured (Table 2). Additionally, it reproduced increased glucose consumption, which is indicated by the drop in glucose level and increased flux in glucose transport [2], as the cell would need to spend more energy to cope with the stress. Minimal changes in trehalose 6-phosphate were also predicted by the model, which is expected as this compound is toxic for the cell in high concentration [40]. It is worth noting though that predictions for glucose 6-phosphate, fructose 6-phosphate and trehalose-6-phosphate were in the opposite direction of those measured by Postmus et al. [8]. The increase in flux as a result of heat stress is at least five-fold [8]; however, our model only predicted a 40% increase in flux for most of the reactions. This discrepancy might be the result of difference in the experimental conditions between the Postmus et al. study [8], which measured the flux values, and the Jarnuczak et al. study [21], which determined the protein concentrations used in this study. Postmus et al. used naïve cells grown to stationary phase to measure fluxes, while Jarnuczak et al. used heat adapted cells to measure protein concentrations, which were generated by moving them to 37 °C during the mid-exponential growth phase. This results in different genomic changes in the cells to better adapt to the increased temperature [30].

Our model was additionally used to test the hypothesis that regulatory effects in the metabolic network contribute to flux increase. While we did not include the activation of phosphofructokinase by fructose 2,6-bisphosphate, we included the feedforward activation on pyruvate kinase by fructose 1,6-bisphosphate. The removal of this activation on pyruvate kinase resulted in drop of flux in lower glycolysis (Figure 4), confirming the hypothesis in silico. However, this also resulted in a slight increase in flux in the trehalose cycle, channeling flux from glycogen instead of upper glycolysis. The regulatory interactions that have the highest impact on fluxes of glycolysis and trehalose cycle are in upper glycolysis. When all three inhibitions on upper glycolysis (glucose transport, hexokinase and phosphofructokinase) were removed, this resulted in a decrease flux within the trehalose cycle with a corresponding higher flux channeled towards the pentose phosphate pathway. Such a phenomenon if it occurred in a living cell would not be sustainable and the cell would not survive the heat stress since it would be unable to generate the trehalose needed for protection. This demonstrates the importance of regulatory control in metabolic pathways.

The prediction made in the enzyme overexpression investigation leads to a high increase in percentage yield of trehalose production (Figure 5), which was also observed experimentally. In the work done by Fonseca et al. [30] it was shown that heat adapted yeast cells have trehalose concentrations that increased from 4 mM up to 100 mM. Additionally, work by Magalhães et al. showed a 5-9-fold increase in trehalose for heat adapted cells [41]. Furthermore, Mahmud et al. [42] have studied S. cerevisiae strains which contained three gene deletions but overexpressed either TPS1 or TPS2. The overexpressed genes were the ones that encode for trehalose 6-phosphatase; they have found that the overexpressed strains had higher trehalose content without stress than their non-overexpressed counterpart. This highlights the utility of systems modelling to show how metabolite levels can be potentially altered via changes in the underlying proteome.

5. Conclusions Our work shows that the use of convenience kinetic for metabolic modelling is an efficient solution for rapid prototyping of metabolic models for hypothesis generation. This methodology allows the use of multi-omics data for parameter estimation. Subsequently, enzyme concentrations can be updated in silico using quantitative proteomics collected from a different condition such as heat stress, which then allows the model to simulate this new condition. The output can then be compared with a different set of multi-omics data, constituting the validation step of the model.

Figure 1. Genetic algorithm decision flowchart utilised by Convenient Modeller.

Figure 2. Trehalose metabolism model, with lower glycolysis simplified to a single reaction from glyceraldehyde 3-phosphate to pyruvate, as well as skipping production of acetaldehyde to ethanol. GlcO: extracellular glucose, GlcI: intracellular glucose, G6P: glucose 6-phosphate, F6P: fructose 6-phosphate, F16BP: fructose 1,6-bisphosphate, DHAP: dihydroxyacetone phosphate, GAP: glyceraldehyde 3-phosphate, G1P: glucose 1-phosphate, UDPG: urine diphosphate glucose, T6P: trehalose 6-phosphate, GLT: glucose transport, HXK: hexokinase, PGI: phosphoglucose isomerase, PFK: phosphofructokinase, ALD: aldolase, TPI: triose-phosphate isomerase, PYK mod: pyruvate kinase, PYC: pyruvate carboxylase, ADH: alcohol dehydrogenase, PGM: phosphoglucomutase, UDPGP: UDPG phosphorylase, GSY: glycogen synthase, GPH: glycogen phosphorylase, T6PP: trehalose phosphatase, T6PS: trehalose-phosphate synthase, G6PDH: glucose 6-phosphate dehydrogenase.

Figure 3. Log2 fold changes of metabolites from 30 °C to 37 °C generated by the trehalose metabolism model. With comparison from data measured experimentally by Puig-Castellví et al. [19] for intracellular glucose (GlcI), trehalose, adenosine triphosphate (ATP), and nicotinamide adenine dinucleotide (NAD) at 37 °C and Postmus et al. [8] for glucose-6-phosphate (G6P), fructose-6-phosphate (F6P), fructose-1,6-bisphosphate (F16P), trehalose-6-phosphate (T6P), and adenosine diphosphate (ADP) at 38 °. Remaining metabolites had no reference experimental data for comparison.

Figure 4. Heat map showing log fold changes for fluxes between the original model and a series of modified models where modifiers in reactions are removed. Each column represents one modified model, each row represents one reaction flux. All these results were generated with the models undergoing heat stress.

Figure 5. Theoretical percentage increase in mol trehalose per mol glucose. Enzymes here are overexpressed two-fold.

Metabolites (mmol/L)	Fitting Data	30 °C	37 °C
GlcI *	6.277	6.39	1.76
G6P *	0.774	0.74	0.14
F6P *	0.235	0.23	1.19
F16P *	4.583	4.52	7.81
DHAP *	1.162	1.15	0.99
GAP *	0.316	0.31	0.36
Pyruvate *	2.107	2.01	4.23
G1P		10.34	5.24
UDPG		0.01	0.01
T6P		0.05	0.07
Trehalose **	0.33	0.33	1.91
ATP **	3.2	5.58	5.57
ADP		0.004	0.01
UDP		0.81	0.97
UTP		0.58	0.42
NAD		1.57	1.62
NADH		0.22	0.17

Blanks means the values were not fitted as they were not available. * metabolites use data from Smallbone et al. [18]. ** metabolites (along with ethanol, glycerol and TCA cycle represented by their by-product succinate and citrate, for the external metabolites) uses data from [19].

Fluxes (mmol/g/h)	Fitting Data	30 °C	37 °C
GLT	16.7	15.52	21.71
HXK	16.7	15.93	22.29
PGI	14.2	13.23	21.44
PFK	14.2	13.23	21.44
ALD	15.2	13.23	21.44
TPI	13.5	11.68	10.87
G6PDH	1.8	1.75	0.12
GlycerolDH	1.7	1.55	10.57
PYK mod	28.6	24.91	32.30
ADH	23.6	18.30	21.58
PYC	5.0	6.62	10.72
PGM	0.8	0.75	0.44
UDPGP		0.85	1.41
GSY		0.65	1.12
GPH	0.10	0.10	0.96
T6PS	0.20	0.20	0.29
T6PP	0.20	0.20	0.29
Trehalase	0.20	0.20	0.29
ATPSyn		5.10	12.83
UDPtoUTP		0.85	1.41

Supplementary Materials

The following are available online at https://www.mdpi.com/2227-9717/9/1/139/s1, Table S1: Kinetic Parameters Values, Table S2: Protein Numbers, Table S3: Regulation Analysis Flux Values output.

Author Contributions

Conceptualization, J.-M.S., S.J.H. and C.F.Y.; methodology, C.F.Y.; software, C.F.Y.; validation, C.F.Y.; formal analysis, C.F.Y.; investigation, A.F.J.; resources, A.F.J. and S.J.H.; data curation, C.F.Y., A.F.J. and M.G.-A.; writing-original draft preparation, C.F.Y.; writing-review and editing, C.F.Y., J.-M.S., S.J.H., A.F.J. and M.G.-A.; visualization, C.F.Y.; supervision, J.-M.S. and S.J.H.; project administration, J.-M.S.; funding acquisition, S.J.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Biotechnology and Biological Sciences Research Council (BBSRC) in the form of grants BB/M025748/1 (Ecoyeast) and BB/G530225/1 (Systems Biology Doctorial Training Centre).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in Table 1 and Table 2 & Table S1. The raw proteomics data of Table S2 are available in PRIDE [43] under identifier PXD006262.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Strassburg, K.; Walther, D.; Takahashi, H.; Kanaya, S.; Kopka, J. Dynamic transcriptional and metabolic responses in yeast adapting to temperature stress. OMICS 2010, 14, 249-259.

2. Mensonides, F.I.C.; Schuurmans, J.M.; Teixeira de Mattos, M.J.; Hellingwerf, K.J.; Brul, S. The metabolic response of Saccharomyces cerevisiae to continuous heat stress. Mol. Biol. Rep. 2002, 29, 103-106.

3. Parrou, J.L.; Teste, M.A.; François, J. Effects of various types of stress on the metabolism of reserve carbohydrates in Saccharomyces cerevisiae: Genetic evidence for a stress-induced recycling of glycogen and trehalose. Microbiology 1997, 143, 1891-1900.

4. Paiva, C.L.A.; Panek, A.D. Biotechnological Applications of the Disaccharide Trehalose. Biotechnol. Annu. Rev. 1996, 2, 293-314.

5. Roser, B. Trehalose, a new approach to premium dried foods. Trends Food Sci. Technol. 1991, 2, 166-169.

6. Schiraldi, C.; Di Lernia, I.; De Rosa, M. Trehalose production: Exploiting novel approaches. Trends Biotechnol. 2002, 20, 420-425.

7. Xu, X.; Wang, R.; Wu, R.; Yan, W.; Shi, T.; Jiang, Q.; Shi, D. Trehalose reduces bone loss in experimental biliary cirrhosis rats via ERK phosphorylation regulation by enhancing autophagosome formation. FASEB J. 2020.

8. Postmus, J.; Canelas, A.B.; Bouwman, J.; Bakker, B.M.; Van Gulik, W.; De Mattos, M.J.T.; Brul, S.; Smits, G.J. Quantitative analysis of the high temperature-induced glycolytic flux increase in Saccharomyces cerevisiae reveals dominant metabolic regulation. J. Biol. Chem. 2008, 283, 23524-23532.

9. Almquist, J.; Cvijovic, M.; Hatzimanikatis, V.; Nielsen, J.; Jirstrand, M. Kinetic models in industrial biotechnology-Improving cell factory performance. Metab. Eng. 2014, 24, 38-60.

10. Costa, R.S.; Hartmann, A.; Vinga, S. Kinetic modeling of cell metabolism for microbial production. J. Biotechnol. 2016, 219, 126-141.

11. Yap, C.F. Convenient Modeller. Available online: https://github.com/chuanfuyap/Convenient-Modeller (accessed on 15 May 2020).

12. Liebermeister, W.; Klipp, E. Bringing metabolic networks to life: Convenience rate law and thermodynamic constraints. Theor. Biol. Med. Model. 2006, 3, 41.

13. Schomburg, I.; Chang, A.; Ebeling, C.; Gremse, M.; Heldt, C.; Huhn, G.; Schomburg, D. BRENDA, the enzyme database: Updates and major new developments. Nucleic Acids Res. 2004, 32, D431-D433.

14. Wittig, U.; Kania, R.; Golebiewski, M.; Rey, M.; Shi, L.; Jong, L.; Algaa, E.; Weidemann, A.; Sauer-Danzwith, H.; Mir, S.; et al. SABIO-RK-database for biochemical reaction kinetics. Nucleic Acids Res. 2012, 40, D790-D796.

15. Morano, K.A.; Grant, C.M.; Moye-Rowley, W.S. The response to heat shock and oxidative stress in Saccharomyces cerevisiae. Genetics 2012, 190, 1157-1195.

16. Le Novère, N.; Hucka, M.; Hoops, S.; Keating, S.; Sahle, S.; Wilkinson, D. Systems biology markup language (SBML) level 2: Structures and facilities for model definitions. Nat. Preced. 2008, 1-38.

17. Holland, J.H. Outline for a logical theory of adaptive systems. J. ACM 1962, 9, 297-314.

18. Smallbone, K.; Messiha, H.L.; Carroll, K.M.; Winder, C.L.; Malys, N.; Dunn, W.B.; Murabito, E.; Swainston, N.; Dada, J.O.; Khan, F.; et al. A model of yeast glycolysis based on a consistent kinetic characterisation of all its enzymes. FEBS Lett. 2013, 587, 2832-2841.

19. Puig-Castellví, F.; Alfonso, I.; Piña, B.; Tauler, R. A quantitative 1H NMR approach for evaluating the metabolic response of Saccharomyces cerevisiae to mild heat stress. Metabolomics 2015, 11, 1612-1625.

20. Blank, L.M.; Kuepfer, L.; Sauer, U. Large-scale 13 C-flux analysis reveals mechanistic principles of metabolic network robustness to null mutations in yeast. Genome Biol. 2005, 6, 49.

21. Jarnuczak, A.F.; Albornoz, M.G.; Eyers, C.E.; Grant, C.M.; Hubbard, S.J. A quantitative and temporal map of proteostasis during heat shock in Saccharomyces cerevisiae. Mol. Omi. 2018, 14, 37-52.

22. Dräger, A.; Rodriguez, N.; Dumousseau, M.; Dörr, A.; Wrzodek, C.; Le Novère, N.; Zell, A.; Hucka, M. JSBML: A flexible java library for working with SBML. Bioinformatics 2011, 27, 2167-2168.

23. Machné, R.; Finney, A.; Müller, S.; Lu, J.; Widder, S.; Flamm, C. The SBML ODE solver library: A native API for symbolic and fast numerical analysis of reaction networks. Bioinformatics 2006, 22, 1406-1407.

24. Yap, C.F. Trehalose SBML Model. Available online: https://github.com/chuanfuyap/Convenient-Modeller/blob/master/example_models/Trehalose Cycle Model.xml (accessed on 11 September 2020).

25. Voit, E.O. Biochemical and genomic regulation of the trehalose cycle in yeast: Review of observations and canonical model analysis. J. Theor. Biol. 2003, 223, 55-78.

26. Fonseca, L.L.; Sánchez, C.; Santos, H.; Voit, E.O. Complex coordination of multi-scale cellular responses to environmental stress. Mol. Biosyst. 2011, 7, 731-741.

27. Rizk, M.L.; Liao, J.C. Ensemble modeling for aromatic production in Escherichia coli. PLoS ONE 2009, 4.

28. Contador, C.A.; Rizk, M.L.; Asenjo, J.A.; Liao, J.C. Ensemble modeling for strain development of l-lysine-producing Escherichia coli. Metab. Eng. 2009, 11, 221-233.

29. Kacser, H.; Burns, J.A. The control of flux. Symp. Soc. Exp. Biol. 1973, 27, 65-104.

30. Heinrich, R.; Rapoport, T.A. A linear steady-state treatment of enzymatic chains. General properties, control and effector strength. Eur. J. Biochem. 1974, 42, 89-95.

31. Ru-Dong, L.; Li, Y.-Y.; Lu, L.-Y.; Ren, C.; Li, Y.-X.; Liu, L. An improved kinetic model for the acetone-butanol-ethanol pathway of Clostr...: EBSCOhost. BMC Syst. Biol. 2011, 5, 1-11.

32. Li, Y.; Wang, J.; Zhou, J.; Wu, H.; Jin, Q. Kinetic model reconstruction in silico of NSD metabolism network. In Proceedings of the 2015 Chinese Automation Congress CAC, Wuhan, China, 27-29 November 2015; pp. 988-992.

33. Chen, N.; Koumpouras, G.C.; Polizzi, K.M.; Kontoravdi, C. Genome-based kinetic modeling of cytosolic glucose metabolism in industrially relevant cell lines: Saccharomyces cerevisiae and Chinese hamster ovary cells. Bioprocess Biosyst. Eng. 2012, 35, 1023-1033.

34. Gasch, A.P.; Spellman, P.T.; Kao, C.M.; Carmel-Harel, O.; Eisen, M.B.; Storz, G.; Botstein, D.; Brown, P.O. Genomic expression programs in the response of yeast cells to environmental changes. Mol. Biol. Cell 2000, 11, 4241-4257.

35. Chu, Y.; Hahn, J. Parameter set selection for estimation of nonlinear dynamic systems. AIChE J. 2007, 53, 2858-2870.

36. Lee, D.; Ding, Y.; Jayaraman, A.; Kwon, J.S. Mathematical modeling and parameter estimation of intracellular signaling pathway: Application to LPS-induced NFκB activation and TNFα production in macrophages. Processes 2018, 6, 1-19.

37. Kiparissides, A.; Koutinas, M.; Kontoravdi, C.; Mantalaris, A.; Pistikopoulos, E.N. "Closing the loop" in biological systems modeling-from the in silico to the in vitro. Automatica 2011, 47, 1147-1155.

38. Brown, K.S.; Sethna, J.P. Statistical mechanical approaches to models with many poorly known parameters. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 2003, 68, 9.

39. Brown, K.S.; Hill, C.C.; Calero, G.A.; Myers, C.R.; Lee, K.H.; Sethna, J.P.; Cerione, R.A. The statistical mechanics of complex signaling networks: Nerve growth factor signaling. Phys. Biol. 2004, 1, 184-195.

40. Fonseca, L.L.; Chen, P.-W.; Voit, E.O. Canonical modeling of the multi-scale regulation of the heat stress response in yeast. Metabolites 2012, 2, 221-241.

41. Magalhães, R.S.S.; Popova, B.; Braus, G.H.; Outeiro, T.F.; Eleutherio, E.C.A. The trehalose protective mechanism during thermal stress in Saccharomyces cerevisiae: The roles of Ath1 and Agt1. FEMS Yeast Res. 2018, 18.

42. Mahmud, S.A.; Hirasawa, T.; Shimizu, H. Differential importance of trehalose accumulation in Saccharomyces cerevisiae in response to various environmental stresses. J. Biosci. Bioeng. 2010, 109, 262-266.

43. Perez-Riverol, Y.; Csordas, A.; Bai, J.; Bernal-Llinares, M.; Hewapathirana, S.; Kundu, D.J.; Inuganti, A.; Griss, J.; Mayer, G.; Eisenacher, M.; et al. The PRIDE database and related tools and resources in 2019: Improving support for quantification data. Nucleic Acids Res. 2019, 47, D442-D450.