Date Published: February 14, 2018
Publisher: Public Library of Science
Author(s): Joydeep Sarkar, Gaurav Dwivedi, Qian Chen, Iris E. Sheu, Mark Paich, Colleen M. Chelini, Paul M. D’Alessandro, Samuel P. Burns, Joseph DiStefano.
A computational model of the physiological mechanisms driving an individual’s health towards onset of type 2 diabetes (T2D) is described, calibrated and validated using data from the Diabetes Prevention Program (DPP). The objective of this model is to quantify the factors that can be used for prevention of T2D. The model is energy and mass balanced and continuously simulates trajectories of variables including body weight components, fasting plasma glucose, insulin, and glycosylated hemoglobin among others on the time-scale of years. Modeled mechanisms include dynamic representations of intracellular insulin resistance, pancreatic beta-cell insulin production, oxidation of macronutrients, ketogenesis, effects of inflammation and reactive oxygen species, and conversion between stored and activated metabolic species, with body-weight connected to mass and energy balance. The model was calibrated to 331 placebo and 315 lifestyle-intervention DPP subjects, and one year forecasts of all individuals were generated. Predicted population mean errors were less than or of the same magnitude as clinical measurement error; mean forecast errors for weight and HbA1c were ~5%, supporting predictive capabilities of the model. Validation of lifestyle-intervention prediction is demonstrated by synthetically imposing diet and physical activity changes on DPP placebo subjects. Using subject level parameters, comparisons were made between exogenous and endogenous characteristics of subjects who progressed toward T2D (HbA1c > 6.5) over the course of the DPP study to those who did not. The comparison revealed significant differences in diets and pancreatic sensitivity to hyperglycemia but not in propensity to develop insulin resistance. A computational experiment was performed to explore relative contributions of exogenous versus endogenous factors between these groups. Translational uses to applications in public health and personalized healthcare are discussed.
Managing the care of people with diabetes is an enormous burden on the US healthcare economy, costing $176 billion per year in direct medical expenses according to estimates from 2012 . The Centers for Disease Control and Prevention (CDC) estimate that in the US diabetes affects nearly 29 million individuals (~9% of US population) . Diabetes significantly increases the risk of several co-morbidities including heart attack and stroke, retinopathy, risk of amputation, nephropathy, neuropathy, hearing loss, and depression among others [3–9]. In addition, an estimated 86 million people have prediabetes, which puts individuals at increased risk of type 2 diabetes (T2D), heart disease, and stroke . While the prevalence of diabetes is stabilizing in the US, the global cost of diabetes is expected to increase by more than 50% by the year 2035 with developing countries being the largest contributors . New approaches are required to curb the increasing global prevalence of T2D and improve care for patients with diabetes to control disease progression.
The transport of glucose into the muscle component is driven by glucose gradient, and the permeability of glucose is defined as a non-linear function of the concentrations of GLUT1 and GLUT4 transporters on the cell surface (v1 in S1 Fig and S3 Table). The GLUT1/GLUT4 cell surface concentrations are inputs from the insulin resistance (ISR) component. As described in the ISR component, GLUT1 acts as the constitutive transporter of glucose. Its trafficking is relatively weakly regulated by insulin and it acts as the insulin-independent glucose transporter. The cellular trafficking of GLUT4 is strongly insulin dependent and it acts as the insulin-dependent glucose transporter of glucose. Glucose transport into cells is assumed to increase linearly with GLUT1 and in a sigmoidal manner with GLUT4, based on previously reported sigmoidal response of GLUT4 to increasing insulin concentration . As described in the section “General mathematical model structure”, the simplest possible model was chosen for GLUT1, and a more complex model was used for GLUT4 in view of the supporting experimental data.
While the liver is the main site for de novo lipogenesis (DNL), in conditions of excess glucose inside the muscle component, there is conversion to fat in this component as well. The conversion of glucose to fat is modeled as a Hill function such that significant lipogenesis only happens in conditions of significant excess of glucose over steady-state conditions (v8 in S2 Fig and S3 Table).
The model represents the dynamic equilibrium between glucose and glycogen in the muscle component. In the model, glycogen is made from 4 glucose molecules and this stoichiometry is maintained throughout the model. Any excess glucose is converted to glycogen until a maximum glycogen threshold is reached. As the simplest possible model, the rate of glycogenesis was assumed proportional to the concentration of glucose in the tissue, insulin sensitivity (IS, an input from the ISR component) and the remaining availability for glycogen storage. The glycogenolysis process is proportional to glycogen concentration but is further controlled by the energy state of the muscle cells measured as the ATP/ADP ratio (AAR; terms v7 and v10 in S2 Fig and S3 Table).
Serum triglycerides (TGs) and chylomicrons serve as the primary sources of fat to the muscle component. Lipoprotein lipase (not explicitly modeled) hydrolyses circulating TGs and chylomicrons (v13 in S2 Fig and S3 Table) into free fatty acids (FFAs) and glycerol in the muscle component. Glycerol transport back into blood is represented as a first order mass action (v4). There is bi-directional transport of FFAs between blood and muscle component (v2, v3).
The amino acid fluxes in and out of the muscle component are represented as first order mass action processes (v5, v6). Inside the muscle component amino acid and protein concentrations are in dynamic equilibrium with a stoichiometric ratio of 500 (v16, v17). Amino acids are converted into ketoacids (v18) for use in energy metabolism (S2 Fig and S3 Table). The contribution of amino acids to energy production is smaller than that of carbohydrates and fats. The main purpose served by amino acids is to maintain the body’s nitrogen balance. Since we have not simulated starvation, the amino acid metabolism component is designed to maintain a steady nitrogen balance and excess protein is excreted .
The muscle component is the main site of macronutrient oxidation and ATP hydrolysis. Glucose, fat, and ketoacids fuel the ATP synthesis process in the model. The model only represents aerobic oxidation of the fuels as a single step reaction. The rate of ATP generation from all three types of fuels is proportional to the concentration of ADP and mitochondrial function (input from ISR component). The rate of oxidation of ketoacids is directly proportional to the concentration of ketoacids in the components while the rates of oxidation of glucose and FFA are represented as saturating processes (terms v19, v9 and v13 in S2 Fig and S3 Table). Since systemic protein levels are maintained at an approximately steady level in the model (as described in “Amino acid update and metabolism”), a simple mass action term was sufficient to capture their nearly constant contribution to energy production under non-starvation conditions. Glucose and FFA are the main contributors to ATP generation regulated by a series of enzymes. The complex process of oxidative phosphorylation was simplified to a single mathematical expression representing the saturable nature of enzyme catalyzed processes. The three fuels lead to the generation of different amounts of ATP and this stoichiometry is also incorporated into the equations. The stoichiometry is tuned such that the respiratory quotient is 0.85 at steady state and ketoacids contribute approximately 20% to ATP synthesis.
The mechanism of glucose uptake by the liver component (v1, v2 in S3 Fig and S4 Table) follows a similar mechanism as that in the muscle component. However, unlike the muscle where there is no net outflow of glucose from the muscle component into blood, there is secretion of glucose from the liver into blood. Hence, the influx and efflux terms are separated. The influx of glucose is represented as a function of GLUT1/GLUT4 and serum glucose. The efflux term (v2 in S3 Fig and v2LVR in S4 Table) is separately represented and depends on liver glucose concentration.
The liver is the primary site for de novo lipogensis (DNL) in the body. The expression for DNL (v12, S3 Fig and S4 Table) in liver is the same as that in the muscle component.
Glycogen-glucose dynamics in the liver component are represented the same way as in the muscle component. During rigorous physical activity liver glycogen is preferentially used as a source of glucose . It has been reported that the suspected mechanism for this enhanced glycogenolysis is increase in levels of catecholamines like epinephrine [45,46]. Hence, the form of the function for change in glycogenolysis follows data for levels of epinephrine with increasing intensity of exercise .
The liver component is the site of gluconeogenesis. The model incorporates gluconeogenesis from glycerol as well as from a glucogenic portion of amino acids. Gluconeogenesis from ketoacids is assumed to be proportional to the concentration of ketoacids in the liver (v17, v21 in S3 Fig and S4 Table). Gluconeogenesis from glycerol is inhibited by increased insulin sensitivity (input from the ISR component).
Unlike in the muscle, in the liver, there is no uptake of FFA from TGs and chylomicrons. There is uptake of FFAs from blood into liver that follows a 1st order process (v3 in S3 Fig and S4 Table). New FFAs that are synthesized through DNL and FFA absorbed from blood are combined with glycerol to form TG that is transported out (v14, v4). Lipolysis of TG in the liver is ignored.
The model also incorporates a simplified representation of the ketogenesis (v15 in S3 Fig and S4 Table) that increases significantly when there is excess FFA in the liver component. Some ketone bodies are reconverted to FFA (v16) while the rest are transported to blood over a gradient (v7, v8).
The adipose component represents the site of both visceral and subcutaneous fat tissue. The adipose component incorporates uptake of fat from blood, and the dynamic equilibrium between lipolysis and esterification of FFA. The fatty acid metabolism in adipose tissue is modeled similarly to the muscle component (Fig 1, S4 Fig), with difference being that adipose tissue takes up more FFA from circulating TG and chylomicrons than muscle. The increase in lipolysis by physical activity [48,49] is captured in the model. It has also been shown that the increase of lipolysis saturates with increased intensity of physical activity ; therefore, a saturating function is used to connect lipolysis with increasing intensity of physical activity. Finally, regulation of lipolysis by insulin is also included in the model such that increased insulin sensitivity (input from ISR component) decreases the rate of lipolysis  (v4 in S4 Fig and S5 Table).
The proposed pancreas model is a modification of that presented by De Gaetano et al. . The model incorporates the long-term dynamics of beta cell mass, beta cell function and insulin production as distinct elements (Fig 1 and S5 Fig). While the fasting plasma insulin is predominantly dictated by the dynamics of the beta cell number, the increased response of insulin secretion to increased serum glucose is predominantly driven by beta cell function. Proliferation of beta cells is represented as a Hill function dependent on blood glucose concentration (v1 in S6 Table). The choice of the Hill function implies that beta cell growth, and hence subsequent insulin production, is low for very low glucose concentrations and increases rapidly toward a saturation point when glucose concentration crosses a certain threshold. Beta cell apoptosis is a function of chronic inflammation (v2). Beta cell function, a value bounded between 0 and 1, decreases at a rate proportional to remaining beta cell function and ROS concentrations. The beta cell function recovers by a first order mass action process. Beta cell function was introduced to take into account the lack of insulin producing capacity in the remaining beta cells in diabetic subjects. The rate of insulin production is proportional to the beta cell mass, beta cell functional capacity, and also increases in a saturating manner (through independent Hill functions) with serum glucose and FFA concentrations. Insulin is removed from blood by a simple first order mass action process.
The insulin resistance (ISR) component (Fig 1, S6 Fig and S7 Table), represents processes that modulate the following: a) the response of cells to insulin, simulating deficiency in insulin action as insulin resistance increases; b) short and long-term regulation of mitochondrial function; and c) the effects of energy depletion on AMPK activity. Binding of insulin to the transmembrane insulin receptors leads to phosphorylation and activation of downstream signaling resulting in translocation of glucose transporters GLUT1 and GLUT4 to the membrane [53,54]. Reactive oxygen species (ROS), FFA and pro-inflammatory cytokines (e.g. IL-6) have been shown to inhibit insulin signaling and contribute to insulin resistance [55–60]. Since the dynamics of insulin receptor phosphorylation and dephosphorylation are significantly faster than the dynamics of metabolic fluxes simulated in the model (such as glycogenesis, lipolysis, insulin secretion etc.), a quasi-steady state approximation was used to convert the differential equations representing the dynamics of insulin receptor activation into algebraic equations. Insulin sensitivity is expressed as the fraction of insulin receptors that are in the phosphorylated (active) state (S6 Fig and S7 Table).
In this study we have described, calibrated and validated a model of the physiological mechanisms of prediabetes development and diabetes onset in an individual. Our model simulates the dynamics of diabetes progression by mechanistically modeling the complex interactions between the following diabetes-related metabolic processes: the dynamics of pancreatic insulin production by beta-cells under normoglycemia and chronic hyperglycemia; pancreatic decompensation by cumulative damage to pancreatic beta cells due to glucotoxicity, and lipotoxicity; development of insulin resistance as a consequence of inhibition of insulin signaling by FFA and ROS; expression of insulin resistance through inhibition of translocation of GLUT1 and GLUT4 receptors in response to plasma insulin; energy production by oxidation of glucose, FFA, and amino acids; glycogenesis; glycogenolysis; lipolysis; de novo lipogenesis; proteolysis; protein synthesis; ketogenesis; gluconeogenesis; dynamics of ATP production and conversion to ADP as function of metabolism and energy expenditure; and hemoglobin glycosylation. Critically, the model is fit at the individual level and requires only common clinical measurements which allows for its translational use in clinical and public health applications.