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.

A Consensus 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 the first consensus 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 equations for insulin secretion coupled to cellular energy state (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 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.

A Graph-Based Framework for Multiscale Modeling of Physiological Transport.

Trillions of chemical reactions occur in the human body every second, where the generated products are not only consumed locally but also transported to various locations in a systematic manner to sustain homeostasis. Current solutions to model these biological phenomena are restricted in computability and scalability due to the use of continuum approaches in which it is practically impossible to encapsulate the complexity of the physiological processes occurring at diverse scales. Here, we present a discrete modeling framework defined on an interacting graph that offers the flexibility to model multiscale systems by translating the physical space into a metamodel. We discretize the graph-based metamodel into functional units composed of well-mixed volumes with vascular and cellular subdomains; the operators defined over these volumes define the transport dynamics. We predict glucose drift governed by advective-dispersive transport in the vascular subdomains of an islet vasculature and cross-validate the flow and concentration fields with finite-element-based COMSOL simulations. Vascular and cellular subdomains are coupled to model the nutrient exchange occurring in response to the gradient arising out of reaction and perfusion dynamics. The application of our framework for modeling biologically relevant test systems shows how our approach can assimilate both multi-omics data from in vitro-in vivo studies and vascular topology from imaging studies for examining the structure-function relationship of complex vasculatures. The framework can advance simulation of whole-body networks at user-defined levels and is expected to find major use in personalized medicine and drug discovery.

Partitioning a reaction-diffusion ecological network for dynamic stability.

The loss of dispersal connections between habitat patches may destabilize populations in a patched ecological network. This work studies the stability of populations when one or more communication links is removed. An example is finding the alignment of a highway through a patched forest containing a network of metapopulations in the patches. This problem is modelled as that of finding a stable cut of the graph induced by the metapopulations network, where nodes represent the habitat patches and the weighted edges model the dispersal between habitat patches. A reaction-diffusion system on the graph models the dynamics of the predator-prey system over the patched ecological network. The graph Laplacian's Fiedler value, which indicates the well-connectedness of the graph, is shown to affect the stability of the metapopulations. We show that, when the Fiedler value is sufficiently large, the removal of edges without destabilizing the dynamics of the network is possible. We give an exhaustive graph partitioning procedure, which is suitable for smaller networks and uses the criterion for both the local and global stability of populations in partitioned networks. A heuristic graph bisection algorithm that preserves the preassigned lower bound for the Fiedler value is proposed for larger networks and is illustrated with examples.

Adaptive Control Based Harvesting Strategy for a Predator-Prey Dynamical System.

This paper deals with designing a harvesting control strategy for a predator-prey dynamical system, with parametric uncertainties and exogenous disturbances. A feedback control law for the harvesting rate of the predator is formulated such that the population dynamics is asymptotically stabilized at a positive operating point, while maintaining a positive, steady state harvesting rate. The hierarchical block strict feedback structure of the dynamics is exploited in designing a backstepping control law, based on Lyapunov theory. In order to account for unknown parameters, an adaptive control strategy has been proposed in which the control law depends on an adaptive variable which tracks the unknown parameter. Further, a switching component has been incorporated to robustify the control performance against bounded disturbances. Proofs have been provided to show that the proposed adaptive control strategy ensures asymptotic stability of the dynamics at a desired operating point, as well as exact parameter learning in the disturbance-free case and learning with bounded error in the disturbance prone case. The dynamics, with uncertainty in the death rate of the predator, subjected to a bounded disturbance has been simulated with the proposed control strategy.

Numerical simulation of a glucose sensitive composite membrane closed-loop insulin delivery system.

Closed-loop insulin delivery system works on pH modulation by gluconic acid production from glucose, which in turn allows regulation of insulin release across membrane. Typically, the concentration variation of gluconic acid can be numerically modeled by a set of non-linear, non-steady state reaction diffusion equations. Here, we report a simpler numerical approach to time and position dependent diffusivity of species using finite difference and differential quadrature (DQ) method. The results are comparable to that obtained by analytical method. The membrane thickness directly determines the concentrations of the glucose and oxygen in the system, and inversely to the gluconic acid. The advantage with the DQ method is that its parameter values need not be altered throughout the analysis to obtain the concentration profiles of the glucose, oxygen and gluconic acid. Our work would be useful for modeling diabetes and other systems governed by such non-linear and non-steady state reaction diffusion equations.

Time scales of the stick-slip dynamics of the peeling of an adhesive tape.

The stick-slip dynamics of the peeling of an adhesive tape is characterized by bifurcations that have been experimentally well studied. In this work, we investigate the time scale in which the the stick-slips happen leading to the bifurcations. This is fundamental to understanding the triboluminescence and acoustic emissions associated with the bifurcations. We establish a relationship between the time scale of the bifurcations and the inherent mathematical structure of the peeling dynamics by studying a characteristic time quantity associated with the dynamics.

Computation of restoration of ligand response in the random kinetics of a prostate cancer cell signaling pathway.

Mutation and/or dysfunction of signaling proteins in the mitogen activated protein kinase (MAPK) signal transduction pathway are frequently observed in various kinds of human cancer. Consistent with this fact, in the present study, we experimentally observe that the epidermal growth factor (EGF) induced activation profile of MAP kinase signaling is not straightforward dose-dependent in the PC3 prostate cancer cells. To find out what parameters and reactions in the pathway are involved in this departure from the normal dose-dependency, a model-based pathway analysis is performed. The pathway is mathematically modeled with 28 rate equations yielding those many ordinary differential equations (ODE) with kinetic rate constants that have been reported to take random values in the existing literature. This has led to us treating the ODE model of the pathways kinetics as a random differential equations (RDE) system in which the parameters are random variables. We show that our RDE model captures the uncertainty in the kinetic rate constants as seen in the behavior of the experimental data and more importantly, upon simulation, exhibits the abnormal EGF dose-dependency of the activation profile of MAP kinase signaling in PC3 prostate cancer cells. The most likely set of values of the kinetic rate constants obtained from fitting the RDE model into the experimental data is then used in a direct transcription based dynamic optimization method for computing the changes needed in these kinetic rate constant values for the restoration of the normal EGF dose response. The last computation identifies the parameters, i.e., the kinetic rate constants in the RDE model, that are the most sensitive to the change in the EGF dose response behavior in the PC3 prostate cancer cells. The reactions in which these most sensitive parameters participate emerge as candidate drug targets on the signaling pathway.

Activator-induced dynamic disorder and molecular memory in human two-pore domain hTREK1 K channel.

UNLABELLED: Ion channels are fundamental molecules in the nervous system that catalyze the flux of ions across the cell membrane. Ion channel flux activity is comparable to the catalytic activity of enzyme molecules. Saturating concentrations of substrate induce "dynamic disorder" in the kinetic rate processes of single-enzyme molecules and consequently, develop correlative "memory" of the previous history of activities. Similarly, binding of ions as substrate alone or in presence of agonists affects the catalytic turnover of single-ion channels. Here, we investigated the possible existence of dynamic disorder and molecular memory in the single human-TREK1-channel due to binding of substrate/agonist using the excised inside-out patch-clamp technique. Our results suggest that the single-hTREK1-channel behaves as a typical Michaelis-Menten enzyme molecule with a high-affinity binding site for K(+) ion as substrate. But, in contrast to enzyme, dynamic disorder in single-hTREK1-channel was not induced by substrate K(+) binding, but required allosteric modification of the channel molecule by the agonist, trichloroethanol. In addition, interaction of trichloroethanol with hTREK1 induced strong correlation in the waiting time and flux intensity, exemplified by distinct mode-switching between high and low flux activities. This suggested the induction of molecular memory in the channel molecule by the agonist, which persisted for several decades in time. Our mathematical modeling studies identified the kinetic rate processes associated with dynamic disorder. It further revealed the presence of multiple populations of distinct conformations that contributed to the "heterogeneity" and consequently, to the molecular memory phenomenon that we observed. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (doi:10.1007/s12154-010-0053-3) contains supplementary material, which is available to authorized users.

A combination of docking, QM/MM methods, and MD simulation for binding affinity estimation of metalloprotein ligands.

To alleviate the problems in the receptor-based design of metalloprotein ligands due to inadequacies in the force-field description of coordination bonds, a four-tier approach was devised. Representative ligand-metalloprotein interaction energies are obtained by subsequent application of (1) docking with metal-binding-guided selection of modes, (2) optimization of the ligand-metalloprotein complex geometry by combined quantum mechanics and molecular mechanics (QM/MM) methods, (3) conformational sampling of the complex with constrained metal bonds by force-field-based molecular dynamics (MD), and (4) a single point QM/MM energy calculation for the time-averaged structures. The QM/MM interaction energies are, in a linear combination with the desolvation-characterizing changes in the solvent-accessible surface areas, correlated with experimental data. The approach was applied to structural correlation of published binding free energies of a diverse set of 28 hydroxamate inhibitors to zinc-dependent matrix metalloproteinase 9 (MMP-9). Inclusion of steps 3 and 4 significantly improved both correlation and prediction. The two descriptors explained 90% of variance in inhibition constants of all 28 inhibitors, ranging from 0.08 to 349 nM, with the average unassigned error of 0.318 log units. The structural and energetic information obtained from the time-averaged MD simulation results helped understand the differences in binding modes of related compounds.

A comparison of the binding sites of matrix metalloproteinases and tumor necrosis factor-alpha converting enzyme: implications for selectivity.

MMPs and TACE (ADAM-17) assume independent, parallel, or opposite pathological roles in cancer, arthritis, and several other diseases. For therapeutic purposes, selective inhibition of individual MMPs and TACE is required in most cases due to distinct roles in diseases and the need to preserve activities in normal states. Toward this goal, we compared force-field interaction energies of five ubiquitous inhibitor atoms with flexible binding sites of 24 known human MMPs and TACE. The results indicate that MMPs 1-3, 10, 11, 13, 16, and 17 have at least one subsite very similar to TACE. S3 subsite is the best target for development of specific TACE inhibitors. Specific binding to TACE compared to most MMPs is promoted by placing a negatively charged ligand part at the bottom of S2 subsite, at the entrance of S1' subsite, or the part of S3' subsite that is close to catalytic zinc. Numerous other clues, consistent with available experimental data, are provided for design of selective inhibitors.

Processing multimode binding situations in simulation-based prediction of ligand-macromolecule affinities.

The linear response (LR) approximation and similar approaches belong to practical methods for estimation of ligand-receptor binding affinities. The approaches correlate experimental binding affinities with the changes upon binding of the ligand electrostatic and van der Waals energies and of solvation characteristics. These attributes are expressed as ensemble averages that are obtained by conformational sampling of the protein-ligand complex and of the free ligand by molecular dynamics or Monte Carlo simulations. We observed that outliers in the LR correlations occasionally exhibit major conformational changes of the complex during sampling. We treated the situation as a multimode binding case, for which the observed association constant is the sum of the partial association constants of individual states/modes. The resulting nonlinear expression for the binding affinities contains all the LR variables for individual modes that are scaled by the same two to four adjustable parameters as in the one-mode LR equation. The multimode method was applied to inhibitors of a matrix metalloproteinase, where this treatment improved the explained variance in experimental activity from 75% for the unimode case to about 85%. The predictive ability scaled accordingly, as verified by extensive cross-validations.

Similarity of binding sites of human matrix metalloproteinases.

Tissue components hydrolyzing matrix metalloproteinases (MMPs) exhibit a high sequence similarity (56-64% in catalytic domains) and yet a significant degree of functional specificity. The hexapeptide-binding sites of 24 known human MMPs were compared in terms of their force field interaction energies with five probes that are most frequently encountered in substrates and inhibitors. The probes moved along a grid enclosing partially flexible binding sites in rigid catalytic domains that were represented by published experimental structures and comparative models and new comparative models for nine most recently characterized MMPs. For individual MMPs, representative interaction energies were obtained as averages for all suitable experimental structures. Correlations of the representative energies for all MMP pairs were succinctly catalogued for individual probes, subsites, and correlation levels. Among the probes (neutral sp(3) carbon and sp(3) oxygen, positive sp(3) nitrogen and hydrogen, and negative carbonyl oxygen), the last probe is least distinctive. Similarities of subsites are decreasing as S1 ' > S2 > S3 ' > S1 approximately S3 > S2 '. Most interesting, occupancies of subsites in published structures of MMP-inhibitor complexes follow an almost parallel trend, alluding to overall low selectivity of known MMP inhibitors. Flexible subsite S1 ' that appears as the specificity pocket in rigid x-ray structures is actually very similar among individual MMPs. Several correlations indicated that MMPs 3, 8, and 12 have similar binding sites. Modeling results are corroborated with published experimental data on MMP inhibition and substrate specificities. The results provide numerous clues for development of specific inhibitors and substrates, as well as for selection of MMPs for testing that provides maximum information without redundant experiments.