A pathway model of glucose-stimulated insulin secretion in the pancreatic β-cell

The pancreas plays a critical role in maintaining glucose homeostasis through the secretion of hormones from the islets of Langerhans. Glucose-stimulated insulin secretion (GSIS) by the pancreatic β-cell is the main mechanism for reducing elevated plasma glucose. Here we present a systematic modeling workflow for the development of kinetic pathway models using the Systems Biology Markup Language (SBML). Steps include retrieval of information from databases, curation of experimental and clinical data for model calibration and validation, integration of heterogeneous data including absolute and relative measurements, unit normalization, data normalization, and model annotation. An important factor was the reproducibility and exchangeability of the model, which allowed the use of various existing tools. The workflow was applied to construct a novel data-driven kinetic model of GSIS in the pancreatic β-cell based on experimental and clinical data from 39 studies spanning 50 years of pancreatic, islet, and β-cell research in humans, rats, mice, and cell lines. The model consists of detailed glycolysis and phenomenological equations for insulin secretion coupled to cellular energy state, ATP dynamics and (ATP/ADP ratio). Key findings of our work are that in GSIS there is a glucose-dependent increase in almost all intermediates of glycolysis. This increase in glycolytic metabolites is accompanied by an increase in energy metabolites, especially ATP and NADH. One of the few decreasing metabolites is ADP, which, in combination with the increase in ATP, results in a large increase in ATP/ADP ratios in the β-cell with increasing glucose. Insulin secretion is dependent on ATP/ADP, resulting in glucose-stimulated insulin secretion. The observed glucose-dependent increase in glycolytic intermediates and the resulting change in ATP/ADP ratios and insulin secretion is a robust phenomenon observed across data sets, experimental systems and species. Model predictions of the glucose-dependent response of glycolytic intermediates and biphasic insulin secretion are in good agreement with experimental measurements. Our model predicts that factors affecting ATP consumption, ATP formation, hexokinase, phosphofructokinase, and ATP/ADP-dependent insulin secretion have a major effect on GSIS. In conclusion, we have developed and applied a systematic modeling workflow for pathway models that allowed us to gain insight into key mechanisms in GSIS in the pancreatic β-cell.

1 Introduction

The pancreas plays a vital role in maintaining glucose homeostasis (1) through the secretion of hormones from the islets of Langerhans. The most important hormones are insulin, secreted by the pancreatic β-cells, and glucagon, secreted by the α-cells, both of which play key roles in regulating glucose homeostasis (2).

Glucose-induced insulin secretion (GSIS) is a physiological process by which the pancreas releases insulin in response to an increase in blood glucose levels. When glucose enters the bloodstream after a meal, it is taken up by β-cells in the pancreas through glucose transporters, primarily GLUT2 (3). Once inside the β-cells, glucose is metabolized via glycolysis, which produces energy in the form of ATP.

The coupling of glycolysis with the insulin secretion mechanism in the β-cell is established by the regulatory effects of glycolytic intermediates on the levels of energy metabolites such as ATP and NADH (4, 5). The rise in ATP levels triggers a series of events that lead to the release of insulin. Specifically, the high ATP levels close ATP-sensitive potassium channels (6), which leads to depolarization of the cell membrane and opening of voltage-gated calcium channels. The influx of calcium triggers the exocytosis of insulin-containing vesicles, leading to the release of insulin into the bloodstream (7, 8). The K_ATP/Ca²⁺ independent signaling mechanisms and the other metabolites besides glucose contribute to the amplification of the signaling events that trigger insulin secretion (9).

GSIS by the pancreatic β-cell is the primary mechanism for lowering elevated plasma glucose levels. The amount of insulin released increases with the glucose in the bloodstream. This process is crucial for the regulation of blood glucose levels by promoting the uptake and use of glucose by cells throughout the body, such as muscle, fat tissue, and the liver (10, 11).

Glycolysis is the primary metabolic pathway responsible for GSIS. It involves the uptake of glucose and its conversion to pyruvate, which is critical for ATP synthesis and maintenance of ATP levels. Experimental data from metabolic profiling studies in islet cells support the key role of glycolysis in GSIS (12–14). As glucose levels increase, glycolytic flux and most glycolytic intermediates increase in a dose-dependent manner. Changes in adenine nucleotide levels due to variations in glycolytic flux lead to changes in nucleotide ratios, with increasing glucose levels resulting in a positive correlation between the ATP/ADP ratio and Ca²⁺ response and insulin release. This trend is consistent across several studies (15–17), including isolated islets perfused with glucose, rat and mouse tissue homogenates, and insulin-secreting cell lines. The increase in ATP/ADP ratio ranges from 2 to 7 when glucose levels are increased from 2.8mM to 30mM, indicating similar behavior in different experimental systems studying insulin secretion by the pancreas (18).

Mathematical models have been developed to investigate the metabolic and signaling mechanisms that trigger and amplify insulin secretion. Early models of β-cells focused on examining the relationship between glycolytic oscillations and pulsatile insulin release to understand GSIS (19, 20). Minimal models of GSIS have examined the effect of dosing patterns such as slow and fast ramps of glucose on the phasic nature of insulin secretion (21–24). Merrins et al. analyzed the oscillations in glycolytic intermediates (i.e. fructose-6-phosphate, fructose-2,6-bisphosphate, and fructose-1,6-bisphosphate) and their effect on pulsatile insulin secretion (25), while other models integrated glycolytic flux with mitochondrial ATP production to study the role of reducing equivalents such as pyridine nucleotides in enhancing insulin secretion (26, 27). Jiang et al. further combined previously developed models of glycolysis, citric acid cycle, β-oxidation, pentose phosphate shunt, and respiratory chain and studied the local and global dynamics of the GSIS mechanism in response to parameter perturbations. These models were coupled with the calcium signaling pathway of Fridyland et al. to create an integrated metabolic model (28, 29).

To investigate the synergistic insulinotropic effect of other nutrient sources, Salvucci et al. (17) developed a model by integrating alanine metabolism with glucose metabolism, the citric acid cycle, and the respiratory chain. Gelbach et al. developed a system of 65 reactions integrating glycolysis, glutaminolysis, the pentose phosphate pathway, the citric acid cycle, the polyol pathway, and the electron transport chain to study the kinetics of insulin secretion (30).

However, the majority of these models are based on earlier models that were developed using kinetic data from organisms other than humans or non-pancreatic tissues, such as a glycolysis model that utilized kinetic data from experiments on yeast cell extract, or a glycolysis model based on kinetic data from mammalian muscle (31). Often, the data used to build these models is limited and comes from a single experimental study. In most models specific to β-cells, reaction kinetics are described by simple mass-action rate laws. There exists no detailed kinetic model of the changes in glycolysis during GSIS that can effectively integrate the observed changes in glycolytic and energy intermediates from a wide range of GSIS experiments.

In systems biology and systems medicine, ensuring the reproducibility of computational models and integrating diverse data from multiple sources into these models are critical challenges. Standards for model description, such as the Systems Biology MarkupLanguage (SBML) (32, 33), have been developed to enable the reusability and reproducibility of existing models, but they have yet to be utilized in the field of pancreatic GSIS modeling. Furthermore, there is a need to address how to integrate heterogeneous data from multiple studies conducted in different organisms and experimental systems in the context of GSIS modeling.

This study aims to develop a detailed kinetic model of GSIS and the associated changes in glycolysis in the pancreatic β-cell. The novel contributions of this work include a systematic curation and integration of changes in glycolytic metabolites from multiple experimental studies across different species and experimental systems to construct a new model of GSIS. Based on this unique data set, a detailed kinetic model of glycolysis and GSIS was constructed using a systematic approach with a focus on reproducibility. This approach allowed the establishment of a consensus model of the changes that occur in insulin secretion with varying glucose concentrations. The overall goal was to provide a better understanding of the mechanisms underlying GSIS and to contribute to the development of improved computational models of these processes.

2 Results

Our study introduces a detailed kinetic model of GSIS in the pancreatic β-cell, which has the ability to simulate alterations in glycolytic intermediates and ATP/ADP ratio due to glucose levels and the effect of change in the energy state of the β-cell on biphasic insulin secretion.

2.1 Systematic curation of data set of changes in GSIS

In the course of this study, we compiled a comprehensive data set (Table 1) of GSIS based on experimental and clinical data from 39 studies spanning half a century of research on pancreatic, islet, and β-cell function in humans, rats, mice, and cell lines. Specifically, we systematically curated metabolomics data from studies conducted between 1970 and 2020, comprising information on the concentration of glycolytic intermediates and cofactors in both time-course and steady-state experiments, as well as the corresponding glucose doses. The data set contains 17 metabolites, comprising 359 data points from steady-state experiments and 249 data points from time-course studies. It includes both absolute and relative measurements of metabolite changes, and an overview of the available information for each metabolite and study is presented in Figure 1.

FIGURE 1

Figure 1 Curated data for model development and evaluation. The data description is detailed from the periphery to the center of the Circos plot. 1. Model elements: The outermost layer provides an overview of the metabolites included in the data set. GLC: glucose, G6P: glucose 6-phosphate, F6P: fructose 6-phosphate, FBP: fructose 1,6-bisphosphate, F26BP: fructose 2,6-bisphosphate, DHAP: dihydroxyacetone phosphate, GRAP: glyceraldehyde 3-phosphate, BPG: 1,3-biphosphoglycerate, 3PG: 3-phosphoglycerate, 2PG: 2-phosphoglycerate, PEP: phosphoenolpyruvate, PYR: pyruvate, LAC: lactate, PHOS: phosphate, NAD: nicotinamide adenine dinucleotide, NADH: reduced nicotinamide adenine dinucleotide, NADH total: NADH + NAD; NADH ratio: NADH/NAD; ATP: adenosine triphosphate, ADP: adenosine diphosphate, ATP total: ATP + ADP, ATP ratio: ATP/ADP, IRS: insulin secretion rate. The metabolites were grouped in the following categories: Color code: • glycolytic intermediates, • cofactors, • cofactor ratio or sum, • insulin secretion rate (IRS); 2. Studies: The second layer depicts the islet-cell specific metabolite profiling studies curated from the literature; 3. Animal species: The third layer indicates the animal species or cell line from which the data was curated. Color code: • Rat, • Human, • Mouse, and • Cell line data; 4. time course data: The fourth layer shows a bar graph illustrating the number of data points collected from studies reporting time course data of metabolites. Color code: • relative (or fold), • concentration, • ratio, • rate measurements; 5. Steady-state data: The fifth layer indicates the number of data points collected from studies reporting steady-state/dose-response data of metabolites. Color code: • relative (or fold), • concentration, • ratio, • rate measurements; 6. Data used for parameter estimation: The innermost layer indicates the subset of data used for parameter fitting.

This data set represents the first open and FAIR (findable, accessible, interoperable, and reusable) large-scale collection of data on changes in glycolysis and insulin secretion in the pancreatic β-cell during GSIS. We used the absolute and relative measurements of glycolysis metabolites and insulin secretion rates in this data set for model calibration and evaluation. To the best of our knowledge, this dataset is the first open-access resource on pancreatic β-cell glycolysis that is easily accessible to the scientific community for further use.

The data set is available under a CC-BY4.0 license from https://github.com/matthiaskoenig/pancreas-model.

2.2 Reproducible modeling workflow

In this study, we describe a comprehensive modeling workflow for building small kinetic pathway models (Figure 2) using SBML (32, 33).

FIGURE 2

Figure 2 Kinetic model development workflow. (A) Initial stoichiometric model in SBML. Glycolytic reactions were collected from VMH database and existing models of glycolysis. (B) Metadata integration. VMH and BiGG database field identifiers were used to retrieve additional metadata such as HMDB, BioCyc, MetaNetX, ChEBI, and SEED database field identifiers. (C) Synonym mapping. The synonyms associated with each metabolite were queried using compound identifier mapping services. (D) Kinetic parameters. EC number and KEGG reaction identifiers were used to query half-saturation/Michaelis-Menten K_M, inhibition K_I, activation K_A, and equilibrium K_eq constants (synonym mapping was applied for all compounds). (E) Model parameters. The parameter values retrieved from different databases were merged and median values were assigned to the model parameters. (F) Data curation. A systematic literature search was performed and metabolite concentrations from islet cell studies were curated. (G) Unit normalization. Absolute and relative quantification of metabolite concentrations reported in heterogeneous units were converted to mM. (H) Data normalization. Systematic bias observed in the unit-normalized data was removed by performing least-squares minimization to minimize the distance between the mean curve of the unit-normalized data curves and the experimental curves of the unit-normalized data. (I) Model inputs. Values of kinetic parameters, initial concentrations, volumes, equations, and annotations have been assigned to the model entities. (J) Model calibration. Time course and steady-state data were used for parameter estimation. (K) Kinetic SBML model. The final kinetic SBML model was generated. (L) Model prediction. Glycolytic intermediates and insulin response were predicted as a function of varying glucose concentrations. Created with BioRender.com.

In our model-building workflow, we followed several steps to construct a novel kinetic SBML model of glycolysis in the pancreatic β-cell. A) First, we built an SBML model based on the stoichiometry of glycolytic reactions and intermediates from existing models and pathway databases. B) We then annotated metabolites and reactions with metadata information which was extended by querying VMH and the BiGG database, resulting in mappings to additional resources such as HMDB, BioCyc, MetaNetX, ChEBI, and SEED. C) and D) We collected and retrieved kinetic parameters such as K_M, K_I, K_A, and K_eq constants from databases and integrated them with synonyms associated with each queried metabolite using compound identifier mapping services. E) We integrated the resulting parameters and assigned median values to the model parameters. F) Next, we curated data from studies reporting metabolite concentrations and changes, and insulin secretion in pancreatic, islet, and β-cell lines through a literature search. G) Unit normalization was then performed to convert reported metabolite concentrations and insulin secretion to mmol/l (mM) and nmol/min/ml (β-cell volume), respectively. H) Data normalization was performed to remove systematic differences between data reported in different studies and experimental systems. I) Next, values for kinetic parameters, initial concentrations, volumes, rate equations, and annotations were integrated into the stoichiometric model. J) We calibrated the model by parameter optimization using time-course and steady-state data and K) generated the final SBML kinetic model using all the information. L) Finally, we performed model predictions of glycolytic intermediates and insulin response as a function of varying glucose concentrations. Steps were performed iteratively to fill gaps and extend the data set and model.

2.3 Computational model

Using the established data set, we utilized the aforementioned workflow to develop a novel data-driven kinetic model of GSIS in the pancreatic β-cell. The model is comprised of detailed glycolysis and equations for insulin secretion which are coupled to the cellular energy state (ATP/ADP ratio) and change in ATP (dATP/dt). The metabolites and reactions incorporated into the kinetic model are depicted in Figure 3, and their biochemical interactions are represented through a system of ordinary differential equations. The model consists of 21 enzyme-catalyzed reactions, 25 metabolites, and 91 parameters, and also includes an empirical model that connects the energy state of the β-cell to insulin secretion.

FIGURE 3

Figure 3 Computational model of glucose-stimulated insulin secretion (GSIS) in the pancreatic β-cell. The model consists of glycolysis and insulin secretion coupled to the energy state (ATP/ADP ratio). The GLUT transporter facilitates the uptake of glucose from the plasma into the cell. Glucose undergoes phosphorylation and the subsequent reactions lead to the production of pyruvate. Pyruvate can either be converted to lactate and exported into blood or transported to the mitochondria where it serves as a fuel source for the production of tricarboxylic acid cycle (TCA) intermediates (the TCA cycle has not been modeled). Depending on the external glucose concentrations, glycolysis intermediates and energy metabolites such as ATP, ADP, NAD, and NADH change. An increase in the ATP/ADP ratio as a result of changes in glucose triggers the cascade of signaling mechanisms that promote insulin secretion by the pancreatic β-cell. Phosphate, water, and hydrogen ions have been omitted from the diagram for clarity (but are included in the model for mass and charge balance). The network diagram was created using CySBML (66). Created with BioRender.com

When glucose levels are high, GLUT transporter allows glucose to enter the cell, and glucokinase converts glucose to glucose-6-phosphate. The upper glycolysis produces fructose-6-phosphate, fructose-1,6-phosphate, and triose phosphates like dihydroxyacetone phosphate and glyceraldehyde phosphate. Lower glycolysis then leads to the creation of 3-phosphoglycerate, 2-phosphoglycerate, phosphoenolpyruvate, and pyruvate. Pyruvate can be transformed into lactate or transported to the mitochondria. For each glucose molecule, two ATP molecules are produced. Changes in ATP/ADP ratio and ATP trigger insulin secretion.

The SBML model is available under a CC-BY4.0 license from https://github.com/matthiaskoenig/pancreas-model.

2.4 Normalization of data

The aim of this study was to investigate variations in glycolysis, glycolytic intermediates, energy metabolites, and insulin secretion during GSIS using the established model. In order to integrate heterogeneous experimental data for each metabolite and insulin secretion rate, we conducted a two-step normalization process to standardize time course and dose-response measurements. The normalization process involved unit normalization (as discussed in Sec. 4.7) and data normalization (as discussed in Sec. 4.8) to normalize the diverse data and eliminate systematic deviations for individual studies. We present the case of glucose 6-phosphate as an example of the normalization process (see Figure 4). The experimental curves were converted to relative (fold) and unit-normalized absolute measurements (Figures 4A, B). To combine the fold data and absolute data, we multiplied the fold values by the basal concentration to obtain absolute values (Figure 4C). If the basal metabolite concentration was not reported, we used the mean curve of the absolute data at the pre-incubation glucose dose of the experiment to determine the basal value. For metabolites consisting of only relative measurements, we used the half-saturation K_m value of the metabolite as an estimate for the basal concentration. Using this strategy, we converted all fold-changes and time courses to absolute data with standardized units, which was then combined with the existing absolute data. However, the variability of the combined measurements was high, and large systematic differences between studies could be observed. We determined scaling factors for every study to minimize the difference between all studies based on least-squares minimization (as discussed in Sec. 4.8.1). The resulting normalized data (Figure 4D) was then used for model calibration. We applied this procedure to all metabolites in the model as well as the insulin secretion rate, reducing the variability in the data substantially.

FIGURE 4

Figure 4 Normalization of steady-state and time course data for glucose 6-phosphate (G6P). (A) Relative data. Experimental curves from β-cell studies reporting relative levels of G6P, expressed as fold to baseline value; (B) Absolute data. Experimental curves from β-cell studies reporting absolute concentrations of G6P, the plot displays the unit-normalized absolute data. (C) Combined data. The relative (fold) measurements were converted to absolute units and combined with the unit-normalized absolute data. (D) Normalized data. Systematic biases between different studies of the combined data were removed by data normalization. Data normalization was performed by minimizing the offset (sum of squared residuals) between the mean curve and the experimental curves. The mean curve was computed as the weighted average of the experimental curves and spline curve is the piecewise-polynomial interpolation of the data points in the mean curve. For steady-state data, the legend indicates studies associated with the experimental curves. For time course data, the legend indicates the pre-incubation glucose dose (☆), incubation glucose dose (□), experimental study, and the value of scale transformation parameter f^α (▭) of experiment α. (top panel) and (bottom panel) show the data of dose-response and time course experiments, respectively. Data from (12, 14, 18, 34, 35, 37, 39, 45, 46, 48, 56–58, 63). For more details, please refer to Sec. 2.1.

2.5 Changes in glycolytic metabolites and insulin secretion in GSIS

Our work has uncovered several key findings related to GSIS. First, we found that almost all glycolytic intermediates increase in a glucose-dependent manner across a wide range of glucose concentrations, as illustrated in Figures 5A–C. This increase in glycolytic intermediates is accompanied by a corresponding increase in energy metabolites, especially ATP and NADH. However, one notable exception is ADP, which decreases with increasing glucose levels. As a result, there is a significant increase in ATP/ADP ratios in β-cells with increasing glucose, a key factor in insulin secretion. This phenomenon is robust across different data sets, experimental systems, and species. An important observation is that not only ATP and NADH increase with increasing glucose, but also the total ATP (ATP + ADP) and total NADH (NAD + NADH).

FIGURE 5

Figure 5 (A) Effect of variations in blood glucose on glycolytic intermediates. (left column) Time course scan. The effect of variation in blood glucose dose on the transient concentration of metabolites; ▪ indicates the mean value of experimental steadystate measurements. (middle column) Dose-response scan. Glucose scan was performed for the calculation of steady-state concentration of metabolites in the model. The steady-state concentrations predicted by the model at various glucose doses were compared with the normalized values of experimental measurements; (right column) Time course. Time course values of glycolytic intermediates and cofactors from multiple experimental studies and the model simulations carried out at the corresponding pre-incubation and incubation doses of glucose. (⭐) in the legend indicates the pre-incubation glucose dose. GLC, glucose; G6P, glucose 6-phosphate; F6P, fructose 6-phosphate; FBP, fructose 1,6-bisphosphate; F26BP, fructose 2,6-bisphosphate; DHAP, dihydroxyacetone phosphate; GRAP, glyceraldehyde 3-phosphate. Data from (12, 14, 18, 34, 35, 37, 39, 45, 46, 48, 49, 56–58, 63). (B) Effect of variations in blood glucose on glycolytic intermediates. The plot is analogous to (A). BPG, 1,3-biphosphoglycerate; 2PG, 2-phosphoglycerate; 3PG, 3-phosphoglycerate; PEP, phosphoenolpyruvate; PYR, pyruvate; LAC, lactate; PHOS, phosphate. Data from (9, 13, 14, 18, 44, 47, 51, 63, 65; 40, 45, 46). (C) Effect of variations in blood glucose on glycolytic cofactors. The plot is analogous to (A). NAD, nicotinamide adenine dinucleotide; NADH, nicotinamide adenine dinucleotide reduced. ATP, adenosine triphosphate; ADP, adenosine diphosphate Data from 9, 12, 14–16, 34, 35, 38, 41–43, 45, 46, 51, 54, 55, 62).

Our model was able to predict the glucose-dependent response of glycolytic intermediates and insulin secretion with good agreement to most experimental measurements, as summarized in Table 1. We observed a dose-dependent increase in glycolytic intermediates when glucose concentrations were increased from 1 mM to 35 mM. The model predicts that steady states of glycolytic metabolites under constant glucose are reached after approximately 20 minutes, which is in good agreement with the data.

TABLE 1

Table 1 Overview of studies reporting concentrations of metabolites used for model calibration.

Figure 6A illustrates the relationship between glucose dose and insulin release, and the time course profiles describe the dynamic first phase and the sustained steady-state release of biphasic insulin secretion. The ATP and ADP concentrations of the β-cell increase and decrease, respectively, with the external glucose dose, resulting in an increased ATP/ADP ratio that triggers insulin release. The model is able to reproduce the fast initial insulin release in the first phase and the steady-state insulin secretion in the second phase depending on glucose concentration. For the second phase, the constants of the Hill function were parameterized to fit the normalized data of adenine nucleotide ratio and steady-state insulin release rates. Experimental observations suggest that the parameters of the response function, such as the slope of the response function and the half-maximal response, can vary between animal species due to differences in the expression levels of glucose transporter (64). In this study, the data corresponds to both human and murine islets.

FIGURE 6

Figure 6 (A) Effect of variations in blood glucose on insulin secretion. The plot is analogous to Figure 5C. Data from 9, 13, 14, 16, 18, 36, 42, 45–47, 50–56, 58–62, 64). (B) Sensitivity analysis indicating the effect of perturbation in model parameters on insulin secretion. Heatmap illustrating the values of scaled local sensitivities illustrating the effect of parameter perturbations on the amount of insulin secretion at varying glucose doses. Highly sensitive values are colored in red and blue. The parameters which cause less than 1% change in insulin response for 10% perturbation were not displayed for clarity. For more details, please refer to Sec. 2.6. (C) Effect of change in model parameters on insulin secretion as a function of glucose dose. The rate of insulin secretion in response to perturbation in the values of ATPconsumption_Vm, HEX1_Vm, IRS_Katp_ratio, IRS_hillKatp_ratio. The vertical line indicates the model value.

2.6 Sensitivity analysis of parameters affecting GSIS

To determine how the model parameters affect the rate of insulin release, we performed a local sensitivity analysis (67). Figure 6B shows the sensitivity of insulin flux to a 10% change in model parameter values at different glucose concentrations. The rate of insulin secretion depends on the ATP/ADP ratio, so perturbing parameters that affect ATP formation and consumption has strong effects. Figure 6C shows the highly sensitive parameters that have positive and negative effects on insulin secretion, including factors affecting ATP consumption, ATP formation, hexokinase, phosphofructokinase, and ATP/ADP-dependent insulin secretion.

In conclusion, our systematic pathway modeling workflow provides insights into the key mechanisms of GSIS in the pancreatic β-cell.

3 Discussion

We have developed a comprehensive kinetic model of GSIS in the pancreatic β-cell that can simulate glucose-dependent changes in glycolytic intermediates, ATP/ADP ratio, and their effect on insulin secretion. The main objective of this study was to establish a standardized workflow for data integration and normalization to construct a tissue-specific model of glycolysis and GSIS in the β-cell. Although we did not model other important pathways related to ATP homeostasis, such as the citric acid cycle, the pentose phosphate pathway, and the respiratory chain, our workflow can be easily extended to include them. Incorporating these pathways into our model will enable us to explicitly model the regulatory effect of downstream metabolites on the ATP/ADP ratio and insulin secretion. Previous studies have shown that fatty acids and amino acids can also induce insulin secretion in addition to glucose (4, 17, 68–71). Therefore, linking glucose metabolism with fatty acid and amino acid metabolism could help in understanding the insulinotropic effects of other fuel sources.

The increase in ATP levels triggers a cascade of events that culminate in the release of insulin from β- cells. Precisely, high ATP levels prompt the closure of ATP-sensitive potassium channels (6). Consequently, the cell membrane depolarizes, opening voltage-gated calcium channels, which allows calcium influx. The influx of calcium triggers exocytosis of insulin-containing vesicles, leading to the release of insulin into the bloodstream (7, 8). A biphasic time course of insulin secretion is observed in vitro and in vivo studies, with a rapid initial phase followed by a sustained steady-state release. The biphasic release is attributed to the time-dependent formation, translocation, exocytosis of insulin granules, and decline in the amplitude of action potential contributing to the termination of first-phase insulin secretion (7, 72). These electrophysiological changes resulting in insulin secretion were not modeled explicitly, but a minimal model was used to capture first-phase insulin secretion and the effect of the ATP/ADP ratio on insulin secretion. The model’s predictive capacity is limited to the biphasic glucose-insulin secretion dynamics. Expanding the model to explicitly describe the β-cell electrophysiology would allow us to study experimentally observed patterns such as pulsatile insulin secretion (23).

To summarize, the advancements presented in our work were employed to study GSIS in the pancreatic β-cell. While the existing models have certain limitations, they also have strengths and features that our model does not include, such as the electrophysiology of calcium handling and insulin granule dynamics (28, 73), the pulsatile insulin response to glucose (23, 74), and the potentiating effects of other fuel sources on insulin secretion (4, 17). Although our model has limitations, it represents the first data-driven approach to integrate information from diverse sources and experimental setups. Moreover, it provides the first systematic analysis of the glycolytic changes that occur during insulin secretion in response to different glucose levels. Our study reveals that in GSIS, almost all glycolytic intermediates increase in a glucose-dependent manner as do total ATP and NADH, which is a significant finding.

Model predictions deviate from the dataset for some glycolytic intermediates, despite incorporating condition-specific experiments with pre-incubation and incubation glucose doses in the model parameterization. Possible reasons for the deviations in the time course and steady-state model predictions of fructose 6-phosphate (F6P), fructose 2,6-bisphosphate (F26BP), glyceraldehyde-3-phosphate (GRAP), and dihydroxyacetone phosphate (DHAP) from the experimental data include the following. For species such as F26BP, the time course data was obtained from a single study at a specific incubation glucose dose. We observed that for the initial concentration specified in the model, the concentration of F26BP and F6P saturates to a value higher than that observed in the normalized dataset of time-course experiments. Since F6P and F26BP are involved in the same reaction, an offset in one metabolite has an effect on the other. Therefore, better initial concentrations can only be defined if additional data are available at different combinations of glucose pre-incubation and incubation doses. In addition, the fit can be improved by improving the kinetics associated with the conversion of F6P to F26BP. The dynamics of these upper glycolytic intermediates may also be influenced by other pathways not modeled in the current study. Lower glycolysis reactions are sequential, and missing data for an intermediate species may affect the fluxes involved in the reaction chain that forms or consumes the intermediate (i.e., a fast first step and a slow second step and vice versa are the same). The fluxes of reactions such as glyceraldehyde-3-phosphate dehydrogenase and phosphoglycerate kinase could potentially be affected due to the lack of steady-state and time-course data for 1,3-biphosphoglycerate (BPG), and the estimates of the associated parameters may not be optimal. Therefore, the uncertainty in the prediction of BPG can lead to a deviation in the prediction of the concentrations of GRAP and DHAP, which are in equilibrium.

The comparative analysis, which shows the dynamics of the model output and the experimental data (Figures 5A–C, right panel), was performed for combinations of the experimental pre-incubation and incubation glucose doses. When the pre-incubation and incubation glucose doses are the same (e.g., in the study by Taniguchi et al. (12), islets were pre-incubated in glucose at 2.8 mM followed by incubation at 2.8 mM, low glucose dose), the system evolves and saturates in the pre-incubation phase. Consequently, steady-state behavior is observed in time-course profiles plotted at incubation conditions. Results are presented for incubation because data points curated from experimental studies correspond to incubation conditions.

The initial concentrations observed in the experimental dataset differ between experiments due to complex experimental designs, such as the islet cells subjected to different pre-incubation and incubation conditions. Moreover, the experimental time course profiles of most of the metabolic intermediates are only available at two incubation conditions, low and high glucose doses (12–14, 18). Therefore, only when experimental time course data at intermediate glucose doses are available, better initial states can be defined to predict the steady-state response of the system accurately. Therefore, in specific experimental cases, appropriate model extensions, improved kinetics, and data integration into our standardized workflows may be needed to refine the model and deal with experimental uncertainties. Otherwise, the simulation explains well the predictions of other metabolic intermediates, thus demonstrating the correctness of our methodology.

Existing models (17, 27, 30, 75) often suffer from several drawbacks such as limited evaluation to a single data set, non-standardized formats of experimental data and kinetic parameters, and non-reproducible formats. To overcome these limitations, we have created open, free, and FAIR assets that can be used for the study of pancreatic physiology and GSIS. These assets include a fully reproducible SBML model of pancreatic β-cell glycolysis, a data curation workflow, strategies for unit and data normalization, and a large database of metabolic data of the pancreatic β-cell. Our model can be extended further to study glucose-insulin regulatory mechanisms by for example integrating calcium handling in β-cell (28, 29, 76) or paracrine regulation of the release of counter-regulatory hormones such as glucagon by gap-junctional coupling of α-, β-, and δ-cells present in the islets (77, 78). Overall, our systematic model-building workflow can be used as a blueprint to construct reproducible kinetic models of cell metabolism.

Computational modeling faces a significant challenge due to the substantial variation in data across different experimental systems, species, and cell lines. Often, relative data instead of absolute data is reported, further complicating the task of data integration. In this study, we developed a reliable data normalization workflow that was applied to experimental and clinical data from 39 studies conducted over the past 50 years on pancreatic, islet, and β-cell function in various species and cell lines. Our approach substantially reduced data heterogeneity and revealed a highly consistent response in glycolytic metabolites and insulin secretion. The high degree of conservation in the system of GSIS may have contributed to the effectiveness of the normalization workflow, as similar mechanisms are at play in different species, and the general changes can be observed across various experimental systems.

The study has laid a strong groundwork for enhancing our comprehension of the underlying reasons behind impaired insulin secretion. By mapping proteomics or transcriptomics data onto specific pathways, the developed model could be utilized to gain further insight into changes in GSIS, for instance in diabetic patients.

Furthermore, this model can serve as a crucial component for physiological whole-body models of glucose homeostasis, allowing researchers to investigate the relationship between the potentiation of insulin release and glucose uptake by insulin-responsive tissues. Evidence suggests that, in addition to nutrient-secretagogues, hormone potentiators such as incretins contribute to 74% (79, 80) of postprandial insulin secretion. The role of gut hormones such as incretins on insulin biosynthesis, insulin secretion and their effect on β-cell mass can be studied by integrating the subcellular model developed in our study with whole-body models (81, 82). Incretins bind to the receptors on β-cells and regulate the ion channels through signaling mechanisms that augment glucose-stimulated insulin secretion (83, 84). For example, the insulinotropic effect of endogenously secreted incretin hormones such as glucose-dependent insulinotropic polypeptide (GIP) and glucagon-like peptide-1 (GLP-1) (83), exogenously administered incretin mimetics such as exenatide and liraglutide on postprandial insulin secretion and renal elimination rates of these antidiabetic drugs can be examined (85–89).

The data presented in this study was obtained from experiments where the incretin effect did not play any role. Specifically, the experiments involved islet or cell studies where glucose was systematically varied and controlled. However, if the model is to be applied in a more physiological context, such as in a physiological-based pharmacokinetics model of glucose regulation, it is essential to extend the model to include the incretins. This is particularly important if the focus of the model is to describe glucose-stimulated insulin secretion (GSIS) in the context of oral glucose tolerance tests or meal challenges.

In conclusion, this study utilized a systematic modeling workflow to gain insight into the key mechanisms involved in glucose-stimulated insulin secretion (GSIS) in pancreatic β-cells. Crucially, by establishing a standardized workflow for data integration, normalization, and data fitting, our approach allows for easy incorporation of additional data sets and re-fitting of the model to extend the scope of the model. When extended for translational purposes in clinical settings, it can serve to create reference models to identify variations in subjects which can lead to useful inferences regarding underlying metabolic conditions with therapeutic relevance.

4 Methodology

The workflow for building the kinetic model is illustrated in Figure 2, with the following sections providing information on the individual steps.

4.1 Stoichiometric model

Chemical formulas and charges were assigned to all metabolites, and reactions were examined to ensure that they maintained mass and charge balance. The kinetic model encompasses glycolytic reactions and correlates the energy status of the β-cell with insulin secretion. sbmlutils (90) was used to create and validate the model, while cy3sbml (66) was used to confirm its coherence. sbmlutils is a collection of python utilities for working with SBML models and cy3sbml is a java-based SBML plugin for Cytoscape (91) used for visualization of SBML models. The mass and charge balance of the system was verified using cobrapy (92).

4.2 Metadata integration

Adding semantic annotations to models is an essential aspect of improving their interoperability and reusability, as well as facilitating data integration for model validation and parameterization (93, 94). To describe the biological and computational significance of models and data in a machine-readable format, semantic annotations are encoded as links to knowledge resource terms. Open modeling and exchange (OMEX) metadata specifications were employed to annotate model compartments, species, and reactions with metadata information (Figure 2B).

4.2.1 Case study: phosphoglycerate kinase

The enzyme phosphoglycerate kinase (PGK) catalyses the conversion of 1,3-biphosphoglycerate (bpg13) and ADP to form 3-phosphoglycerate (pg3) and ATP.

$$adp+bpg13\rightleftharpoons atp+pg3$$

In our model, PGK is described by the following annotations: SBO:0000176, vmhreaction/PGK, bigg.reaction/PGK, kegg.reaction/R01512, ec-code/2.7.2.3, biocyc/META : PHOSGLYPHOS-RXN, uniprot:P00558, uniprot:P07205.

The model components, including physical volumes, reactions, metabolites, and kinetic-rate laws, were annotated using Systems Biology Ontology (SBO) terms, which describe the computational or biological meaning of the model and data (95). Biomedical ontology services such as Ontology Lookup Service (OLS) (96), VMH (97), and BiGG (98) were used to collect these terms. Additional information for species and reactions were gathered from various databases such as HMDB, BioCyc, MetaNetX, ChEBI, and SEED. For instance, the model’s metabolites were annotated with identifiers from VMH, BiGG, KEGG, HMDB, BioCyc, ChEBI, MetaNetX, and SEED, while reactions were annotated with VMH, Rhea, MetaNetX, SEED, BiGG, BioCyc, and KEGG identifiers (99). Enzymes catalyzing reactions were annotated with identifiers from enzyme commission (EC) numbers, UniProt (100), and KEGG. Finally, the annotations were incorporated into the SBML file using 707 sbmlutils (90) and pymetadata (101).

4.2.2 Case study: 1,3-biphosphoglycerate

There is currently a bottleneck in data integration due to the use of multiple synonyms to refer to a single compound in data repositories. For instance, bpg13 is identified by different names in SABIO-RK (Glycerate 1,3-bisphosphate, 3-phospho-D-glyceroyl phosphate) and BRENDA (3-phospho-D-glyceroyl phosphate). Additionally, the labeling of 1,3-biphosphoglycerate, abbreviated as DPG, varies across existing β-cell models (e.g., 1,3-bisphospho-D-glycerate in (75) and 1,3-biphosphoglycerate in (17). Overall, bpg13 is associated with seven synonyms: 1,3-Bisphospho-D-glycerate, 13dpg, 3-Phospho-D-glyceroylphosphate, Glycerate 1,3-bisphosphate, 3-phospho-d-glyceroyl-phosphate, 1,3-diphosphoglyceric acid, 3-Phospho-D-glyceroyl phosphate. This issue makes it difficult to integrate data and information from different resources, highlighting the need to link chemical entities in the model to knowledge resource terms.

In our model, bpg13 is clearly described by the following metadata annotations: SBO:0000247, vmhmetabolite/13dpg, bigg.metabolite/13dpg, biocyc/METADPG, kegg.compound/C00236, CHEBI:16001, inchikey:LJQLQCAXBUHEAZ-UWTATZPHSA-N.

The formula and charge of bpg13 are C3H4O10P2 and -4, respectively.

4.3 Kinetic parameters

Kinetic parameters, such as half-saturation constants (K_M), inhibition constants (K_I), activation constants (K_A), and equilibrium constants (K_eq), were gathered from literature and a variety of databases (see Figure 2C). Values were programmatically accessed from UniProt (100), BRENDA (102) using brendapy (103), and SABIO-RK (104). These databases were searched using an organism’s NCBI taxonomy identifier and reaction EC number as input search terms. Various parameters, including measurement type (K_m, Ki, and K_a), experimental conditions (pH, temperature), KEGG reaction identifiers, enzyme type (wildtype or mutant), associated metabolite identifiers (SABIO compound name or BRENDA ligand id), UNIPROT identifiers associated with the isoforms of an enzyme, source tissue, and details of data source (PubMed identifier) were obtained. Since there is limited availability of kinetic data for Homo sapiens, we also searched for parameter values reported in studies of animal species that are closely related to humans and utilized them if no data were available for humans.

4.4 Synonym mapping

To map compound synonyms associated with each queried metabolite, we utilized compound identifier mapping services and available metadata annotations. First, we associated the name of each compound with internal database identifiers, such as the internal identifier of Glycerone-phosphate in SABIO, which is 28. Then, we linked the internal identifiers to external identifiers, such as those from ChEBI and KEGG. The external identifiers associated with the SABIO ligand identifier were obtained from cross-ontology mappings available in SABIO-RK. Similarly, we queried the REST API of UniChem to obtain the external identifiers associated with the BRENDA ligand identifier. By doing so, we were able to map most of the kinetic parameters to their respective compounds (Figure 2D).

4.5 Model parameters

For each parameter in the model, the median value was calculated after synonym mapping and the values were assigned to the model parameters, see Figure 2E. This was performed for initial concentrations, equilibrium K_eq constants, half-saturation constants K_m, inhibition K_i, and activation K_a constants.

4.6 Data curation

The next step involved curating data from studies that reported metabolite values, insulin secretion, or maximal velocities of glycolytic reactions V_max in pancreatic, islet, and β-cell lines (Figure 2F). Our search for the studies used in model development was performed by using any combination of the following words: “glycolytic intermediates”, “metabolite profiling”, “concentration measurements”, “time course”, “glucose-dependence”, “pancreatic β-cell”, “pancreatic islets”, “endocrine pancreas”, “glucose-stimulated insulin secretion”, “fuel-stimulated insulin secretion”, “insulin response” and the name of the metabolite or the name of adenine and pyridine nucleotides in the search string. Relevant studies were identified through a literature search in PubMed, with a focus on time course and dose-response profiles of metabolite concentrations for metabolites and insulin secretion. Tissue homogenates were prepared by isolating islets from rodents, humans, or insulin-secreting cell lines (see Table 1). Assays were performed by stimulating the medium with various pre-incubation and incubation concentrations of glucose. To curate the data, established curation workflows from PK-DB (105), which were applied in a recent meta-analysis (106), were used. The numerical data was digitized by extracting the data points from the figures and tables using WebPlotDigitizer (107). The incubation time and glucose concentration of the stimulation medium were recorded for all measurements, and meta-information such as organism and tissue type were documented.

The data is available under a CC-BY 4.0 license from https://github.com/matthiaskoenig/pancreas-model. In this study, version 0.9.6 of the data set is used (108, 109).

4.7 Unit normalization

The data measured in different studies is often reported in different units. Therefore, unit normalization was performed to integrate the data and convert metabolite concentrations and insulin secretion to standardized units of mmole/l (mM) and nmole/min/ml (β-cell volume), respectively (Figure 2G).

Absolute measurements reported in metabolic profiling studies were found in various units such as per gram DNA, per gram wet weight or dry weight of the islet tissue, per cell, per islet, etc. To use these values for model calibration, both the absolute and relative measurements were first converted to concentration units in mM. The absolute values were converted to model units by multiplying the raw values with appropriate unit conversion factors. For instance, the islet content of glucose 6-phosphate, G6P, (pmol/islet) was converted to concentration units (mM) using the distribution volume of water in the islet (2nl/islet) (35) as the conversion factor. Relative measurements were mainly reported with reference to a basal concentration. These relative measurements were converted to absolute quantity by multiplying the fold values with the respective metabolite concentration at the basal or pre-incubation concentration of glucose.

4.8 Data normalization and integration

Data collected from experiments performed in different laboratories, under different experimental conditions, and with different animal species showed significant variability after unit normalization. Therefore, data normalization was performed to eliminate systematic discrepancies between data reported in different studies (as shown in Figure 2H). To achieve this, least squares approach was used to minimize the distance between individual experimental curves and the mean curve, which is the weighted average of all curves for a given metabolite. The data normalization process involved a two-step procedure in which the steady-state data were first normalized for each metabolite. The resulting steady-state normalization was then used to normalize the time course data for that metabolite (see Figure 4 for the example of glucose-6 phosphate).

4.8.1 Steady-state data normalization

Steady-state (ss) experiments consisted of pre-incubation with one glucose dose followed by incubation with another glucose dose. The steady state data of the experiment α, $(c_0^{\alpha },c_1^{\alpha },\dots ,c_n^{\alpha })$ observed at n incubation glucose doses $(d_0^{\alpha },d_1^{\alpha },\dots ,d_n^{\alpha })$ is expressed by the piecewise linear-interpolation function ${\mathcal{C}}^{ss}$. Here, α belongs to the set of steady-state experiments $1\le {\alpha }^{}$ $\le N^{{\alpha }_{ss}}$ with $N^{{\alpha }_{ss}}$ being the number of steady-state experimental curves of the metabolite s.

Mean curve. The mean steady-state curve $\overline{{\mathcal{C}}^{ss}}$ of each metabolite s is calculated as the weighted average of all experimental curves. The data points of the mean curve were interpolated using a piecewise smooth spline function. For data sets consisting of 2 data points, a linear interpolation was used.

We formulate a least-squares problem to minimize the distance between the individual experimental curves and the mean curve $\overline{{\mathcal{C}}^{ss}}$. The cost function F of the minimization problem is given by,

$\begin{array}{ll}F(f^{\alpha })={\displaystyle \sum _{i=1}^n}{(f^{\alpha }·{\mathcal{C}}^{ss}(d_i^{\alpha })-\overline{{\mathcal{C}}^{ss}}(d_i^{\alpha }))}^2& (1)\end{array}$

In Eq. 1, ${\mathcal{C}}^{ss}(d_i^{\alpha })$and $\overline{{\mathcal{C}}^{ss}}(d_i^{\alpha })$ are the function values of the individual and mean interpolation function at the i^th value of the glucose dose. N is the number of glucose values in the dose-response curve of the experiment α.

For each experimental curve, the factor $f^{\alpha }$ was determined so that the residual error in Eq. 1 is minimized. The residual error is minimum at the point where the derivative of the cost function F is zero. Taking the partial derivative of Eq. 1 with respect to the scale transformation parameter gives factor f^α of the experimental curve α (Eq. 2).

$\begin{array}{ll}{f}^{\alpha }=\frac{{\sum }^{ }{\mathcal{C}}^{ss}(d_i^{\alpha })\overline{{\mathcal{C}}^{ss}}(d_i^{\alpha })}{{\sum }^{ }{({\mathcal{C}}^{ss}(d_i^{\alpha }))}^2}& (2)\end{array}$

The scale factors of all steady state curves $(f^1,\dots f^{N^{{\alpha }_{ss}}})$ were determined by minimizing the respective cost functions $(F(f^1),\dots F(f^{N^{{\alpha }_{ss}}})$. Multiplying the experimental curve C^α by the scaling factor $f^{\alpha }$ shifts the experimental curve towards the mean curve. A new mean curve can be calculated with the scaled data. The curves were scaled iteratively until all $f^{\alpha }$ converged.

The scale transformation factors are chosen by minimizing the variance with respect to the mean of observations (Eq. 1), which is the conditional expectation given a set of observations (110). The minimal variance estimate is the optimal estimate given a set of observations and this results in smoothing the noise in the data.

4.8.2 Time course data normalization

Time course (tc) experiments consisted of pre-incubation with one glucose dose followed by incubation with another glucose dose. The time-dependent data of the time course experiment $\beta (c_0^{\beta },c_1^{\beta },\dots ,c_m^{\beta })$ observed at m time points $(t_0^{\beta },t_1^{\beta },\dots ,t_m^{\beta })$ is expressed by the piecewise linear-interpolation function C^β. Here, β belongs to the set of time course experiments $1\le \beta \le N^{{\beta }_{tc}}$ with $N^{{\beta }_{tc}}$ being the number of time course experimental curves of the metabolite s. For normalization, each time course was scaled by a factor f^β.

For a given incubation glucose dose d^β, the metabolite concentration at the last time point ${\mathcal{C}}^{tc}(t_m)$ corresponds to the steady state value reached for the given d^β:

$\begin{array}{ll}{f}^{\beta }·{\mathcal{C}}^{tc}(t_m)-\overline{{\mathcal{C}}^{ss}}(d^{\beta })=0& (3)\end{array}$

The scaling factor for the time course experiment follows as:

$\begin{array}{ll}{f}^{\beta }=\frac{\overline{{\mathcal{C}}^{ss}}(d^{\beta })}{{\mathcal{C}}^{tc}(t_m)}& (4)\end{array}$

4.9 Model inputs

The SBML model was generated by specifying initial concentrations, rate expressions, parameter values, and compartmental volumes as the model inputs, see Figure 2I.

Volume. The physical volume of the cytoplasmic compartment and the β-cell volume were obtained from the values reported in a morphometric study of β-cells (111).

Initial concentrations. The initial concentrations of glycolytic intermediates and adenine nucleotides were obtained from the mean curve $\overline{{\mathcal{C}}^{ss}}$ (Sec. 2.1) at a basal glucose concentration of 3 mM. The initial value of glucose in the external/blood compartment is 3 mM.

The initial concentrations of cofactors, phosphate and pyridine nucleotides, were expressed as polynomial functions passing through the data points of the mean curve, which is computed as the weighted average of data normalized experimental curves (Sec. 2.1). In the SBML model, the polynomial expressions were defined using assignment rules.

Kinetic constants. The median values of the half-saturation or Michaelis-Menten constants K_m (Sec. 4.5), were assigned to the model parameters.

Equilibrium constants. The values of the equilibrium constants K_eq were collected from NIST (112) and EQUILIBRATOR (113).

Model equations. For all the glycolytic reactions, the biochemical interactions were expressed using modular rate laws (114) of the form Eq. 5.

$\begin{array}{ll}\nu =\frac{V_{max}{\displaystyle \prod _i}{a}_i(1-\frac{\Gamma }{Keq})}{{\displaystyle \prod _i}(1+a_i)+{\displaystyle \prod _j}(1+b_j)-1}& (5)\end{array}$

Here, $a_i$ is $S_i/K{m}_s$, $b_i$ is $P_i/K{m}_p$, S refers to the substrate and P refers to the product. $K_{eq}$ is the equilibrium constant and $\text{Γ}$ is the mass-action ratio (114). The use of detailed mechanistic rate laws was avoided due to the challenges associated with finding a large number of parameter values.

Biphasic insulin secretion in response to elevated glucose levels and change in the energy state of the β-cell was modeled as the sum of two components, a dynamic first phase and a static second-phase insulin profile (21, 23, 115). In Eq. 6, the first phase ${\text{IRS}}_{fp}$ accounts for the rapid rise in insulin. We model this using a function proportional to the rate of change in ATP. The second phase ${\text{IRS}}_{sp}$, which captures the sustained steady-state release, was modeled via a phenomenological equation depending on ATP/ADP ratio. The insulin release flux given by Eq. 6, is characterized by five parameters, the proportionality constant $k_{d_{fp}}$ of the first phase insulin release, the maximal rate of the second phase $V_{max}^{IR{S}_{sp}}$ insulin release, $K_{m_{fp}}$ half-maximal constant of the first phase and $K_{m_{sp}}$ the ratio of ATP/ADP that results in half-maximal insulin release of the second phase, and the Hill coefficient $n_{sp}$ of second phase insulin release.

$\begin{array}{ll}{\nu }^{IRS}=\underset{{\text{IRS}}_{fp}(\text{t})}{\underbrace{k_{d_{fp}}\frac{\text{max} (\frac{dATP}{dt},0)}{\text{max} (\frac{dATP}{dt},0)+K_{m_{fp}}}}}+\underset{{\text{IRS}}_{sp}(\text{t})}{\underbrace{V_{max}^{IR{S}_{sp}}\frac{{(\frac{ATP}{ADP})}^{n_{sp}}}{{(\frac{ATP}{ADP})}^{n_{sp}}+K_{m_{sp}}^{n_{sp}}}}}& (6)\end{array}$

Boundary metabolites and reactions. Species in the external and mitochondrial compartments were assumed to be boundary species with constant concentrations, i.e. glucose and lactate in the external compartment and pyruvate in the mitochondrial compartment were held constant. Some boundary reactions were modeled as irreversible reactions, i.e. the export of lactate and the transport of pyruvate in the mitochondrion.

Metabolites determined by rate rules. To account for glucose-dependent changes in the concentrations of phosphate, NAD, and NADH, polynomial functions were used to express the concentrations as rate rules. This approach ensured that the concentration of fixed metabolites in the system increased as a function of glucose dose.

Changes in total adenine nucleotides. The sum of adenine nucleotides ($ATP+ADP=AT{P}_{tot}$) changes with glucose. To account for these changes, a reaction $\text{Δ}$ATP was added that changes the total ATP according to the observed steady-state data for a given glucose value (Eq. 7).

$\begin{array}{ll}\text{Δ}ATP=f(AT{P}_{tot}(glc)-(ATP+ADP))& (7)\end{array}$

The $AT{P}_{tot}(glc)$ values are determined by the interpolating polynomial of the mean steady-state glucose dose response of the ATP+ADP data.

4.10 Model calibration

The normalized time-course data was used for model calibration and parameter estimation (Figure 2J). An overview of the subset of data used for model calibration is shown in Figure 1. The following data were not used: NADH and NAD were fixed metabolites in the model, with NAD/NADH and NADH+NAD calculated from the metabolites. Total ATP was calculated by summing ATP and ADP, and ATP ratio was calculated by finding the ratio. The insulin secretion rate (IRS) was used to derive the parameters of the IRS function.

A subset of the $V_{max}$ parameters was optimized to minimize the error between model predictions and experimental observations. The cost function is given by the sum of squares of residuals

$\begin{array}{ll}F(\mathit{P})={\displaystyle \sum _{\alpha ,s}}{({\mathit{c}}_s^{\alpha }-{\mathit{c}}_s^{ℳ}(\mathit{P}))}^2& (8)\end{array}$

In Eq. 8, ${\mathit{c}}_s^{\alpha }$ is the concentration of the metabolite $s$ in the experiment $\alpha$ and ${\mathit{c}}_s^{ℳ}$ is the concentration of the metabolite $s$ predicted by the model $ℳ$. P is the set of 16 parameters of maximum reaction rates $V_{max}$. The experimental data of all transient metabolites in the model were stored in spreadsheets. The parameter estimation simulation experiments were set up using basiCO (116), the Python interface of COPASI (117).

To enable the simulation of experimental setups such as pre-incubation and incubation conditions, the corresponding glucose doses were curated from experimental studies. We perform condition-specific model simulations by running pre-simulations using the pre-incubation glucose dose for 60 minutes. Pre-simulation or pre-equilibration at given conditions is a task often performed during model simulation or parameter optimization (118, 119). Following pre-simulation, the system was subjected to the simulation phase at the incubation glucose dose for the duration indicated in the experimental studies. We set up the pre-simulation and simulation phases for the parameter estimation task using Events. The pre-incubation and incubation glucose concentrations were mapped to the independent variable ($gl{c}_{ext}$, glucose in the external compartment), and incubation time was mapped to model time. The transient metabolites were assigned to the model elements as dependent variables. The mean values of $V_{max}$calculated from the curated values of the enzyme activities were assigned as initial values. The lower and upper bounds specified for the reaction rates $V_{max}$ were set to 1e-2 and 5000, respectively. When zero was used as the lower bound, the global optimization resulted in parameter sets for which reaction fluxes were close to equilibrium (i.e., zero or negligible flux). For a high upper bound value (10000), we observed that the concentration profiles rise to saturation faster, possibly due to the high $V_m$ values of the reactions.

The calculations were performed using Cloud-COPASI, the front-end to a computer cluster at the Centre for Cell Analysis and Modelling. Cloud-COPASI is an extension of Condor-COPASI (120). We carried out a hybrid optimization approach (121), following the global optimization a local optimization was performed. 100 iterations of parameter estimation were performed with random initial guesses on Cloud-COPASI using Evolutionary Strategy (SRES), a global optimization method (121–124). The parameter set obtained from the iteration that yielded the minimum objective value and steady-state was updated in the model. The system was then subjected to a local optimization run using the Hooke and Jeeves algorithm to obtain the optimal estimate.

4.11 Kinetic model and model predictions

All information was written into the model, validation was performed using sbmlutils, and model simulations were performed, see Figures 2K, L.

Finally, we performed model predictions of glycolytic intermediates and insulin response as a function of varying glucose concentrations. The set of differential equations was numerically integrated using basiCO (116) based on COPASI (117) and sbmlsim (125) based on libroadrunner (126, 127). Pre-simulations were performed by simulating the model with optimal parameter values at a pre-incubation glucose dose of 3 mM for 60 minutes. For the time course simulations, glucose was varied as linspace (1, 35, num=11), and simulations were run for 60 minutes. For the glucose dose-response, glucose was varied identically, and the model was simulated to steady-state. To compare the dynamics of the model predictions and the experimental data, simulations were performed using the combinations of the experimental pre-incubation and incubation glucose doses. The time course predictions presented in Sec. 3.5 correspond to the simulation phase. Simulations were performed either with COPASI or independently using libroadrunner to ensure reproducibility of key model results.

The model is available in SBML (33, 128) under a CC-BY 4.0 license from https://github.com/matthiaskoenig/pancreas-model. In this study, version 0.9.6 of the model is presented (108, 109).

Data availability statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.

Author contributions

DM, SR, MK, and DP conceived and designed the study. DM and MK developed and implemented the computational model and data normalization workflow, and performed the analysis. DM curated the experimental data, performed parameter estimation, and drafted the initial version of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

Research of DM was supported by the Senior Research Fellowship from the Ministry of Human Resource Development (MHRD), Government of India. MK was supported by the Federal Ministry of Education and Research (BMBF, Germany) within the research network Systems Medicine of the Liver (LiSyM, grant number 031L0054) and ATLAS (grant number 031L0304B) and by the German Research Foundation (DFG) within the Research Unit Program FOR 5151 “QuaLiPerF (Quantifying Liver Perfusion-Function Relationship in Complex Resection - A Systems Medicine Approach)” by grant number 436883643 and grant number 465194077 (Priority Programme SPP 2311, Subproject SimLivA). This work was supported by the BMBF-funded de.NBI Cloud within the German Network for Bioinformatics Infrastructure (de.NBI) (031A537B, 031A533A, 031A538A, 031A533B, 031A535A, 031A537C, 031A534A, 031A532B).

Acknowledgments

DM thanks Dr. Murthy Madiraju S.R., Montreal Diabetes Research Center, CRCHUM, Montréal, Canada, Dr. Pedro Mendes, University of Connecticut, and Dr. Frank Bergmann, University of Heidelberg for the invaluable discussions and incredible support with basiCO. Access to Cloud-COPASI is supported by NIH Grant R24 GM137787 from the National Institute for General Medical Sciences.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Woods SC, Lutz TA, Geary N, Langhans W. Pancreatic signals controlling food intake; insulin, glucagon and amylin. philosophical transactions of the royal society of London. Ser B Biol Sci (2006) 361:1219–35. doi: 10.1098/rstb.2006.1858

CrossRef Full Text | Google Scholar

2. König M, Bulik S, Holzhütter H-G. Quantifying the contribution of the liver to glucose homeostasis: a detailed kinetic model of human hepatic glucose metabolism. PloS Comput Biol (2012) 8:e1002577. doi: 10.1371/journal.pcbi.1002577

PubMed Abstract | CrossRef Full Text | Google Scholar

3. MacDonald PE, Joseph JW, Rorsman P. Glucose-sensing mechanisms in pancreatic beta-cells. philosophical transactions of the royal society of London. Ser B Biol Sci (2005) 360:2211–25. doi: 10.1098/rstb.2005.1762

CrossRef Full Text | Google Scholar

4. Prentki M, Matschinsky FM, Madiraju SM. Metabolic signaling in fuel-induced insulin secretion. Cell Metab (2013) 18:162–85. doi: 10.1016/j.cmet.2013.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Newsholme P, Cruzat V, Arfuso F, Keane K. Nutrient regulation of insulin secretion and action. J Endocrinol (2014) 221:R105–20. doi: 10.1530/JOE-13-0616

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ashcroft FM. K(ATP) channels and insulin secretion: a key role in health and disease. Biochem Soc Trans (2006) 34:243–6. doi: 10.1042/BST20060243

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Rorsman P, Braun M. Regulation of insulin secretion in human pancreatic islets. Annu Rev Physiol (2013) 75:155–79. doi: 10.1146/annurev-physiol-030212-183754

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Guerrero-Hernandez A, Verkhratsky A. Calcium signalling in diabetes. Cell Calcium (2014) 56:297–301. doi: 10.1016/j.ceca.2014.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Guay C, Joly E, Pepin E, Barbeau A, Hentsch L, Pineda M, et al. A role for cytosolic isocitrate dehydrogenase as a negative regulator of glucose signaling for insulin secretion in pancreatic ss-cell. PloS One (2013) 8:e77097. doi: 10.1371/journal.pone.0077097

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Fritsche L, Weigert C, Häring H-U, Lehmann R. How insulin receptor substrate proteins regulate the metabolic capacity of the liver–implications for health and disease. Curr Medicinal Chem (2008) 15:1316–29. doi: 10.2174/092986708784534956

CrossRef Full Text | Google Scholar

11. Di Camillo B, Eduati F, Nair SK, Avogaro A, Toffolo GM. Leucine modulates dynamic phosphorylation events in insulin signaling pathway and enhances insulin-dependent glycogen synthesis in human skeletal muscle cells. BMC Cell Biol (2014) 15:9. doi: 10.1186/1471-2121-15-9

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Taniguchi S, Okinaka M, Tanigawa K, Miwa I. Difference in mechanism between glyceraldehyde- and glucose-induced insulin secretion from isolated rat pancreatic islets. J Biochem (2000) 127:289–95. doi: 10.1093/oxfordjournals.jbchem.a022606

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Spégel P, Sharoyko VV, Goehring I, Danielsson APH, Malmgren S, Nagorny CLF, et al. Time-resolved metabolomics analysis of β-cells implicates the pentose phosphate pathway in the control of insulin release. Biochem J (2013) 450:595–605. doi: 10.1042/BJ20121349

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Spégel P, Andersson LE, Storm P, Sharoyko V, Göhring I, Rosengren AH, et al. Unique and shared metabolic regulation in clonal β-cells and primary islets derived from rat revealed by metabolomics analysis. Endocrinology (2015) 156:1995–2005. doi: 10.1210/en.2014-1391

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Malaisse WJ, Hutton JC, Kawazu S, Sener A. The stimulus-secretion coupling of glucose-induced insulin release. metabolic effects of menadione in isolated islets. Eur J Biochem (1978) 87:121–30. doi: 10.1111/j.1432-1033.1978.tb12357.x

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Detimary P, Van den Berghe G, Henquin JC. Concentration dependence and time course of the effects of glucose on adenine and guanine nucleotides in mouse pancreatic islets. J Biol Chem (1996) 271:20559–65. doi: 10.1074/jbc.271.34.20559

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Salvucci M, Neufeld Z, Newsholme P. Mathematical model of metabolism and electrophysiology of amino acid and glucose stimulated insulin secretion: In vitro validation using a β-cell line. PloS One (2013) 8:e52611. doi: 10.1371/journal.pone.0052611

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Huang M, Joseph JW. Assessment of the metabolic pathways associated with glucose-stimulated biphasic insulin secretion. Endocrinology (2014) 155:1653–66. doi: 10.1210/en.2013-1805

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Tornheim K. Are metabolic oscillations responsible for normal oscillatory insulin secretion? Diabetes (1997) 46:1375–80. doi: 10.2337/diabetes.46.9.1375

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Bertram R, Sherman A, Satin LS. Metabolic and electrical oscillations: partners in controlling pulsatile insulin secretion. Am J Physiol-Endocrinol Metab (2007) 293:E890–900. doi: 10.1152/ajpendo.00359.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Toffolo G, Breda E, Cavaghan MK, Ehrmann DA, Polonsky KS, Cobelli C. Quantitative indexes of beta-cell function during graded up&down glucose infusion from c-peptide minimal models. American journal of physiology. Endocrinol Metab (2001) 280:E2–10. doi: 10.1152/ajpendo.2001.280.1.E2

CrossRef Full Text | Google Scholar

22. Toffolo G, Campioni M, Basu R, Rizza RA, Cobelli C. A minimal model of insulin secretion and kinetics to assess hepatic insulin extraction. American journal of physiology. . Endocrinol Metab (2006) 290:E169–76. doi: 10.1152/ajpendo.00473.2004

CrossRef Full Text | Google Scholar

23. Pedersen MG, Corradin A, Toffolo GM, Cobelli C. A subcellular model of glucose-stimulated pancreatic insulin secretion. philosophical transactions. Ser A Mathematical Physical Eng Sci (2008) 366:3525–43. doi: 10.1098/rsta.2008.0120

CrossRef Full Text | Google Scholar

24. Buchwald P. A local glucose-and oxygen concentration-based insulin secretion model for pancreatic islets. Theor Biol Med Model (2011) 8:20. doi: 10.1186/1742-4682-8-20

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Merrins MJ, Bertram R, Sherman A, Satin LS. Phosphofructo-2-kinase/Fructose-2,6-bisphosphatase modulates oscillations of pancreatic islet metabolism. PloS One (2012) 7:e34036. doi: 10.1371/journal.pone.0034036

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Bertram R, Pedersen MG, Luciani DS, Sherman A. A simplified model for mitochondrial ATP production. J Theor Biol (2006) 243:575–86. doi: 10.1016/j.jtbi.2006.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Westermark PO, Kotaleski JH, Björklund A, Grill V, Lansner A. A mathematical model of the mitochondrial NADH shuttles and anaplerosis in the pancreatic beta-cell. Am J Physiol Endocrinol Metab (2007) 292:E373–93. doi: 10.1152/ajpendo.00589.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Fridlyand LE, Philipson LH. Glucose sensing in the pancreatic beta cell: a computational systems analysis. Theor Biol Med Model (2010) 7:15. doi: 10.1186/1742-4682-7-15

PubMed Abstract | CrossRef Full Text | Google Scholar

29. McKenna JP, Dhumpa R, Mukhitov N, Roper MG, Bertram R. Glucose oscillations can activate an endogenous oscillator in pancreatic islets. PloS Comput Biol (2016) 12:e1005143. doi: 10.1371/journal.pcbi.1005143

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Gelbach PE, Zheng D, Fraser SE, White KL, Graham NA, Finley SD. Kinetic and data-driven modeling of pancreatic β-cell central carbon metabolism and insulin secretion. PloS Comput Biol (2022) 18:e1010555. doi: 10.1371/journal.pcbi.1010555

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Smolen P. A model for glycolytic oscillations based on skeletal muscle phosphofructokinase kinetics. J Theor Biol (1995) 174:137–48. doi: 10.1006/jtbi.1995.0087

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hucka M, Nickerson DP, Bader GD, Bergmann FT, Cooper J, Demir E, et al. Promoting coordinated development of community-based information standards for modeling in biology: the COMBINE initiative. Front Bioengineering Biotechnol (2015) 3:19 Hucka2015. doi: 10.3389/fbioe.2015.00019

CrossRef Full Text | Google Scholar

33. Keating SM, Waltemath D, König M, Zhang F, Dräger A, Chaouiya C, et al. SBML level 3: an extensible format for the exchange and reuse of biological models. Mol Syst Biol (2020) 16:e9110. doi: 10.15252/msb.20199110

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Matschinsky FM, Ellerman JE. Metabolism of glucose in the islets of langerhans. J Biol Chem (1968) 243:2730–6. doi: 10.1016/S0021-9258(18)93432-0

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Ashcroft SJ, Hedeskov CJ, Randle PJ. Glucose metabolism in mouse pancreatic islets. Biochem J (1970) 118:143–54. doi: 10.1042/bj1180143

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ashcroft SJ, Capito K, Hedeskov CJ. Time course studies of glucose-induced changes in glucose-6-phosphate and fructose-1,6-1171 diphosphate content of mouse and rat pancreatic islets. Diabetologia (1973) 9:299–302. doi: 10.1007/BF01221858

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Ashcroft SJ, Weerasinghe LC, Randle PJ. Interrelationship of islet metabolism, adenosine triphosphate content and insulin release. 1186 Biochem J (1973) 132:223–31. doi: 10.1042/bj1320223

CrossRef Full Text | Google Scholar

38. Matschinsky FM, Pagliara AS, Stillings SN, Hover BA. Glucose and ATP levels in pancreatic islet tissue of normal and diabetic rats. J Clin Invest (1976) 58:1193–200. doi: 10.1172/JCI108572

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Akhtar MS, Verspohl E, Hegner D, Ammon HP. 6-Phosphogluconate/glucose-1141 6-phosphate ratio in rat pancreatic islets during1142 inhibition of insulin release by exogenous insulin. Diabetes (1977) 26:857–63. doi: 10.2337/diab.26.9.857

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Sugden MC, Ashcroft SJ. Phosphoenolpyruvate in rat pancreatic islets:1676 a possible intracellular trigger of insulin release? Diabetologia (1977) 13:481–6. doi: 10.1007/BF01234500

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Malaisse WJ, Sener A. Glucose-induced changes in cytosolic ATP content in pancreatic islets. Biochim Et Biophys Acta (1987) 927:190–5. doi: 10.1016/0167-4889(87)90134-0

CrossRef Full Text | Google Scholar

42. Sener A, Hutton JC, Kawazu S, Boschero AC, Somers G, Devis G, et al. The stimulus secretion coupling of glucose-induced insulin release. metabolic and functional effects of NH4+ in rat islets. J Clin Invest (1978) 62:868–78. doi: 10.1172/JCI109199

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ammon HP, Hoppe E, Akhtar MS, Niklas H. Effect of leucine on the pyridine nucleotide contents of islets and on the insulin released–interactions in vitro with methylene1159 blue, thiol oxidants, and p-chloromercuribenzoate. Diabetes (1979) 28:593–9. doi: 10.2337/diab.28.6.593

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Ashcroft SJ, Christie MR. Effects of glucose on the cytosolic ration of reduced/oxidized nicotinamide-adenine dinucleotide phosphate in rat islets of langerhans. Biochem J (1979) 184:697–700. doi: 10.1042/bj1840697

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Trus MD, Hintz CS, Weinstein JB, Williams AD, Pagliara AS, Matschinsky FM. A comparison of the effects of glucose and acetylcholine on insulin release and intermediary1712 metabolism in rat pancreatic islets. J Biol Chem (1979) 254:3921–9. doi: 10.1016/S0021-9258(18)50675-X

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Trus M, Warner H, Matschinsky F. Effects of glucose on insulin release and on intermediary metabolism of isolated perifused1706 pancreatic islets from fed and fasted rats. Diabetes (1980) 29:1–14. doi: 10.2337/diab.29.1.1

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Ewart RB, Yousufzai SY, Bradford MW, Shrago E. Rat islet mitochondrial adenine 1296 nucleotide translocase and the regulation of insulin secretion. Diabetes (1983) 32:793–7. doi: 10.2337/diab.32.9.793

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Giroix MH, Sener A, Pipeleers DG, Malaisse WJ. Hexose metabolism in pancreatic islets. inhibition of hexokinase. Biochem J (1984) 223:447–53. doi: 10.1042/bj2230447

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Sener A, Van Schaftingen E, Van de Winkel M, Pipeleers DG, Malaisse-Lagae F, Malaisse WJ, et al. Effects of glucose and glucagon on the fructose 2,6-bisphosphate content of pancreatic islets 1645 and purified pancreatic b-cells. a comparison with isolated hepatocytes. Biochem J (1984) 221:759–64. doi: 10.1042/bj2210759

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Meglasson MD, Matschinsky FM. Pancreatic islet glucose metabolism and regulation of insulin secretion. Diabetes/Metabolism Rev (1986) 2:163–214. doi: 10.1002/dmr.5610020301

CrossRef Full Text | Google Scholar

51. Hedeskov CJ, Capito K, Thams P. Cytosolic ratios of free [NADPH]/[NADP+] and 1367 [NADH]/[NAD+] in mouse pancreatic islets, and nutrient-induced insulin secretion. Biochem J (1987) 241:161–7. doi: 10.1042/bj2410161

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Corkey BE, Glennon MC, Chen KS, Deeney JT, Matschinsky FM, Prentki M. A role for malonyl-CoA in glucose-stimulated insulin secretion from clonal pancreatic beta-cells. J of1228 Biol Chem (1989) 264:21608–12. doi: 10.1016/S0021-9258(20)88227-1

CrossRef Full Text | Google Scholar

53. Brun T, Roche EF, Corkey BE, Kim KH, Prentki M. Evidence for an Anaplerotic/Malonyl-Pathway in Pancreatic Beta-Cell Nutrient Signaling. Diabetes (1996) 45(2):190–8. doi: 10.2337/diab.45.2.190

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Ammon HP, Bacher M, Brändle WF, Waheed A, Roenfeldt M, El-Sayed ME, et al. Effect of forty-eight-hour glucose infusion into rats on islet ion fluxes, ATP/ADP ratio and redox ratios of pyridine nucleotides. J Endocrinol (1998) 156:583–90. doi: 10.1677/joe.0.1560583

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Detimary P, Dejonghe S, Ling Z, Pipeleers D, Schuit F, Henquin JC. The changes in adenine nucleotides measured in glucose-stimulated rodent islets occur in beta cells but not in alpha cells and are also observed in human islets. The Journal of Biological Chemistry (1998) 273(51):33905–8. doi: 10.1074/jbc.273.51.33905

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Liu YQ, Tornheim K, Leahy JL. Fatty acid-induced beta cell hypersensitivity to glucose. increased phosphofructokinase activity and lowered glucose-6-phosphate content. J Clin Invest (1998) 101:1870–5. doi: 10.1172/JCI1211

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Miwa I, Ichimura N, Sugiura M, Hamada Y, Taniguchi S. Inhibition of glucose-induced insulin secretion by 4-hydroxy-2-nonenal and other1529 lipid peroxidation products. Endocrinology (2000) 141:2767–2772. doi: 10.1210/endo.141.8.7614

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Liu YQ, Moibi JA, Leahy JL. Chronic 1466 high glucose lowers pyruvate dehydrogenase activity in islets through enhanced production of long chain acyl-CoA: prevention of impaired glucose oxidation by enhanced pyruvate recycling through the malate pyruvate shuttle. J Biol Chem (2004) 279:7470–5. doi: 10.1074/jbc.M307921200

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Johnson D, Shepherd RM, Gill D, Gorman T, Smith DM, Dunne MJ. Glucose dependent modulation of insulin secretion intracellular calcium ions by GKA50, a glucokinase activator. Diabetes (2007) 56:1694–702. doi: 10.2337/db07-0026

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Xu J, Han J, Long YS, Epstein PN, Liu YQ. The role of pyruvate carboxylase in insulin secretion and proliferation in rat pancreatic beta cells. Diabetologia (2008) 51:2022–30. doi: 10.1007/s00125-008-1130-9

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Xu J, Han J, Long YS, Lock J, Weir GC, Epstein PN, et al. Malic enzyme is present in mouse islets and modulates insulin secretion. Diabetologia (2008) 51:2281–9. doi: 10.1007/s00125-008-1155-0

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Lamontagne J, Pepin E, Peyot M-L, Joly E, Ruderman NB, Poitout V, et al. Pioglitazone acutely reduces insulin secretion and causes metabolic deceleration of the pancreatic beta-cell at submaximal glucose concentrations. Endocrinology (2009) 150(8):3465–74. doi: 10.1210/en.2008-1557

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Malmgren S, Spégel P, Danielsson APH, Nagorny CL, Andersson LE, Nitert MD, et al. Coordinate changes in histone modifications, mRNA levels, and metabolite profiles in clonal INS-11501 832/13 β-cells accompany functional adaptations to lipotoxicity. J Biol Chem (2013) 288:11973–87. doi: 10.1074/jbc.M112.422527

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Alcazar O, Buchwald P. Concentration-dependency and time profile of insulin secretion: dynamic perifusion studies with human and murine islets. Front Endocrinol (2019) 10:680. doi: 10.3389/fendo.2019.00680

CrossRef Full Text | Google Scholar

65. Malinowski RM, Ghiasi SM, Mandrup-Poulsen T, Meier S, Lerche MH, Ardenkjær-Larsen JH, et al. Pancreatic β-cells respond to fuel pressure with an early metabolic switch. Sci Rep (2020) 10:15413. doi: 10.1038/1496s41598-020-72348-1

PubMed Abstract | CrossRef Full Text | Google Scholar

66. König M, Dräger A, Holzhütter H-G. Cysbml: a cytoscape plugin for SBML. Bioinformatics (2012) 28:2402–3. doi: 10.1093/bioinformatics/bts432

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Sauro HM. Systems biology: introduction to pathway modeling. Ambrosius Publishing (2020).

Google Scholar

68. Brim T, Roche E, Assimacopoulos-Jeannet F, Corkey BE, Kim K-H, Prentki M. Evidence for an Anaplerotic/Malonyl-CoA pathway in pancreatic p-cell nutrient signaling. Diabetes (1996) 45:9. doi: 10.2337/diab.45.2.190

CrossRef Full Text | Google Scholar

69. Itoh Y, Kawamata Y, Harada M, Kobayashi M, Fujii R, Fukusumi S, et al. Free fatty acids regulate insulin secretion from pancreatic beta cells through GPR40. Nature (2003) 422:173–6. doi: 10.1038/nature01478

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Li C, Najafi H, Daikhin Y, Nissim IB, Collins HW, Yudkoff M, et al. Regulation of leucine-stimulated insulin secretion and glutamine metabolism in isolated rat islets. J Biol Chem (2003) 278:2853–8. doi: 10.1074/jbc.M210577200

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Roduit R, Nolan C, Alarcon C, Moore P, Barbeau A, Delghingaro-Augusto V, et al. A role for the malonyl-CoA/long-chain acyl-CoA pathway of lipid signaling in the regulation of insulin secretion in response to both fuel and nonfuel stimuli. Diabetes (2004) 53:1007–19. doi: 10.2337/diabetes.53.4.1007

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Bertuzzi A, Salinari S, Mingrone G. Insulin granule trafficking in beta-cells: mathematical model of glucose-induced insulin secretion. American journal of physiology. Endocrinol Metab (2007) 293:E396–409. doi: 10.1152/ajpendo.00647.2006

CrossRef Full Text | Google Scholar

73. Félix-Martínez GJ, Godínez-Fernández JR. Mathematical models of electrical activity of the pancreatic β-cell: a physiological review. Islets (2014) 6:e949195. doi: 10.4161/19382014.2014.949195

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Nesher R, Cerasi E. Modeling phasic insulin release: immediate and time-dependent effects of glucose. Diabetes (2002) 51(Suppl 1):S53–59. doi: 10.2337/diabetes.51.2007.s53

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Jiang N, Cox RD, Hancock JM. A kinetic core model of the glucose-stimulated insulin secretion network of pancreatic β cells. Mamm Genome (2007) 18:508–20. doi: 10.1007/s00335-007-9011-y

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Pedersen MG, Bertram R, Sherman A. Intra- and inter-islet synchronization of metabolically driven insulin secretion. Biophys J (2005) 89:107–19. doi: 10.1529/biophysj.104.055681

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Watts M, Ha J, Kimchi O, Sherman A. Paracrine regulation of glucagon secretion: the β/α/δ model. American journal of physiology. Endocrinol Metab (2016) 310:E597–611. doi: 10.1152/ajpendo.00415.2015

CrossRef Full Text | Google Scholar

78. Briant LJB, Reinbothe TM, Spiliotis I, Miranda C, Rodriguez B, Rorsman P. δ-cells and β-cells are electrically coupled and regulate α-cell activity via somatostatin. J Physiol (2018) 596:197–215. doi: 10.1113/JP274581

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Muoio DM, Newgard CB. Mechanisms of disease:Molecular and metabolic mechanisms of insulin resistance and beta-cell failure in type 2 diabetes. nature reviews. Mol Cell Biol (2008) 9:193–205. doi: 10.1038/nrm2327

CrossRef Full Text | Google Scholar

80. Holst JJ, Gasbjerg LS, Rosenkilde MM. The role of incretins on insulin function and glucose homeostasis. Endocrinology (2021) 162:bqab065. doi: 10.1210/endocr/bqab065

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Dalla Man C, Rizza RA, Cobelli C. Meal simulation model of the glucose-insulin system. IEEE Trans bio-medical Eng (2007) 54:1740–9. doi: 10.1109/TBME.2007.893506

CrossRef Full Text | Google Scholar

82. Piccinini F, Dalla Man C, Vella A, Cobelli C. A model for the estimation of hepatic insulin extraction after a meal. IEEE Trans bio-medical Eng (2016) 63:1925–32. doi: 10.1109/TBME.2015.2505507

CrossRef Full Text | Google Scholar

83. Girard J. The incretins: from the concept to their use in the treatment of type 2 diabetes. part a: incretins: concept and physiological functions. Diabetes Metab (2008) 34:550–9. doi: 10.1016/j.diabet.2008.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Riz M, Braun M, Pedersen MG. Mathematical modeling of heterogeneous electrophysiological responses in human β-cells. PloS Comput Biol (2014) 10:e1003389. doi: 10.1371/journal.pcbi.1003389

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Greenbaum CJ, Prigeon RL, D’Alessio DA. Impaired beta-cell function, incretin effect, and glucagon suppression in patients with type 1 diabetes who have normal fasting glucose. Diabetes (2002) 51:951–7. doi: 10.2337/diabetes.51.4.951

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Nauck MA, Vilsbøll T, Gallwitz B, Garber A, Madsbad S. Incretin-based therapies: viewpoints on the way to consensus. Diabetes Care (2009) 32(Suppl 2):S223–231. doi: 10.2337/dc09-S315

PubMed Abstract | CrossRef Full Text | Google Scholar

87. König M, Holzhütter H-G. Kinetic modeling of human hepatic glucose metabolism in type 2 diabetes mellitus predicts higher risk of hypoglycemic events in rigorous insulin therapy. J Biol Chem (2012) 287:36978–89. doi: 10.1074/jbc.M112.382069

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Pacini G, Ahrén B, Göbl C, Tura A. Assessing the effect of incretin hormones and other insulin secretagogues on pancreatic beta-cell function: review on mathematical modelling approaches. Biomedicines (2022) 10:1060. doi: 10.3390/biomedicines10051060

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Tura A, Göbl C, Vardarli I, Pacini G, Nauck M. Insulin clearance and incretin hormones following oral and "isoglycemic" intravenous glucose in type 2 diabetes patients under different antidiabetic treatments. Sci Rep (2022) 12:2510. doi: 10.1038/s41598-022-06402-5

PubMed Abstract | CrossRef Full Text | Google Scholar

90. König M. Sbmlutils: Python utilities for SBML. (2022). doi: 10.5281/zenodo.7462781sbmlutils-0.8.1.

CrossRef Full Text | Google Scholar

91. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

92. Ebrahim A, Beber ME, Mandal S, König M, Redestig H, Diener C, et al. Opencobra/cobrapy: 0.26.2. (2023). doi: 10.5281/zenodo.7502407.

CrossRef Full Text | Google Scholar

93. Neal ML, König M, Nickerson D, Mısırlı G, Kalbasi R, Dräger A, et al. Harmonizing semantic annotations for computational models in biology. Briefings Bioinf (2019) 20:540–50. doi: 10.1093/bib/bby087

CrossRef Full Text | Google Scholar

94. Neal ML, Gennari JH, Waltemath D, Nickerson DP, König M. Open modeling and exchange (OMEX) metadata specification version 1.0. J Integr Bioinf (2020) 17:20200020. doi: 10.1515/jib-2020-0020

CrossRef Full Text | Google Scholar

95. Courtot M, Juty N, Knüpfer C, Waltemath D, Zhukova A, Dräger A, et al. Controlled vocabularies and semantics in systems biology. Mol Syst Biol (2011) 7:543. doi: 10.1038/msb.2011.77

PubMed Abstract | CrossRef Full Text | Google Scholar

96. Cote R, Reisinger F, Martens L, Barsnes H, Vizcaino JA, Hermjakob H. The ontology lookup service: bigger and better. Nucleic Acids Res (2010) 38:W155–60. doi: 10.1093/nar/gkq331

PubMed Abstract | CrossRef Full Text | Google Scholar

97. Noronha A, Modamio J, Jarosz Y, Guerard E, Sompairac N, Preciat G, et al. The virtual metabolic human database: integrating human and gut microbiome metabolism with nutrition and disease. Nucleic Acids Res (2019) 47:D614–24. doi: 10.1093/nar/gky992

PubMed Abstract | CrossRef Full Text | Google Scholar

98. King ZA, Lu J, Dräger A, Miller P, Federowicz S, Lerman JA, et al. BiGG models: a platform for integrating, standardizing and sharing genome-scale models. Nucleic Acids Res (2016) 44:D515–522. doi: 10.1093/nar/gkv1049

PubMed Abstract | CrossRef Full Text | Google Scholar

99. Hari A, Lobo D. Mergem: merging and comparing genome-scale metabolic models using universal identifiers. bioRxiv (2022). doi: 10.1101/2022.07.14.499633

PubMed Abstract | CrossRef Full Text | Google Scholar

100. The UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res (2017) 45:D158–69. doi: 10.1093/nar/gkw1099

PubMed Abstract | CrossRef Full Text | Google Scholar

101. König M. Pymetadata: Python utilities for SBML. (2022). doi: 10.5281/zenodo.7432576pymetadata-v0.3.8.

CrossRef Full Text | Google Scholar

102. Placzek S, Schomburg I, Chang A, Jeske L, Ulbrich M, Tillack J, et al. BRENDA in 2017: new perspectives and new tools in BRENDA. Nucleic Acids Res (2017) 45:D380–8. doi: 10.1093/nar/gkw952

PubMed Abstract | CrossRef Full Text | Google Scholar

103. König M. Brendapy: BRENDA parser in python. (2022). doi: 10.5281/zenodo.6555202brendapy-0.5.0.

CrossRef Full Text | Google Scholar

104. Wittig U, Rey M, Weidemann A, Kania R, Müller W. SABIO-RK: an updated resource for manually curated biochemical reaction kinetics. Nucleic Acids Res (2018) 46:D656–60. doi: 10.1093/nar/gkx1065

PubMed Abstract | CrossRef Full Text | Google Scholar

105. Grzegorzewski J, Brandhorst J, Green K, Eleftheriadou D, Duport Y, Barthorscht F, et al. PK-DB: pharmacokinetics database for individualized and stratified computational modeling. Nucleic Acids Res (2021) 49:D1358–64. doi: 10.1093/nar/gkaa990

PubMed Abstract | CrossRef Full Text | Google Scholar

106. Grzegorzewski J, Bartsch F, Köller A, König M. Pharmacokinetics of caffeine: a systematic analysis of reported data for application in metabolic phenotyping and liver function testing. Front Pharmacol (2021) 12:752826 Grzegorzewski2021a. doi: 10.3389/fphar.2021.752826

CrossRef Full Text | Google Scholar

107. Rohatgi A. Webplotdigitizer: version 4.5. (2021).

Google Scholar

108. Deepa Maheshvare M, König M. Model of glucose-stimulated insulin secretion in the pancreatic β-cell. (2023). doi: 10.5281/zenodo.7932991

PubMed Abstract | CrossRef Full Text | Google Scholar

109. Deepa Maheshvare M, Raha S, Konig M, Pal D. A consensus model of glucose-stimulated insulin secretion in the pancreatic beta-cell. bioRxiv (2023). doi: 10.1101/2023.03.10.532028

PubMed Abstract | CrossRef Full Text | Google Scholar

110. Cohn SE. An introduction to estimation theory. J Meteorological Soc Japan (1997).

Google Scholar

111. Dean PM. Ultrastructural morphometry of the pancreatic β-cell. Diabetologia (1973) 9:115–9. doi: 10.1007/BF01230690

PubMed Abstract | CrossRef Full Text | Google Scholar

112. Goldberg R, Tewari Y. Thermodynamics of enzyme-catalyzed reactions. New York: McGraw Hill (2003).

Google Scholar

113. Noor E, Haraldsdóttir HS, Milo R, Fleming RMT. Consistent estimation of Gibbs energy using component contributions. PloS Comput Biol (2013) 9:e1003098. doi: 10.1371/journal.pcbi.1003098

PubMed Abstract | CrossRef Full Text | Google Scholar

114. Liebermeister W, Uhlendorf J, Klipp E. Modular rate laws for enzymatic reactions: thermodynamics, elasticities and implementation. Bioinf (Oxford England) (2010) 26:1528–34. doi: 10.1093/bioinformatics/btq141

CrossRef Full Text | Google Scholar

115. De Gaetano A, Gaz C, Palumbo P, Panunzi S. A unifying organ model of pancreatic insulin secretion. PloS One (2015) 10:e0142344. doi: 10.1371/journal.pone.0142344

PubMed Abstract | CrossRef Full Text | Google Scholar

116. Bergmann FT. Basico: a simplified python interface to COPASI. (2023). doi: 10.5281/zenodo.7541481.

CrossRef Full Text | Google Scholar

117. Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, et al. COPASI–a complex pathway simulator. Bioinformatics (2006) 22:3067–74. doi: 10.1093/bioinformatics/btl485

PubMed Abstract | CrossRef Full Text | Google Scholar

118. Fröhlich F, Weindl D, Schälte Y, Pathirana D, Paszkowski Ł., Lines GT, et al. AMICI: high-performance sensitivity analysis for large ordinary differential equation models. Bioinf (Oxford England) (2021) 37:3676–7. doi: 10.1093/bioinformatics/btab227

CrossRef Full Text | Google Scholar

119. Schmiester L, Schälte Y, Bergmann FT, Camba T, Dudkin E, Egert J, et al. Petab–interoperable specification of parameter estimation problems in systems biology. PloS Comput Biol (2021) 17:1–10. doi: 10.1371/journal.pcbi.1008646

CrossRef Full Text | Google Scholar

120. Kent E, Hoops S, Mendes P. Condor-COPASI: high-throughput computing for biochemical networks. BMC Syst Biol (2012) 6:91. doi: 10.1186/1752-0509-6-91

PubMed Abstract | CrossRef Full Text | Google Scholar

121. Rodriguez-Fernandez M, Mendes P, Banga JR. A hybrid approach for efficient and robust parameter estimation in biochemical pathways. Bio Syst (2006) 83:248–65. doi: 10.1016/j.biosystems.2005.06.016

CrossRef Full Text | Google Scholar

122. Angeline PJ. Evolutionary optimization versus particle swarm optimization: philosophy and performance differences. (2005), 601–10.

Google Scholar

123. Dana S, Nakakuki T, Hatakeyama M, Kimura S, Raha S. Computation of restoration of ligand response in the random kinetics of a prostate cancer cell signaling pathway. Comput Methods Programs Biomed (2011) 101:1–22. doi: 10.1016/j.cmpb.2010.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

124. El Rassy E, Delaroque A, Sambou P, Chakravarty HK, Matynia A. On the potential of the particle swarm algorithm for the optimization of detailed kinetic mechanisms. comparison with the genetic algorithm. J Phys Chem (2021) 125:5180–9. doi: 10.1021/acs.jpca.1c02095

CrossRef Full Text | Google Scholar

125. König M. Sbmlsim: SBML simulation made easy. (2021). doi: 10.5281/zenodo.5531088sbmlsim-0.2.2.

CrossRef Full Text | Google Scholar

126. Somogyi ET, Bouteiller J-M, Glazier JA, König M, Medley JK, Swat MH, et al. libRoadRunner: a high performance SBML simulation and analysis library. Bioinf (Oxford England) (2015) 31:3315–21. doi: 10.1093/bioinformatics/btv363

CrossRef Full Text | Google Scholar

127. Welsh C, Xu J, Smith L, König M, Choi K, Sauro HM. libRoadRunner 2.0: a high performance SBML simulation and analysis library. Bioinf (Oxford England) (2023) 39:btac770. doi: 10.1093/bioinformatics/btac770

CrossRef Full Text | Google Scholar

128. Hucka M, Bergmann FT, Chaouiya C, Dräger A, Hoops S, Keating SM, et al. The systems biology markup language (SBML): language specification for level 3 version 2 core release 2. J Integr Bioinf (2019) 16:20190021. doi: 10.1515/jib-2019-0021

CrossRef Full Text | Google Scholar