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.


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 b-cells, and glucagon, secreted by the a-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 b-cells in the pancreas through glucose transporters, primarily GLUT2 (3). Once inside the b-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 b-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 2+ 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 b-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)(13)(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 2+ response and insulin release. This trend is consistent across several studies (15)(16)(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 b-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,6bisphosphate) 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, b-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 b-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 b-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.

Results
Our study introduces a detailed kinetic model of GSIS in the pancreatic b-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 b-cell on biphasic insulin secretion.

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 bcell 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 timecourse 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.
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 bcell 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 b-cell glycolysis that is easily accessible to the scientific community for further use.

Reproducible modeling workflow
In this study, we describe a comprehensive modeling workflow for building small kinetic pathway models ( Figure 2) using SBML (32, 33).
In our model-building workflow, we followed several steps to construct a novel kinetic SBML model of glycolysis in the pancreatic b-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 b-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 (b-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 timecourse 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.

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 b-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 b-cell to insulin secretion.
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.

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 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 Frontiers in Endocrinology frontiersin.org 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 leastsquares 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.

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 glucosedependent 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 b-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).
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. 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 b-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.

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). Computational model of glucose-stimulated insulin secretion (GSIS) in the pancreatic b-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 b-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 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 b-cell.

Discussion
We have developed a comprehensive kinetic model of GSIS in the pancreatic b-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 b-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 bcells. Precisely, high ATP levels prompt the closure of ATP-sensitive potassium channels (6). Consequently, the cell membrane depolarizes, opening voltagegated 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). 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 b-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 b-cell. While the existing   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-3phosphate (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-3phosphate dehydrogenase and phosphoglycerate kinase could potentially be affected due to the lack of steady-state and timecourse 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 preincubation 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)(13)(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, nonstandardized 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 b-cell glycolysis, a data curation workflow, strategies for unit and data normalization, and a large database of metabolic data of the pancreatic b-cell. Our model can be extended further to study glucose-insulin regulatory mechanisms by for example integrating calcium handling in b-cell (28, 29,76) or paracrine regulation of the release of counterregulatory hormones such as glucagon by gap-junctional coupling of a-, b-, and d-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 b-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 b-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 b-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 glucosedependent 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)(86)(87)(88)(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 b-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.

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

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 b-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).

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).
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).
The formula and charge of bpg13 are C3H4O10P2 and -4, respectively.

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.

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).

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.

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 b-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", "glucosedependence", "pancreatic b-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 metaanalysis (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).

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 (b-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.

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 twostep 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).

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 a, (c a 0 , c a 1 , …, c a n ) observed at n incubation glucose doses (d a 0 , d a 1 , …, d a n ) is expressed by the piecewise linear-interpolation function C ss . Here, a belongs to the set of steady-state experiments 1 ≤ a ≤ N a ss with N a ss being the number of steady-state experimental curves of the metabolite s. Mean curve. The mean steady-state curve 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 C ss . The cost function F of the minimization problem is given by, In Eq. 1, C ss (d a i )and C ss (d a i ) 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 a.
For each experimental curve, the factor f a 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 a of the experimental curve a (Eq. 2).
The scale factors of all steady state curves (f 1 , …f N ass ) were determined by minimizing the respective cost functions (F(f 1 ), … F(f N ass ). Multiplying the experimental curve C a by the scaling factor f a 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 a 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. For a given incubation glucose dose d b , the metabolite concentration at the last time point C tc (t m ) corresponds to the steady state value reached for the given d b :

Time course data normalization
The scaling factor for the time course experiment follows as:

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 b-cell volume were obtained from the values reported in a morphometric study of b-cells (111).
Initial concentrations. The initial concentrations of glycolytic intermediates and adenine nucleotides were obtained from the mean curve 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.
Here, a i is S i =Km s , b i is P i =Km p , S refers to the substrate and P refers to the product. K eq is the equilibrium constant and G 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 b-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 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 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 IRS sp max 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.
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 glucosedependent 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 = ATP tot ) changes with glucose. To account for these changes, a reaction DATP was added that changes the total ATP according to the observed steady-state data for a given glucose value (Eq. 7) .
The ATP tot (glc) values are determined by the interpolating polynomial of the mean steady-state glucose dose response of the ATP+ADP data.

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 In Eq. 8, c a s is the concentration of the metabolite s in the experiment a and c M s is the concentration of the metabolite s predicted by the model M. 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 preincubation 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. Presimulation 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 (glc 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 frontend 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)(122)(123)(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.

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). Presimulations 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 preincubation 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.