From steady-state to synchronized yeast glycolytic oscillations II: model validation

Kinetic modelling of glycolytic oscillations

Glycolytic oscillations have been studied for well over 60 years, but aspects of their function, and mechanisms of regulation and synchronisation remain unclear. Glycolysis is amenable to mechanistic mathematical modelling, as its components have been well characterised, and the system can be studied at many organisational levels: in vitro reconstituted enzymes, cell free extracts, individual cells, and cell populations. In recent years, the emergence of individual cell analysis has opened new ways of studying this intriguing system.

Exploring the evolution of biochemical models at the network level

The evolution of biochemical models is difficult to track. At present, it is not possible to inspect the differences between model versions at the network level. Biochemical models are often constructed in a distributed, non-linear process: collaborators create model versions on different branches from novel information, model extensions, during curation and adaption. To discuss and align the versions, it is helpful to abstract the changes to the network level. The differences between two model versions can be detected by the software tool BiVeS. However, it cannot show the structural changes resulting from the differences. Here, we present a method to visualise the differences between model versions effectively. We developed a JSON schema to communicate the differences at the network level and extended BiVeS accordingly. Additionally, we developed DiVil, a web-based tool to represent the model and the differences as a standardised network using D3. It combines an automatic layout with an interactive user interface to improve the visualisation and to inspect the model. The network can be exported in standardised formats as images or markup language. Our method communicates the structural differences between model versions. It facilitates the discussion of changes and thus supports the collaborative and non-linear nature of model development.

Availability and implementation: DiVil prototype: https://divil.bio.informatik.uni-rostock.de, Code on GitHub: https://github.com/Gebbi8/DiVil, licensed under Apache License 2.0.

Contact:url="tom.gebhardt@uni-rostock.de.

Differential equation-based minimal model describing metabolic oscillations in Bacillus subtilis biofilms

Biofilms offer an excellent example of ecological interaction among bacteria. Temporal and spatial oscillations in biofilms are an emerging topic. In this paper, we describe the metabolic oscillations in Bacillus subtilis biofilms by applying the smallest theoretical chemical reaction system showing Hopf bifurcation proposed by Wilhelm and Heinrich in 1995. The system involves three differential equations and a single bilinear term. We specifically select parameters that are suitable for the biological scenario of biofilm oscillations. We perform computer simulations and a detailed analysis of the system including bifurcation analysis and quasi-steady-state approximation. We also discuss the feedback structure of the system and the correspondence of the simulations to biological observations. Our theoretical work suggests potential scenarios about the oscillatory behaviour of biofilms and also serves as an application of a previously described chemical oscillator to a biological system.

Effect of macromolecular crowding on the kinetics of glycolytic enzymes and the behaviour of glycolysis in yeast

Water is involved in all aspects of biological activity, both as a solvent and as a reactant. It is hypothesized that intracellular water is in a highly structured state due to the high concentrations of macromolecules in the cell and that this may change the activity of intracellular enzymes due to altered binding affinities and allosteric regulations. Here we first investigate the kinetics of two glycolytic enzymes in artificially crowded aqueous solutions and show that crowding does indeed change their kinetics. Based on our kinetic measurements we propose a new model of oscillating glycolysis that instead of Michaelis-Menten or Monod-Wyman-Changeux kinetics uses the Yang-Ling adsorption isotherm introduced by G. Ling in the frame of the Association-Induction (AI) hypothesis. Using this model, we can reproduce previous experimental observations of the coupling of glycolytic oscillations and intracellular water dynamics, e.g., (i) during the metabolic oscillations, the latter variable oscillates in phase with ATP activity, and (ii) the emergence of glycolytic oscillations largely depends on the extent of intracellular water dipolar relaxation in cells in the resting state. Our results support the view that the extent of intracellular water dipolar relaxation is regulated by the ability of cytoplasmic proteins to polarize intracellular water with the assistance of ATP, as suggested in the AI hypothesis. This hypothesis may be relevant to the interpretation of many other biological oscillators, including cell signalling processes.

Modeling studies of heterogeneities in glycolytic oscillations in HeLa cervical cancer cells

Previous experiments demonstrated that a population of HeLa cells starved of glucose or both glucose and serum exhibited a strong heterogeneity in the glycolytic oscillations in terms of the number of oscillatory cells, periods of oscillations, and duration of oscillations. Here, we report numerical simulations of this heterogeneous oscillatory behavior in HeLa cells by using a newly developed mathematical model. It is simple enough that we can apply a mathematical analysis, but capture the core of the glycolytic pathway and the activity of the glucose transporter (GLUT). Lognormal distributions of the values of the four rate constants in the model were obtained from the experimental distributions in the periods of oscillations. Thus, the heterogeneity in the periods of oscillations can be attributed to the difference in the rate constants of the enzymatic reactions. The activity of GLUT is found to determine whether the HeLa cells were oscillatory or non-oscillatory under the same experimental conditions. Simulation with the log-normal distribution of the maximum uptake velocity of glucose and the four randomized rate constants based on the log-normal distributions successfully reproduced the time-dependent number of oscillatory cells (oscillatory ratios) under the two starving conditions. The difference in the initial values of the metabolites has little effect on the simulated results.

Phosphofructokinase controls the acetaldehyde-induced phase shift in isolated yeast glycolytic oscillators

The response of oscillatory systems to external perturbations is crucial for emergent properties such as synchronisation and phase locking and can be quantified in a phase response curve (PRC). In individual, oscillating yeast cells, we characterised experimentally the phase response of glycolytic oscillations for external acetaldehyde pulses and followed the transduction of the perturbation through the system. Subsequently, we analysed the control of the relevant system components in a detailed mechanistic model. The observed responses are interpreted in terms of the functional coupling and regulation in the reaction network. We find that our model quantitatively predicts the phase-dependent phase shift observed in the experimental data. The phase shift is in agreement with an adaptation leading to synchronisation with an external signal. Our model analysis establishes that phosphofructokinase plays a key role in the phase shift dynamics as shown in the PRC and adaptation time to external perturbations. Specific mechanism-based interventions, made possible through such analyses of detailed models, can improve upon standard trial and error methods, e.g. melatonin supplementation to overcome jet-lag, which are error-prone, specifically, since the effects are phase dependent and dose dependent. The models by Gustavsson and Goldbeter discussed in the text can be obtained from the JWS Online simulation database: (https://jjj.bio.vu.nl/models/gustavsson5 and https://jjj.bio.vu.nl/models/goldbeter1).

Effect of the periodic properties of toxic stress on the oscillatory behaviour of glycolysis in yeast-evidence of a toxic effect frequency

Starving and nondividing yeast cells induce changes in the electron donor nicotinamide adenine dinucleotide (NADH) levels in a cyclic and wave-like manner for over 90min. Yeast suspensions were used to examine the toxic effects of contaminants on the cyclic behaviour of metabolite changes during anaerobic glycolysis. The cyclic behaviour NADH levels in yeast cell suspensions starved for 2 to 5h was studied after the addition of 10mM glucose for 5min followed by 10mM KCN to block aerobic glycolysis. The effects of three toxic elements (CuSO₄, silver nanoparticles-nAg, and GdCl₃), known for their potential to alter glycolsysis, on NADH levels over time were examined during the 3-h starvation step. The data were analyzed using spectral analysis (Fourier transformation) to characterize the cyclic behaviour of NADH levels during anaerobic glycolysis. Increasing the starvation time by 3h increased the amplitude of changes in NADH levels with characteristic periods of 3 to 8min. Longer starvation times decreased the amplitude of oscillations during these periods, with the appearance of NADH changes at higher frequencies. Moreover, the amplitude changes in NADH were proportional to the starvation time. Exposure to the above chemicals during the 3-h starvation time led to the formation of higher frequencies with concentration-dependent amplitude changes. In comparison with nAg and Gd³⁺, Cu²⁺ was the most toxic (decreased viability the most) and produce changes at higher frequencies as well. It is noteworthy that each element produced a characteristic change in the frequency profiles, which suggests different mechanisms of action in which the severity of toxicity shifted NADH changes at higher frequencies. In conclusion, the appearance of synchronized oscillations in dense yeast populations following synchronization stress could be induced by starvation and exposure to chemicals. However, synchronicity could be abolished when cells desynchronize as a result of loss of cell viability, which contributes to heterogeneity in yeast populations, translating into NADH changes at higher frequencies. This is the first report on the influence of environmental contaminants on the cyclic or wave-like behaviour of biochemical changes in cells.

Validation and selection of ODE based systems biology models: how to arrive at more reliable decisions

Background

Most ordinary differential equation (ODE) based modeling studies in systems biology involve a hold-out validation step for model validation. In this framework a pre-determined part of the data is used as validation data and, therefore it is not used for estimating the parameters of the model. The model is assumed to be validated if the model predictions on the validation dataset show good agreement with the data. Model selection between alternative model structures can also be performed in the same setting, based on the predictive power of the model structures on the validation dataset. However, drawbacks associated with this approach are usually under-estimated.

Results

We have carried out simulations by using a recently published High Osmolarity Glycerol (HOG) pathway from S.cerevisiae to demonstrate these drawbacks. We have shown that it is very important how the data is partitioned and which part of the data is used for validation purposes. The hold-out validation strategy leads to biased conclusions, since it can lead to different validation and selection decisions when different partitioning schemes are used. Furthermore, finding sensible partitioning schemes that would lead to reliable decisions are heavily dependent on the biology and unknown model parameters which turns the problem into a paradox. This brings the need for alternative validation approaches that offer flexible partitioning of the data. For this purpose, we have introduced a stratified random cross-validation (SRCV) approach that successfully overcomes these limitations.

Conclusions

SRCV leads to more stable decisions for both validation and selection which are not biased by underlying biological phenomena. Furthermore, it is less dependent on the specific noise realization in the data. Therefore, it proves to be a promising alternative to the standard hold-out validation strategy.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-015-0180-0) contains supplementary material, which is available to authorized users.

Multiplicity of steady states in glycolysis and shift of metabolic state in cultured mammalian cells

Cultured mammalian cells exhibit elevated glycolysis flux and high lactate production. In the industrial bioprocesses for biotherapeutic protein production, glucose is supplemented to the culture medium to sustain continued cell growth resulting in the accumulation of lactate to high levels. In such fed-batch cultures, sometimes a metabolic shift from a state of high glycolysis flux and high lactate production to a state of low glycolysis flux and low lactate production or even lactate consumption is observed. While in other cases with very similar culture conditions, the same cell line and medium, cells continue to produce lactate. A metabolic shift to lactate consumption has been correlated to the productivity of the process. Cultures that exhibited the metabolic shift to lactate consumption had higher titers than those which didn't. However, the cues that trigger the metabolic shift to lactate consumption state (or low lactate production state) are yet to be identified. Metabolic control of cells is tightly linked to growth control through signaling pathways such as the AKT pathway. We have previously shown that the glycolysis of proliferating cells can exhibit bistability with well-segregated high flux and low flux states. Low lactate production (or lactate consumption) is possible only at a low glycolysis flux state. In this study, we use mathematical modeling to demonstrate that lactate inhibition together with AKT regulation on glycolysis enzymes can profoundly influence the bistable behavior, resulting in a complex steady-state topology. The transition from the high flux state to the low flux state can only occur in certain regions of the steady state topology, and therefore the metabolic fate of the cells depends on their metabolic trajectory encountering the region that allows such a metabolic state switch. Insights from such switch behavior present us with new means to control the metabolism of mammalian cells in fed-batch cultures.

Understanding bistability in yeast glycolysis using general properties of metabolic pathways

Glycolysis is the central pathway in energy metabolism in the majority of organisms. In a recent paper, van Heerden et al. showed experimentally and computationally that glycolysis can exist in two states, a global steady state and a so-called imbalanced state. In the imbalanced state, intermediary metabolites accumulate at low levels of ATP and inorganic phosphate. It was shown that Baker's yeast uses a peculiar regulatory mechanism--via trehalose metabolism--to ensure that most yeast cells reach the steady state and not the imbalanced state.

RESULTS

Here we explore the apparent bistable behaviour in a core model of glycolysis that is based on a well-established detailed model, and study in great detail the bifurcation behaviour of solutions, without using any numerical information on parameter values.

CONCLUSION

We uncover a rich suite of solutions, including so-called imbalanced states, bistability, and oscillatory behaviour. The techniques employed are generic, directly suitable for a wide class of biochemical pathways, and could lead to better analytical treatments of more detailed models.

How informative is your kinetic model?: using resampling methods for model invalidation

Background

Kinetic models can present mechanistic descriptions of molecular processes within a cell. They can be used to predict the dynamics of metabolite production, signal transduction or transcription of genes. Although there has been tremendous effort in constructing kinetic models for different biological systems, not much effort has been put into their validation. In this study, we introduce the concept of resampling methods for the analysis of kinetic models and present a statistical model invalidation approach.

Results

We based our invalidation approach on the evaluation of a kinetic model’s predictive power through cross validation and forecast analysis. As a reference point for this evaluation, we used the predictive power of an unsupervised data analysis method which does not make use of any biochemical knowledge, namely Smooth Principal Components Analysis (SPCA) on the same test sets. Through a simulations study, we showed that too simple mechanistic descriptions can be invalidated by using our SPCA-based comparative approach until high amount of noise exists in the experimental data. We also applied our approach on an eicosanoid production model developed for human and concluded that the model could not be invalidated using the available data despite its simplicity in the formulation of the reaction kinetics. Furthermore, we analysed the high osmolarity glycerol (HOG) pathway in yeast to question the validity of an existing model as another realistic demonstration of our method.

Conclusions

With this study, we have successfully presented the potential of two resampling methods, cross validation and forecast analysis in the analysis of kinetic models’ validity. Our approach is easy to grasp and to implement, applicable to any ordinary differential equation (ODE) type biological model and does not suffer from any computational difficulties which seems to be a common problem for approaches that have been proposed for similar purposes. Matlab files needed for invalidation using SPCA cross validation and our toy model in SBML format are provided at http://www.bdagroup.nl/content/Downloads/software/software.php.

Allosteric regulation of phosphofructokinase controls the emergence of glycolytic oscillations in isolated yeast cells

Oscillations are widely distributed in nature and synchronization of oscillators has been described at the cellular level (e.g. heart cells) and at the population level (e.g. fireflies). Yeast glycolysis is the best known oscillatory system, although it has been studied almost exclusively at the population level (i.e. limited to observations of average behaviour in synchronized cultures). We studied individual yeast cells that were positioned with optical tweezers in a microfluidic chamber to determine the precise conditions for autonomous glycolytic oscillations. Hopf bifurcation points were determined experimentally in individual cells as a function of glucose and cyanide concentrations. The experiments were analyzed in a detailed mathematical model and could be interpreted in terms of an oscillatory manifold in a three-dimensional state-space; crossing the boundaries of the manifold coincides with the onset of oscillations and positioning along the longitudinal axis of the volume sets the period. The oscillatory manifold could be approximated by allosteric control values of phosphofructokinase for ATP and AMP.

DATABASE

The mathematical models described here have been submitted to the JWS Online Cellular Systems Modelling Database and can be accessed at http://jjj.mib.ac.uk/webMathematica/UItester.jsp?modelName=gustavsson5. [Database section added 14 May 2014 after original online publication].

Heterogeneity of glycolytic oscillatory behaviour in individual yeast cells

There are many examples of oscillations in biological systems and one of the most investigated is glycolytic oscillations in yeast. These oscillations have been studied since the 1950s in dense, synchronized populations and in cell-free extracts, but it has for long been unknown whether a high cell density is a requirement for oscillations to be induced, or if individual cells can oscillate also in isolation without synchronization. Here we present an experimental method and a detailed kinetic model for studying glycolytic oscillations in individual, isolated yeast cells and compare them to previously reported studies of single-cell oscillations. The importance of single-cell studies of this phenomenon and relevant future research questions are also discussed.

An experimental study of the regulation of glycolytic oscillations in yeast

We have studied oscillating glycolysis in the strain BY4743 and isogenic strains with deletions of genes encoding enzymes in glycolysis, mitochondrial electron transport and ATP synthesis. We found that deletion of the gene encoding the hexokinase 1 isoform does not affect the oscillations while deletion of the gene encoding the hexokinase 2 isoform results in oscillations with smaller amplitude. The latter is associated with an almost 50% decrease in hexokinase activity. Deletions in the genes encoding the α- and β-subunits of phosphofructokinase abolish the oscillations entirely. This loss in oscillatory activity is associated with a fourfold decrease in phosphofructokinase activity. Deletions of genes encoding subunits of the F1F0 ATPase also inhibit the oscillations in accordance with earlier studies using for example inhibitors. Finally, we identified an apparently new control point involving the mitochondrial cytochrome c oxidase. The latter is difficult to explain as oscillatory activity entails 100% inhibition of this enzyme. The mitochondria of this strain seem to have normal F1F0 ATPase activity. Overall these results support earlier experimental and model studies suggesting that in addition to processes within glycolysis also processes outside this pathway contribute to the control of the oscillatory behaviour.

Oscillatory enzyme reactions and Michaelis-Menten kinetics

Oscillations occur in a number of enzymatic systems as a result of feedback regulation. How Michaelis-Menten kinetics influences oscillatory behavior in enzyme systems is investigated in models for oscillations in the activity of phosphofructokinase (PFK) in glycolysis and of cyclin-dependent kinases in the cell cycle. The model for the PFK reaction is based on a product-activated allosteric enzyme reaction coupled to enzymatic degradation of the reaction product. The Michaelian nature of the product decay term markedly influences the period, amplitude and waveform of the oscillations. Likewise, a model for oscillations of Cdc2 kinase in embryonic cell cycles based on Michaelis-Menten phosphorylation-dephosphorylation kinetics shows that the occurrence and amplitude of the oscillations strongly depend on the ultrasensitivity of the enzymatic cascade that controls the activity of the cyclin-dependent kinase.

Programming strategy for efficient modeling of dynamics in a population of heterogeneous cells

MOTIVATION

Heterogeneity is a ubiquitous property of biological systems. Even in a genetically identical population of a single cell type, cell-to-cell differences are observed. Although the functional behavior of a given population is generally robust, the consequences of heterogeneity are fairly unpredictable. In heterogeneous populations, synchronization of events becomes a cardinal problem-particularly for phase coherence in oscillating systems.

RESULTS

The present article presents a novel strategy for construction of large-scale simulation programs of heterogeneous biological entities. The strategy is designed to be tractable, to handle heterogeneity and to handle computational cost issues simultaneously, primarily by writing a generator of the 'model to be simulated'. We apply the strategy to model glycolytic oscillations among thousands of yeast cells coupled through the extracellular medium. The usefulness is illustrated through (i) benchmarking, showing an almost linear relationship between model size and run time, and (ii) analysis of the resulting simulations, showing that contrary to the experimental situation, synchronous oscillations are surprisingly hard to achieve, underpinning the need for tools to study heterogeneity. Thus, we present an efficient strategy to model the biological heterogeneity, neglected by ordinary mean-field models. This tool is well posed to facilitate the elucidation of the physiologically vital problem of synchronization.

AVAILABILITY

The complete python code is available as Supplementary Information.

CONTACT

bjornhald@gmail.com or pgs@kiku.dk

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.