A novel fast-slow model of diabetes progression: Insights into mechanisms of response to the interventions in the Diabetes Prevention Program

Summary

Several models for the long-term development of T2DM already exist, focusing on the dynamics of the interaction between glycemia, insulinemia and β-cell mass. Current models consider representative (fasting or daily average) glycemia and insulinemia as characterizing the compensation state of the subject at some instant in slow time. This implies that only these representative levels can be followed through time and that the role of fast glycemic oscillations is neglected. An improved model (DPM15) for the long-term progression of T2DM is proposed, introducing separate peripheral and hepatic (liver and kidney) insulin actions. The DPM15 model no longer uses near-equilibrium approximation to separate fast and slow time scales, but rather describes, at each step in slow time, a complete day in the life of the virtual subject in fast time. The model can thus represent both fasting and postprandial glycemic levels and describe the effect of interventions acting on insulin-enhanced tissue glucose disposal or on insulin-inhibited hepatic glucose output, as well as on insulin secretion and β-cell replicating ability. The model can simulate long-term variations of commonly used clinical indices (HOMA-B, HOMA-IR, insulinogenic index) as well as of Oral Glucose Tolerance or Euglycemic Hyperinsulinemic Clamp test results. The model has been calibrated against observational data from the Diabetes Prevention Program study: it shows good adaptation to observations as a function of very plausible values of the parameters describing the effect of such interventions as Placebo, Intensive LifeStyle and Metformin administration.

Keywords:

insulin – Glucose – Insulin secretion – toxicity – Glomerular filtration rate – glucagon – Nephrons

Introduction

Mathematical modelling is being increasingly used in diabetology, in order to help explain the mechanisms of normal and diseased control of the glucose-insulin system, both in short-term dynamical perturbation experiments and in the long-term development of the disease [1–6]. In particular, in order to understand quantitatively the interplay between insulin sensitivity, pancreatic β-cell responsiveness and β-cell population dynamics, several mathematical models of the long-term development of Type 2 Diabetes Mellitus (T2DM) have been formulated [7–14]. Besides the Topp [7] and deWinter [9] models, which were extensively commented upon in comparison with our previous model [10], in the past few years contributions on long-term diabetes progression modelling by Bagust [8], Ribbing [11], Boutayeb [12], Palmer [13] and Ha [14] have appeared in the literature.

Bagust et al. [8] developed a first-principles spreadsheet model linking several clinically observable physiologic variables, essentially centered on HOMA. No equations or computational formulas are reported in this publication, and the model cannot therefore be reproduced in any way nor can the assumptions be quantitatively tested. However, these Authors show plausible curves of evolution of several clinical indices over long times, specifically highlighting the eventual therapeutic failure of the sequence of progressively more intensive treatment regimens they simulate (sulfonilureas, metformin, insulin etc.).

Ribbing et al. [11] used the original Topp [7] model to represent the relationship between fasting glycemia, fasting insulinemia and β-cell mass, as these evolve over the years. In order to adapt model forecasts to patient subgroups from the GLAD and GALLANT clinical studies, these Authors introduce at some point in the life of each subject a discrete step-up (OFFSET) in the fasting glycemia “set-point” for β-cell mass dynamics. Following this regime shift the patients, over time, become diabetic. In fact, these Authors state that “…Reduced insulin sensitivity alone does not cause diabetes…”.

Boutayeb et al. [12] introduce yet again another modification to the Topp [7] model by adding an ϵ factor (taking values between 0 and 1) multiplying the glucose toxicity term in Topp’s model, arguing therefore that with ϵ = 1 their β-cell dynamics equation reduces to that of Topp, while with ϵ = 0 no genetic predisposition to diabetes exists, no glucose toxicity occurs and no diabetes develops. They also replace the mass balance equations for glucose and insulin with more general equations involving Michaelis-Menten terms, without however clarifying the physiological basis of the new formulation. They finally conduct a local stability study of the qualitative behavior of the solutions of the model thus defined.

Palmer et al. [13] used our previously published model [10] as a basis for focusing on the effect of IL-1β (Inter-Leukin -1 β) inhibitors such as anakinra [15, 16]. They derived their parameter values from available literature sources and came to the conclusion that most of the effect of IL-1β blockers is likely due to improvement of insulin secretion by existing β-cells, whereas appreciable changes in β-cell mass could take several years to occur. This is consistent with the interpretation of the findings by Sloan-Lancaster et al. [17].

Ha et al. [14] again started from Topp’s model[7] adding to it a dose-response shift in glucose-stimulated insulin secretion (governed by a dynamically varying coefficient γ), and an increase in maximal insulin secretory capacity under persistent hyperglycemia (governed by a dynamically varying coefficient σ). In this way they introduced an intermediate time-scale between fast glucose-insulin equilibration and slow diabetes evolution. They calibrated parameter values for their model in order to represent the time-course of the disease in experimental ZDF (Zucker Diabetic Fatty) rats. These Authors reached the conclusion that an emergent threshold in glycemia separates two stability basins, one (at lower glycemia) leading to compensation, the other (at higher glycemia) leading to manifest diabetes. There are clear similarities between this work and our previous model [10], in that both predict the eventual fast acceleration of the development of diabetes once insulin hypersecretion slows down, and both account for the bistability of achievable compensation (in case of normoglycemia and/or maintained pancreatic reserve) or development of disease (in case of prolonged hyperglycemia with reserve exhaustion due to glucotoxicity).

We previously introduced a mathematical model of T2DM [10], which was able to replicate acceptably well the observed time courses of fasting glycemia and diabetes incidence [18], when compared quantitatively with the results of the non-intervention group in a large study of individuals at-risk for T2DM [the Diabetes Prevention Program (DPP) study [19–21]]. In addition to exploring physiological processes and the natural history of T2DM, such models can be utilized to predict the long-term effects of pharmacological or non-pharmacological interventions. To this end, our model was also able to replicate the effects of metformin, troglitazone, and intensive LifeStyle modification arms in the DPP study [18].

One limiting feature of our previous model, as well as other models which have appeared in the literature [7–14], was that the considered model structure did not allow for an independent assessment of post-prandial or post-OGTT (Oral Glucose Tolerance Test) glycemia, and lacked the ability to represent separately hepatic and peripheral insulin sensitivity. It was felt that this limitation was of practical importance in the use of the model for physiologically-based clinical trial simulation in diabetes and the exploration of potential effects of diverse pharmacologic mechanisms. We have therefore developed a new version of the Diabetes Progression Model (DPM). In this context, the mathematical approach used to reconcile daily glucose homeostasis (fast) with compensatory evolution of β-cell population and pancreatic reserve (slow) in the previous version of DPM (i.e. the consideration that in slow time fast variables are essentially at equilibrium, while in fast time slow variables are essentially constant) proved to be insufficient in capturing the possible effects of variable daily glycemia on long-term compensation mechanisms. Also, in the course of the study leading to the final form chosen for the new version of the model (DPM15), collateral issues arose, such as the representation of renal glucose reabsorption, or the possible modeling of Type 1 Diabetes Mellitus (T1DM) and its therapy. All of the above considerations determined a substantial restructuring of the whole modeling approach, which had the additional benefit of making the model versatile, allowing us to simulate with it other experimental conditions and perturbations (e.g. OGTT and Euglycemic Clamp).

The goal of the present work is to detail the assumptions underlying the functional form of the new model, named DPM15; to justify the numerical values assigned to its parameters; and to show model forecasts corresponding to all of the endpoints that were recorded in the DPP study (fasting glycemia and insulinemia, 30-min. glycemia and insulinemia during OGTT, 2-hr glycemia during OGTT).

Materials and methods

The model to be described is named DPM15 since it is the 15th version of Diabetes Progression Model we built in the ongoing effort to capture the relevant features of the development of this disease. In the following, the model equations defining the model variables will be introduced and discussed. The model variables and parameters are summarized in Tables 1 and 2 respectively.

In the Parameters table (Table 2) the Value reported refers to the calibrated “baseline” value, the hypothetical value in the untreated DPP cohort. All DPP treatment arms shared these same general parameter values except for those specific parameters (Table 3) embodying the differences between groups induced by the different treatments.

Aspects of the model that are retained from the previous model are discussed here only briefly. Further details can be found in DeGaetano et al. [10]. This holds in particular for the choice of parameter values as distilled from the literature or from the adaptation of specific sub-models to available observations.

In the present work, a rather mechanistic, physiologic approach has been followed, explicitly computing the daily time-course of glycemia, insulinemia and related variables, at each step of the numerical integration of the slow system (say, every month), using a model of fast glucose homeostasis. From this daily portrait, the desired target descriptors of glucose/insulin control at that moment in slow time may be computed (e.g. average daily insulinemia or glycemic Area Under the Curve AUC after meals), and their value used to affect the further evolution of the slow pancreatic compensation system. In this way, a kind of alternating-step solution is obtained: the overall compensation status as reflected by the current values of the slow model variables (e.g. β-cell mass) is used as framework for the reconstruction of the daily profiles of fast model variables. In turn, fast variables can impact slow variables (e.g. potential effects of daily glycemia on β-cell replication) and these effects are allowed to integrate over another slow time interval. The model structure described above portrays explicitly the interplay of a slow model for the evolution of pancreatic compensation and changes in insulin sensitivity, with a fast model for immediate glucose/insulin homeostasis. The possibility of availing oneself of a (simplified but) complete description of glucose homeostasis at any time in the life of the subject allows the investigator to study new issues, which a solely slow model could not address: such are, for instance, the implications of assuming β-cell glucose toxicity to be due to elevated peak rather than average glycemias. Another new issue that can be addressed by this change of strategy is that of explicitly modeling, at any given time in the patients life, some clinical indices (e.g. HOMA-IR and HOMA-B) and the expected response of the subject to some commonly used experimental perturbation procedures, such as the Intra-Venous Glucose Tolerance Test (IVGTT), the OGTT or the Euglycemic Hyperinsulinemic Clamp (EHC). In this way the investigator has the possibility of correlating directly the hypothesized lifetime evolution with observed experimental measures and commonly employed clinical indices.

Given the slow-fast structure of the model, in order to avoid confusion when referring to “time”, the letter t (months) has been reserved for slow time, indicating the evolution of overall compensation on a scale of months to years, while the letter τ (minutes) has been used to indicate the evolution of the glucose homeostasis mechanism after acute perturbations such as meals.

Slow model

β-cell mass (B)

Eq 1 defines the variation of the β-cell population B as depending on a (variable) net replication coefficient k_BB and on a possible additional coefficient of β-cell mortality μ, through which it is possible, for instance, to represent the early auto-immune development of T1DM or the postulated cytotoxic effects of cytokines or lipid species in T2DM. Besides the introduction of μ, the form of the equation differs from the previous model in that a limiting ‘carrying capacity’ for B has been introduced, transforming the previous exponential model [10] into the present logistic model. The value of the carrying capacity B_max has been set at 4 billion cells, i.e. four times the normal value of approximately 1 billion β-cells previously estimated on the basis of several literature sources [22–27] and assumptions about cell size. Pregnancy and obesity may be associated with a doubling of β-cell mass [28, 29]. A maximal four-fold increase does not seem unreasonable and may allow ample space for normal variation.

β-cell net replication rate (k_BB)

A fundamental assumption of the present model is the dual effect of glycemia on β-cell replication: β-cell population dynamics is affected by hyperglycemia through both a direct short-term stimulation [30–32] and an indirect longer-term inhibition of net replication, possibly due to glucose toxicity “exhausting” β-cell replication reserve [22, 33–36]. The same assumption underlies our previous model [10] and is similar to what was postulated by Topp et al. [7].

Eq 2 relates the replication coefficient k_BB (β-cell net replication rate) with (short-term) glycemic stimulation:

The net replication rate is increased by hyperglycemia (with a nonlinear, saturating mechanism depending on pancreatic replication reserve η) above a minimum rate k B B m i n, taken to be negative. In this way, allowance is made for both positive and negative oscillations of β-cell net replication rate, translating into increments and decrements of β-cell mass. Notice that in this formulation G_B is some function assumed to best describe the aggregated effect of the daily glycemic variations in stimulating β-cell replication: for the current implementation of the model, G_B has been taken simply as average daily glycemia.

β-cell replication reserve (η)

The possible excursion of the β-cell net replication rate k_BB has been termed η = k B B m a x - k B B m i n. It is non-negative and is governed by Eq 3:

This excursion or range is a measure of the maximal replication rate of the pancreatic β-cells as depending on the current state of pancreatic “health”: in other words, it represents pancreatic β-cell replication reserve. It has some starting value η₀ and then increases to a maximum η_max or decreases towards zero depending on glycemia levels. Hyperglycemia is supposed to be toxic to β-cell replication [33, 34], hence sustained hyperglycemia will lead to a decrease of η. Notice that setting k_η at zero, we would assert that pancreatic reserve necessarily decreases with age.

The function G_η is computed as the integrated mean over 24 hours (computed from the fast daily model) of the glucose toxicity produced by the varying glucose concentrations throughout the day. Glucose toxicity is monotonically increasing with glycemia, has been calibrated on TUNEL (Terminal deoxynucleotidyl transferase dUTP nick end labeling) data [37, 38] and indicates η suppression as a fraction of normal (normal = 1 = 100% at fixed 5.5 mM glycemia). In other words, at each instant in fast time the current glycemia determines the current glucose toxicity following Maedler et al. [37, 38] (according to an increasing, saturating Hill function with 0 toxicity at 0 mM glycemia, toxicity 1 at 5.5 mM, toxicity 3.5 at 30 mM and asymptotical toxicity 4 at infinite glycemia); toxicities throughout the day are integrated and then divided by the day’s duration in order to obtain the integrated average toxicity G_η for that day.

Pancreatic glucose toxicity

In Eq 3, k_XηG is the coefficient expressing the intensity of pancreatic glucose toxicity. It may be allowed to vary over time, starting at some value k_XηG0 and being possibly modified by therapy.

Glycemia-independent β-cell mortality rate (μ)

Independent factors, such as inflammation or auto-immune processes, may contribute to β-cell mortality independently of glycemic levels. Eq 4 allows the expression of this excess mortality rate μ over time, starting from a level μ₀ (assumed to be zero in health) and progressing sigmoidally towards a maximum additional mortality rate μ_max:

This equation is used, for instance, when using the model to represent the time course of Type 1 Diabetes Mellitus (T1DM), when additional β-cell mortality leads to rapid disappearance of most of the β-cell population over a relatively short time period.

Spontaneous pancreatic recovery rate

In Eq 5, k_η indicates the current ability of the pancreas to increase (recover) its β-cell proliferation rate. This ability is assumed to vary linearly between t₀ and t_end (from young age to end of life) so as to allow the possibility of representing non-constant spontaneous pancreatic recovery rate throughout the subject’s lifetime. This is of interest when considering that natural aging may reduce, over time, β-cell proliferative capacity [39, 40].

Fasting plasma glucose concentration

Given the current slow state, a Daily run of the fast model, inclusive of meals etc., when starting at glycemia G_f calculates glycemia G_f24 exactly 24 hours later. In the long run therefore, G_f will tend to the value G_f24 over slow time, and this is assumed in the model to occur at a rate k_GG. In general, we do not expect G_f to be equilibrated within a day, in particular when glycemia is in a rising phase (e.g. in the pre-diabetic state). Writing Eq 6 we recognize that slow G_f(t) may differ from the corresponding slow G_f24(t) and assume that G_f tends to G_f24 (in the pre-diabetic state example, we assume it increases towards G_f24) with rate k_GG. Equilibrium may in fact never be attained, because as G_f converges to G_f24, G_f24 itself may shift due to the concurrent change in slow model variables. It is to be noticed here how the concept of the convergence of G_f to G_f24 provides the link between fast time and slow time glycemia, so that changes in G_f24 due to fast-time dynamics in fact drive the evolution of the whole system of glucose homeostasis in slow time.

Fasting serum insulin concentration

The fasting serum insulin concentration is computed as the fast-equilibrium value at current β-cell mass and fasting glycemia values. The function ψ I g l u c ( B , G ) indicates glucose-driven pancreatic insulin secretion at given β-cell mass B and driving glycemia G:

In Eq 8 insulin secretion depends therefore on current β-cell mass B, on k I B m a x, the glucose sensitivity of the existing β-cells, and on a saturating stimulus provided by increasing glucose concentrations.

Difference of exponentials

In the following, we use a standard difference-of-exponentials functional form to express the time-course of the action of a given therapy on some control variable, i.e. on some key state variable determining the evolution of the whole system:

The difference of exponentials defines a curve starting from 0 at time 0, rising to a peak amplitude Ampl and redescending to zero (α > β) as t → ∞, with respective rates of increase and decrease defined by the relative values of α and β, the larger the α the faster the rise, the larger the β the faster the fall.

Using the difference of exponential formula, it is possible to specify a progressive rise of the effect of therapy on some control variable, starting at some therapy initiation time t_thx.

Notice how, by changing the values of α and β appropriately, it is possible to specify very different time courses of the effect: for instance, with α very large and β very small, the difference of exponentials can approximate a step function, with immediate increase and essentially permanent effect.

Insulin secretory function

We indicate with k ˜ I B m a x the “natural” maximal insulin production rate per million β-cells, which can be made to vary (decrease) over time in order to express a possible decay of functional secreting ability:

Insulin elimination rate

Over time, the apparent elimination rate of insulin from serum may be supposed to vary. In fact, it is presumed to be slowly decreasing over lifetime [41] (i.e. k_XIend < k_XIstart):

Peripheral insulin sensitivity

In order to represent variable maximal insulin-dependent tissue glucose uptake (Peripheral Glucose Disposition, PGD) rate, we may take the “natural” value of k X G I m a x to be represented by:

Expressing the possible decay of insulin sensitivity over time in the same form as that used above for the time-course of insulin secretory function offers a large flexibility in the kind of behavior that can be represented (early or late start, fast or slow decrement, large or small eventually preserved function etc.). The “natural” value of peripheral insulin sensitivity k ˜ X G I m a x, determined by a combination of genetic factors, lifestyle etc., is then possibly modified by therapy, with effect expressed by f_exp(LyKxgiCurr, α_LyKxgi, β_LyKxgi, t − t_thx) at any given time after the initiation of therapy, so as to yield the actual current value k X G I m a x.

We note that hepatic glucose uptake is here comprised in the total Peripheral Glucose Disposition, and is considered to be insulin sensitive in the normal individual. There has been debate about this point. It is true that the hepatic glucose transporter (GLUT2) is not responsive to insulin [42]. It is also true that some studies have not shown stimulation of hepatic glucose uptake by hyperinsulinemia. Nevertheless, compelling evidence exists for insulin stimulation of hepatic glucose uptake [43–45]. While the assumption, that the effect of insulin on hepatic (or splanchnic) glucose uptake is similar to insulin-dependent glucose disposition in other tissues, is clearly an oversimplification of the differential mechanisms of insulin stimulated glucose uptake in liver and peripheral tissues, the above literature supports similar effects and overlapping concentration-responses in liver and other insulin-responsive tissues. The assumption, that we may model with a single overall effect the composition of peripheral and hepatic glucose uptake, seems not unreasonable in the light of the above considerations.

Hepatic insulin sensitivity

In the present version of the model, “Hepatic Insulin Sensitivity” refers to the action of insulin to decrease hepatic glucose output (HGO) and also kidney glucose output to the extent that it is significant. More specifically, hepatic insulin sensitivity (λ_GI) is formalized as the rate of the exponential decay of hepatic glucose output with increasing glycemia and insulinemia (see Eq 20). This relationship is supported by experimental observations [46, 47]. Note that this expression only reflects the effect of insulin on hepatic glucose production and not on hepatic glucose uptake (so that HGO is always non-negative in our model). Instead, the effect of insulin on hepatic (or splanchnic) glucose uptake is assumed to be similar to insulin-dependent glucose disposition in other tissues and is therefore incorporated in the k_XGI (peripheral insulin sensitivity) term. The study by Basu and colleagues [45] supports this assumption.

In order to represent variable hepatic insulin sensitivity (decreasing over slow time), we may take its “natural” value λ ˜ G I to be represented by:

Stomach emptying rate

k_JS is the rate expressing apparent first-order stomach emptying. In order to allow for the possibility that its value changes over time (e.g. it may be affected by drugs slowing gastric emptying), it is represented as a variable rather than a parameter. Its healthy constant average value is indicated by

We therefore endowed the model with the ability to represent in general, for each treatment and for each affected control variable, arbitrary rates of onset of effect and arbitrary rates of loss of the same. In the current version of the model, the five control variables for which this scheme has been implemented are η (β-cell replication reserve restoration), k_XGI (peripheral insulin sensitivity), λ_GI (hepatic insulin sensitivity), k_IB (insulin secretory function)and k_JS (glucose absorption rate from the GI tract).

Daily model

In order to describe the fast variations of glycemia over minutes, across a time span typically not exceeding the length of one day, a fast, “Daily” model is utilized. The present version of the fast model borrows heavily from previous rapid glucose-insulin control models, which proved to be effective and parsimonious representations of fast dynamics during acute perturbation experiments (both Intra-Venous and Oral Glucose Tolerance Tests, IVGTT and OGTT) [48–51]. As mentioned above, we denote “fast” time with the greek letter τ in order to underscore the difference between “fast” variables and slowly varying phenomena, which change over the months of “slow” time t, as in the previous equations.

Gastrointestinal glucose transit and absorption

The variation of S(τ), glucose content in the stomach, is described as a simple first-order linear elimination after impulsive loading corresponding to the meals:

Similarly, the rate of change of glucose content in the intestine (indicated as “Jejunum”, J), based on delivery from the stomach and exit into the circulation, is expressed as a simple linear ODE:

The rate of glucose appearance in the systemic circulation, following intestinal absorption, is a fraction of the disappearance rate of glucose from the “jejunum” (since some of the absorbed glucose is either utilized or stored by the gut and/or liver):

Fast plasma glucose concentration

It is first of all to be clarified that “Fast” is here meant as the opposite of “Slow” and not as characterizing an early morning or otherwise “fasted” state. The point here is that while some characteristic daily glucose concentrations, such as for instance early-morning fasting glycemia G_f(t), vary over the months of a lifetime, glucose concentration G(τ) can also be thought of as varying, minute-by-minute, over the period of an experiment or over a given day in the life of the individual. The mass balance considerations determining the evolution of plasma glucose concentration over a given day are described by the following Equation, contemplating insulin-dependent glucose elimination (saturating), possible insulin-independent linear glucose elimination, and contributions to the variations of plasma glycemia deriving from liver gluconeogenesis and glycogen lysis, from the action of glucagon (L), from brain glucose oxidation, from the appearance of foodstuff-derived glucose, and from the loss of glucose in the urine:

At any slow time t, we assume we start our daily evolution from a fasting glycemia value G_f(t). Peripheral insulin sensitivity depends on the prevailing levels of insulin (in the sense that the marginal utility of an increase in insulin decreases as levels of the hormone increase). When considering tissue insulin sensitivity to be nonlinearly increasing with increasing insulin concentration, we wish in any case to retain the common notion of a (linear) insulin sensitivity index at zero insulinemia, k X G I 0, with the same measurement units and with comparable numerical values with respect to previously published insulin sensitivity indices (see Panunzi et al. [50] for a comparative discussion of two such indices). In the present formulation, k X G I 0 expresses the slope of the non-linearly increasing peripheral tissue insulin sensitivity at I = 0. Its value should be around 1.e − 4[/min/pM] for normal individuals [48, 50] and somewhat higher for fit, athletic individuals. At a value k X G I 0 = 1 . e - 4 / m i n / p M and with a reasonable insulinemia at half-effect I K X G I 50 = k X G I 0 m a x k X G I 0 = 500 p M [52] we find a current value k X G I 0 m a x = I K X G I 50 k X G I 0 = 500 · 1 . e - 4 = 0 . 05 / m i n, which seems very reasonable (i.e. maximum tissue uptake of glucose at hyperinsulinization, or in other words maximal clearance of glucose from plasma at extremely high levels of insulin, is approximately 5%/min).

We could directly extend the model by considering different meal compositions at breakfast, lunch and dinner, with ensuing differences in gastric emptying rate and glucose absorption rate. This could be attained by considering N_foods types of foodstuff indexed by j at each meal m: different foodstuffs would load separate S_j and J_j compartments, with transit governed by M j g l u c, f_GJj, k_GJj, τ_j, and k_JSj parameters. For simplicity, the uniform meal composition formulation is retained. The ur renal glucose elimination rate is provided by the concurrently running nephron model (see below).

Fast serum glucagon concentration

The fast model contemplates glucagon as an index of the overall counterregulatory response.

Glucagon plasma concentration (L) is supposed to undergo first-order linear elimination [53, 54]. This has in fact been demonstrated to be true in man [53] as well as in dogs, with a suggestion that the process may be saturable at pharmacologic concentrations [55]. Glucagon plasma concentration has also been shown to increase above some minimum (determined by some measure of continuous production), in an exponentially increasing fashion as glycemia decreases: a large set of experimental data were fitted to show this relationship [56]. We note that there is literature support for the idea that insulin cannot override hypoglycemia in suppressing glucagon secretion [57] (Cavallo-Perin et al. [58] showed only 20 − 40% suppression of fasting glucagon levels with euglycemic clamps when plasma insulin was increased to 350 pM, and similarly Elahi et al. when using 700 pM [59]). A determinant role for insulin would seem to be indicated by the observation of hyperinsulinemic coma with low glucose and glucagon concentrations: the present form of the model has however been considered synthetic and sufficient to describe relative insulin deficiency situations, such as occur with the development of T2DM. Clearly, insulinomas and similar conditions would be poorly represented by the current form of Eq 21. Notice also that we assume the dynamics of glucagon to remain unchanged over slow time. Even so, there are situations [60] where glucagon release is substantially reduced in T1DM, though the mechanisms for this are not entirely clear [61].

Fast serum insulin concentration

Insulin kinetics in the short period may be approximated, following mass-balance considerations as

Insulin is assumed to be eliminated from plasma in a first-order, linear fashion. The entry of insulin into plasma derives from glucose-driven, saturable pancreatic secretion. For the purpose of the present model we do not distinguish peripheral from portal insulin concentration, therefore hepatic insulin action will be made to depend on peripheral insulin concentration itself (whose effect will be apparently larger with larger hepatic insulin extraction). The equations above express insulin secretion as depending linearly on the currently available β-cell mass and nonlinearly, in a saturable way, on glycemia. Notice that, in order to represent incretin contribution to the stimulation of insulin secretion, the effectiveness of glycemia is supplemented (saturably) by current intestinal glucose contents, via a proportionality constant f_IJ.

Nephron model

The concept of Renal Glucose Threshold is well rooted in common medical and diabetological practice. According to this concept, renal elimination of glucose occurs when, glycemia having exceeded some threshold G_thresh, glucose delivery to the nephrons exceeds their reabsorptive ability. While physiologically very plausible, the concept has unfortunately been translated, typically, into a mathematical formulation stating that if at some time τ glycemia G(τ)>G_thresh, then at that time there will be glucose loss in the urine at a rate proportional to the difference G(τ) − G_thresh. The problems arising from this simplistic interpretation of the glucose threshold principle have been treated in detail elsewhere [62]: in the same work, an alternative, partial differential formulation of the principle, has been proposed, which does not suffer from the problems described. In the current context we therefore use the same partial differential equations approach used in that publication in order to build a simple nephron sub-model, able to realistically represent glucose elimination in the urine produced by glycemia oscillations over time. For the reasons explained in detail in the work referred to above [62], we believe that the simpler renal-threshold sub-model would predict glucose renal elimination inappropriately over the very glycemic range around the supposed threshold, range that is repeatedly entered by glucose concentrations in the course of the day, particularly in insulin-resistant subjects. Having a more advanced nephron sub-model available, we took advantage of it, incorporating it in the overall disease progression model. It is to be noticed that in the nephron model the considered time τ in minutes is the same “fast” time as in the Daily model above: in fact, the Nephron and Daily sub-models progress in parallel through a typical short time period (a few hours), with the Daily model determining at each discretization step the current glycemia, used by the Nephron model to determine at each step the corresponding urinary elimination, used then again by the Daily sub-model to determine glycemia variations and the resulting glycemia, and so on.

Density q of glucose at time τ at nephron level z

The nephron model considers amounts of glucose in small (infinitesimal) segments of an idealized single tubule, representing the collection of glucose-reabsorbing proximal and distal nephron tubular sections. Henceforth, all quantities are to be understood as referred to Kg body weight. The density q of glucose quantity Q with respect to normalized tubule length z (z ∈ [0, 1]) varies over time and along the tubule, q = q(τ, z):

The variation of glucose density over time at some point z in the tubule is given in general by a transport equation with diffusion D_U along the tubule, with advection driven by the flow-rate ϕ_U (variable along the tubule), and by saturable extraction from the tubule lumen operated by lining cells. All three effects are expressed in terms of concentration C = Q/V = q/v. The notation k G U m a x 1 indicates that the maximum total transport by the tubule is to be divided by the total z-length of the tubule itself, which in the present normalized case equals exactly 1. We further assume that glucose reabsorptive capacity is essentially zero at the end of the tubule, so that e−λQZ⋅1≈0 and ∫ 0 1 λ Q Z e - λ Q Z z d z ≈ 1, so that k G U m a x does indeed represent the maximum glucose absorptive capacity of the entire tubule system.

In this formalization, the volume density v (akin to the total tubular cross-sectional area at some level in the tubular system) is defined as the density of water amount with respect to normalized tubule length, reflecting progressive water reabsorption along the nephron, and is proportional to 180 L/day ultrafiltration at z = 0 and to 2 L/day urine output at z = 1 (in a 70-kg person). An approximation to the volume density profile v(z) at given depths down the renal tubules is obtained by hypothesizing an exponential volume decay over the length of the tubule, with entry proportional to the flow of ultrafiltrate and exit proportional to the flow of urine:

The initialization profile q(0, z) is only a first rough (exponentially decaying) approximation to the quantity of glucose in tubular water: before starting with the Daily / Nephron numerical integration, the Nephron model is run at G = G₀ for as long as necessary to reach convergence in the q profile. The boundary condition q(τ, 0), on the other hand, is given by the glucose concentration in plasma at time τ times the volume density at z = 0.

Given quantity and volume, the concentration C of glucose in the pre-urine is algebraically determined:

Finally, the rate of glucose elimination in the individual at time t is given by the most distal/caudal (pre-)urine glucose concentration (i.e. glucose concentration at z = 1) times the total urinary flow rate:

Index variables

One fundamental assumption of the present model is that the variation over months and years of the “slow” homeostasis variables (such as β-cell population size or pancreatic β-cell replication reserve) may depend not only on monthly averages of fast variables (such as glycemia or insulinemia), but also on other slow functions of these same fast variables, possibly derived from the explicit computation of the time course of the fast variables over a representative day. Examples of these slow indices derived from fast dynamics and possibly affecting disease progression are average glycemia, glycemic variability, or post-meal glycemic peaks. Commonly used clinical indices such as the HOMA-IR or the EHC M-values are also computed as slow index variables. The weighted daily glucose toxicity determining η suppression, as fraction of normal (normal = 100% at fixed 5.5 mM glycemia) derived from the TUNEL [37] study, also belongs to this class. The definition of some such index variables is straightforward, assuming the daily time course of the necessary fast variables to be available: mean daily glycemia and insulinemia; their standard deviations, minima and maxima; their value at relevant times (e.g. baseline, 30 min, 60 min, 2 hours after administration of the glucose bolus during an Oral Glucose Tolerance Test). Given the availability of the necessary ingredients, the computation of some clinically relevant slow indices (HOMA-IR, HOMA-B, insulinogenic Index) is also straightforward. Finally, a two-step Euglycemic Hyperinsulinemic Clamp experiment (120 min at 100 pM insulinemia, followed by 120 min at 420 pM insulinemia, replicating conditions that have been used previously to assess respectively hepatic and peripheral insulin sensitivity) can be conducted on the virtual subject at any time, and the corresponding low- and high-insulinization Glucose Infusion Rate (GIR) values, expressed in mg/kgBW/min are also slow indices of interest and are indicated as ClampM1 and ClampM2.

Parameter assessment

The model proposed here combines disease progression over slow time with daily absorption, metabolism and renal elimination in fast time. Parameter values for the slow time component have been derived from the in-depth assessment undertaken with the publication of the previous model [10] wherever applicable. This includes fasting glycemias and insulinemias; deduced characteristics of β-cell replication (but see also [18] and above comments on the switch from linear to saturable model of β-cell population dynamics); pancreatic β-cell replication reserve capacity; glucose toxicity; insulin secretion (per millon β-cells) and elimination rates; glucose effectiveness; production and decay rates of glycosylated haemoglobin. However, the model introduces separate insulin sensitivity components for hepatic (and possibly renal) and peripheral (mainly referred to muscle and adipose tissue) insulin sensitivity, and assumes that both components may vary over slow time (years or decades) depending on genetics, lifestyle, intervening diseases etc.

Furthermore, insulin independent glucose clearance (mostly due to brain glucose consumption) has been explicitly introduced in the fast glycemia equation.

Parameter value assessment, unless specifically discussed below, follows the same approach as delineated for the previous model [10]. The majority of parameters were obtained from in vivo studies on adult subjects, pediatric or developmental data were not considered. Typical values and ranges at time t₀ were generally taken from data obtained in young, healthy adults (age 18-30 years): while t₀ is defined as the time of birth for the present work, we make no attempt to actually replicate the physiological variations of the glucose-insulin system along infancy and childhood. Rather, we let simulations through the first 18 years of age run on early adulthood parameters, thereby reaching a normal or healthy steady state at age 18, and start perturbing the model (e.g by imposing a progressive decrement in insulin sensitivity) after age 18 years.

Parameters describing the average evolution of insulin resistance in the whole patient population were calibrated in order to achieve starting conditions at age 50 consistent with the overall average starting conditions in the DPP study, where average age at enrollment was approximately 50 years. More specifically, we simulated an average subject presenting as pre-diabetic at age 50 by introducing a progressive decline (after age 18) in peripheral and hepatic insulin sensitivity, as well as a modest, progressive decline in insulin secretory function (see the functional description of k_XGI, λ_GI and k_IB in Eqs 15, 14 and 12, Table 2).

The next step was to calibrate those parameters describing the rates of onset and decay and the size of the effects produced by the three experimental maneuvers considered (Placebo, Metformin, LifeStyle), based on our understanding of the likely physiological effect of the three treatments, in order to reproduce observed time courses for fasting and post-prandial glycemias (Table 3). We assumed, in particular, that the DPP treatments would likely impact hepatic and/or peripheral insulin sensitivity and would not have a direct effect on insulin secretion; moving from this assumption, by changing the rapidity of onset, rapidity of decay and size of the hepatic and peripheral insulin sensitivity improvements for each treatment, we explored the ability of the model to reproduce simultaneously the size and shape of the observed time courses of fasting and post-prandial glycemias and insulinemias.

The selection of parameter values was subjective (visual assessment). While a large number of simulations were performed, no systematic exploration of the whole (high-dimensional) parameter space was conducted, but rather the parameter value combinations reflecting our understanding of the likely effects of the various interventions were marginally adjusted to better approximate the data points. For example, when considering the adaptation of the model predictions to the observed DPP averaged data points, we had to improve peripheral insulin sensitivity, over predicted no-intervention levels, by a maximum of about 10% for Placebo and Metformin and 22% for the Intensive LifeStyle group; further, while rates of decay of effect were apparently similar for the three groups, onset of effect was much faster for Intensive LifeStyle than for Placebo or Metformin.

In the present work the emphasis was on showing the reasonableness and robustness of the model rather than on estimating effect size. Since the model is relatively large with respect to the independent sources of information from the DPP study (which, together with likely correlations among parameters, would have made the model a-posteriori unidentifiable), we eschewed the use of formal optimization of some classical loss function for statistical parameter estimation. While the model structure is apparently correct in that it can replicate observations, we cannot therefore assess the variability of the parameter estimates, cannot construct confidence intervals around them and cannot offer measures of Goodness-of-Fit.

DPP dataset

The objectives and results of the Diabetes Prevention Program (DPP) study have been extensively described elsewhere [20, 21] Briefly, in this study the primary aim was that of evaluating the incidence of T2DM in an at-risk population randomized to placebo (n = 1082), intensive lifestyle modification (n = 1079), metformin (n = 1073), or troglitazone (n = 585). From 1996 to 1999, the study enrolled adult subjects with elevated Fasting Plasma Glucose (FPG from 5.6 to 7.7 mmol/L before June 1997; FPG from 5.3 to 6.9 mmol/L after June 1997) and 2-hr-Glucose (G2h from 7.8 to 11.0 mmol/L) during a 75g Oral Glucose Tolerance Test (OGTT), as well as elevated Body Mass Index (BMI ≥ 22 kg/m in Asians, BMI ≥ 24 kg/m otherwise). We obtained the original dataset from the DPP study through application from NIH-NIDDK (February 2008 Full Scale data release; data request No.608, see Acknowledgements). The data were fully anonymized as supplied by NID-NIDDK, not only was identifying information eliminated from the data set, but also subjects’ ages were not reported as recorded, but only in five-year classes. We excluded from this dataset a small number (59) of patients who presented with FPG in the diabetic range at entry into the study (“Entry Diabetics”) and proceeded then to compute averages of FPG (G_f) at entry and at each six-monthly observation thereafter, as well as averages of 30-minutes and 2-hour glycemia (G_30m, G_2h) and baseline (FSI) and 30-minutes insulinemia (I_f, I_30m) after Oral Glucose Tolerance Test (OGTT), at entry into the study and at each yearly interval thereafter.

We did not receive any special access or privileges to the data: interested researchers will be able to access the data in the same manner as we did. Interested researchers can replicate our study findings exactly and in their entirety by implementing the equations constituting the model described in the Methods section, populating the implementation with the parameter values reported in the Tables; and finally plotting the resulting model predictions together with the averages of the DPP study data at each time point.

Implementation

The computational engine of the model has been implemented in C++ (Microsoft^® Visual Studio Community Edition 2017), with a MATLAB^® graphical front-end (MATLAB version R2009b, The MathWorks Inc.). The model engine is also accessible both for guest researchers use (through a browser HTML interface) and as a web-service for guest machine-to-machine use (via a WSDL) at the address biomatlab.iasi.cnr.it/models/login.php.

Results

Fig 1 shows the hypothesized disease progression in terms of the evolution of peripheral insulin sensitivity (Panel a, substantially decreasing with age), hepatic insulin sensitivity (Panel b, substantially decreasing around ages 40 to 60), and insulin secretory ability per unit β-cell mass (Panel c, mildly decreasing with age). This is a combination of original modifications (attributable to advancing age, lifestyle factors and dietary habits), which can explain, through the model, the observed average changes in all three treatment arms.

Figs 2 through 6 show the time course, over an interval of slow time spanning the study period, of the endpoints measured in the DPP study, together with the computed time-courses of the corresponding model variables.

Fig 7 shows daily and OGTT time courses of both glycemia and insulinemia as predicted by the model for a representative virtual DPP subject. From these, relevant clinical indices can be computed.

Figs 8 through 12 show the time course, over the interval of slow time spanning the study period, of some commonly employed clinical indices of diabetic (de-)compensation: HOMA-IR, HOMA-B, insulinogenic index and 1^st and 2^nd phase clamp M indices.

Fig 13 finally shows the model-predicted time course of the five DPP-observed endpoints as well as of glycated hemoglobin, over the whole adult life span of representative virtual subjects, without treatment and undergoing each of the three examined DPP treatment protocols.

We would like to underscore that, of the many variables whose time course is predicted by the model, only five (fasting glycemia, OGTT glycemia at 30 and 120 minutes, fasting insulinemia and OGTT insulinemia at 30 minutes) were actually observed within the DPP study. All of the other variables were not observed (including the complete daily variations of glycemia and insulinemia) and their time courses are inferred by the structure of the model and the fit of the model, with given parameter values, with the actually observed variables.

In the figures, numerical simulations obtained with the DPM15 model are shown together with corresponding sample averages obtained from the Placebo, LifeStyle and Metformin arms of the DPP study. In order to make predictions and average observations comparable, we simulated a representative (“average”) subject, reconstructing the subject’s life-trajectory from age 18 (assuming a completely normal insulin sensitivity and secretion profile at this age) until age 50 (approximately the average age of entry into the DPP study [20].

In the last column of Table 2, the values attributed to the parameters under the assumption of no intervention are reported. Table 3 reports those calibrated model parameter values which had to be set differently for the four treatment options considered (no intervention, Placebo, LifeStyle and Metformin). No change in insulin secretory capacity was hypothesized. Improvements in peripheral and hepatic insulin sensitivity are not apparent from the table, as they are determined by relative values of rate of onset, rate of decay and intensity of effect. The maximal improvements are zero for “No treatment”; approximately + 10% and <+ 1% for Placebo; approximately + 22% and + 7% for LifeStyle; approximately + 10% and + 7% for Metformin, all expressed as as percent of normal levels, respectively for peripheral and hepatic insulin sensitivity.

By comparing the values reported in Table 3 with the time-courses of the DPP endpoints in Figs 2 through 6 it can be appreciated how changes in peripheral insulin sensitivity translated mainly into changes in 30 min and 2 hr OGTT glycemia, while changes in hepatic insulin sensitivity translated mainly into changes in fasting glycemia, consistently with common physiological understanding.

The calibration exercise itself provided meaningful information in several ways.

First of all, it was supposed that patient’s glycemias would be progressively increasing, on the average, right before the time of entry into the study: were this not the case, it would be very difficult to imagine that patients could have arrived to pre-diabetic levels from normal levels during young adulthood and that they would then have developed increasing glycemias during the study.

Given the above assumption of increasing glycemias just before enrolment into the study, since even for the Placebo group glycemias were not observed to increase immediately after enrolment (indeed, for the placebo group the fasting glycemia at 6 months was on the average smaller, even if minimally so, than the average fasting glycemia at entry into the study), some clinically significant, even if small, Placebo effect on some of the control variables (peripheral or hepatic insulin sensitivity, gastric emptying rate) must have been present, so we had to conclude that a positive Placebo effect existed.

Secondly, it was not possible to reproduce the observed data assuming the same rapidity of onset and the same rapidity of progressive loss of efficacy of the treatments on each control variable, even assuming possibly different effect sizes. In other words, we could not fit the data assuming that say, metformin and placebo have the same rates of onset and loss of efficacy on peripheral insulin sensitivity, or that metformin and intensive lifestyle differ only in effect size with respect to their impact on hepatic insulin sensitivity. The converse was also true: it was not possible to explain the observations for some treatment, say metformin, by assuming same rapidity of onset and loss of efficacy on different control variables, e.g. on both hepatic and peripheral insulin sensitivity. We had therefore to hypothesize differences in rate of onset and rate of decay of effects for different treatments on the same control variable as well as for the same treatment on different control variables.

The best fits with the observed data for metformin treatment were obtained when we assumed that this intervention acts primarily on hepatic insulin sensitivity. This determined a glucose lowering effect largely limited to fasting glycemia. Interpretation and observation are consistent with studies showing that metformin acts to reduce fasting glucose mainly by decreasing hepatic glucose output and has little effect on peripheral insulin sensitivity [63–65]. Metformin’s specific effect on hepatic insulin sensitivity appeared to be transient; its more modest, but somewhat longer-lasting effects appeared indistinguishable from placebo. The transient nature of metformin’s effects is consistent with previous reports [66–68].

The observed time-course of the recorded variables for the ILS arm was best replicated by the model when assuming an effect on hepatic insulin sensitivity of about the same order of magnitude, and with the same rates of onset and decay, as that hypothesized for Metformin. In addition, however, the observations on the ILS arm were consistent with a very substantial effect on peripheral insulin sensitivity, increasing progressively from the start of treatment, peaking in about two years on average, and, while decreasing thereafter, being substantially maintained over the course of a few more years.

Over the 5-year period of the study the predicted change in β-cell mass in a typical subject was greater than 10%, the model predicting in fact the fastest decrease in β-cell mass right at this period in the lifetime of the kind of pre-diabetic subjects enrolled in the DPP study. The treatments studied did not appreciably influence, over this time bracket, this apparent rate of loss: model simulations suggest that earlier, more aggressive treatment aimed at drastically reducing insulin resistance and possible additional causes of β-cell mass decrement, such as systemic inflammation, would be beneficial.

Finally, while it may be imagined that model parameters could easily be varied so as to arbitrarily change model output, the present model’s structure did not actually allow a wide variety of possible forecasts. Model predictions for the five DPP endpoints are tightly connected: the physiologic assumptions incorporated in the model equations and the relatively small number of free parameters (when compared with the large number of different observations) coerce the predicted curves into rather rigid patterns. Still, parameter values had to be found such that all five DPP endpoints were simultaneously matched when attributing to the different treatments plausible effect sizes and temporal evolutions. This turned out to be in fact possible. In other words, the model incorporates current knowledge in a coherent mathematical structure that cannot be bent to reproduce whatever arbitrarily specified behavior. Such mathematical structure, however, does translate physiologically plausible parameter value changes into realistic predicted time-courses, matching the observed averages in the DPP patient sample.

Discussion

Reasons for modelling diabetes progression

It is problematic to test alternative T2DM therapies over the long time-scale of the disease. What data sets exist cover at most a few years for each subject. Mathematical models of disease progression allow meaningful quantitative extrapolation beyond commonly observed time intervals: such models can be used to predict outcomes for novel therapies with anticipated disease-modifying properties and to help guide the design of clinical trials. Simulations can be used to guide sample size and treatment duration based on hypotheses about mechanisms of action (e.g. on insulin resistance or β-cell replication), about the magnitude of the effects, about the variability in both disease severity and treatment effects.

The validity of the predictions depends crucially on the robustness of model assumptions and on the mechanistic structure that results from them. The forecasts obtained are the direct expression of what diabetological knowledge is incorporated in the equations. We describe here a model, which we believe is robust and which has enough complexity to provide multiple different outputs allowing the study of multiple, diverse mechanisms of intervention.

Changes with respect to the previous model

The present work describes a radically improved version of our previous model [10] for the long-term development of Type 2 Diabetes Mellitus and its validation against observations on three patient groups from the Diabetes Prevention Program study [20, 21].

The consideration of a new model was prompted by the need to forecast multiple realistic, testable endpoints (fasting and post-prandial glycemia and insulinemia) rather than restraining oneself to discussing fasting,“representative” or “prevailing” glycemia and insulinemia. This need was met by abandoning the quasi-steady state assumption and by introducing hepatic as well as peripheral insulin sensitivity.

In the previous version of the DPM [10] a classical mathematical approach (singularly perturbed ordinary differential equations by means of small ϵ-parameter [69–71]) was employed to harmonize the slow progression of the disease (portrayed by β-cell mass or insulin sensitivity) with fast glycemia compensation mechanisms. The classical method considers that, over the few hours of fast dynamics, slow variables are essentially constant, and that, seen from the perspective of decades, whatever fast dynamics occur within a given day can be considered as having converged to equilibrium. Transients are however important in their own right, and, due to meals, the subject is never at equilibrium: the necessary mathematical paradigm shift consisted in discarding the near-equilibrium approach and computing instead, at each slow time step, a complete fast dynamics. This approach not only allows us to separately assess quantities (like glycemia at given times of the day or relative to given meal perturbations), which were lumped together with the previous approach; it also offers us the opportunity of simulating diverse perturbation maneuvers at each slow time and follow over the years some of the most commonly used clinical indices of disease.

It should be noticed that “hepatic insulin sensitivity”, as is commonly understood, refers to the effects of insulin on both liver glucose production and liver glucose uptake. In the context of the present discussion, however, this term only refers to insulin-mediated suppression of Hepatic Glucose Output (HGO): the action of insulin in promoting Peripheral Glucose Disposal (PGD) has been termed “peripheral insulin sensitivity” and refers to muscle, adipose tissue and the liver itself. The second major change with respect to the previous version of the model was then the explicit separation of the action of insulin in suppressing HGO from the action of insulin in promoting PGD.

These two major changes are logically related to each other: in order to express the effects of therapeutic regimens, potentially affecting to different degrees the action of insulin on different target tissues, we need the separate representation of hepatic and peripheral insulin sensitivities, as well as the computation of both fasting levels and post-prandial profiles. Daily oscillations of glycemia and insulinemia were felt to carry important information, which is lost when neglecting transient behavior and considering only equilibrium values. Other recent literature points in the same direction [72].

Aging and body changes

Once the physiologic part of the model was finalized, we coupled it with equations that describe naturally varying insulin sensitivity and secretion, as well as the time course of the effect of some typical treatment regimens. Peripheral and hepatic insulin sensitivity, as well as functional insulin secretion (maximal insulin secretory ability per Million β-cells) are here assumed to decrease gradually over time, with rates depending on genetic and lifestyle factors [41, 73]. At the same time, insulin clearance may slightly decrease in the elderly [41], and the natural restoration or healing capacity of the β-cell population (which opposes toxicity to the β-cells from whatever cause) may also decline with age [38, 40]. Toxicity to β-cells is represented in the current model as depending on (hyper-)glycemic values: while toxicity may well be attributed to lipid products [74], we simplified model structure using glycemia as an indicator of possibly several causes of metabolic β-cell impairment.

Therapy mechanisms

In the present work, all DPP treatments were assumed to affect hepatic and/or peripheral insulin sensitivity (λ_GI or k_XGI, respectively). No treatment effect was hypothesized on insulin secretion ability (k_IB) for any of the treatment arms: in fact, it might be postulated that no such direct effect would have been attributed to any of these treatments, unless it was secondary and due to improvement in glucose toxicity.

An initial assumption of identical rates of onset and rates of decay for all effects caused by a given therapy could not be sustained: simulated results based on this assumption could not match the observed data. Conversely, good fits could be obtained by assuming different dynamics for the several effects, which is plausible based on biological and pharmacodynamical principles. In particular, the action of metformin on hepatic insulin sensitivity appeared to be fast and relatively quickly disappearing (whether due to diminished compliance or actual pharmacological loss of effect); instead, metformin appeared to have a slowly appearing and longer lasting action on peripheral insulin sensitivity (in this case indistinguishable from Placebo). This behavior would be plausible if the pharmacologic action of Metformin only affected HGO, while its effect on PGD were due to generic lifestyle changes due to entry into the study).

Conversely, it proved to be reasonable to assume similar rates of onset and decay of effect for both placebo and for the more substantial intensive lifestyle changes: there may be some accumulation of effect over several months as the lifestyle modifications (great or small) produce results, but the positive changes eventually tend to wane as subjects become less compliant.

Some specific conclusions can be drawn by observing what effects should be postulated for each treatment option in order for the corresponding predicted curves to simultaneously approximate all endpoints observed in the DPP study.

The magnitude and time course of the effect on insulin sensitivity produced by the Intensive LifeStyle (ILS) intervention (magnitude that it was necessary to assume in order to adapt model predictions to the observed data) is consistent with the literature: an approximately 20% increase in peripheral insulin sensitivity with ILS is in the range of what has been observed with similar interventions [75–77]. We would expect both hepatic and peripheral insulin sensitivity to improve with weight loss or exercise, with peripheral insulin sensitivity being especially sensitive to exercise and perhaps hepatic sensitivity being at least as responsive to diet. Van der Heijden et al. [78] reported 59% and 23% improvement in peripheral and hepatic insulin sensitivity, respectively, in obese adolescents after initiating exercise. Winnick et al. [79] assessed changes in peripheral and hepatic insulin sensitivity with exercise: they reported a 63% increase in peripheral insulin sensitivity and no significant change in hepatic insulin sensitivity. Malin et al. [80] support the same findings as Winnick. We further note that Ross et al. [75], in a group with diet-induced weight loss of approximately 8% (similar to peak mean weight loss in the DPP ILS group), report a 43% increase in clamp M value. The exercise-induced weight loss group in the same study had a 64% increase in M (with approximately 6% weight loss). Malin et al. [80] also reported a 53% increase in M after an exercise regimen (without weight loss). Interestingly, exercise did not impact hepatic insulin sensitivity in the Malin study (even though this study may have been underpowered because small). Based on this information from the literature, the increase in peripheral insulin sensitivity, which our model would predict as necessary to reproduce the DPP ILS observations, seems to be very realistic.

The calibrated effects of the considered treatments could reproduce simultaneously the time course of the several observed DPP variables (Figs 2 through 6). From these and the corresponding time courses of common derived clinical indices (Figs 8 through 12) a few interesting considerations emerge. First of all, the effect of metformin on HOMA-IR, HOMA-B and insulinogenic index appears very similar to that of Intensive LifeStyle: this is not surprising since these three indices all reflect either only fasting, or fasting and (rather variable) 30 min glycemias and insulinemias. The role of peripheral insulin sensitivity in attenuating post-prandial glycemic excursions is not fully captured by these indices and is evident only when considering 2 hr OGTT values. Another interesting observation concerns the fact that the model predicts a substantial improvement due to metformin (over and above placebo effect) only in low-insulinization, but not in high-insulinization EHC results: this agrees with the interpretation that the low-insulinization clamp does reflect HGO variability / hepatic insulin sensitivity, while at high serum insulin concentrations HGO is effectively suppressed in any case, and the high-insulinization clamp only reflects the current state of peripheral insulin sensitivity.

The trend of all modelled variables points to a continuing progression of the disease. This is plausible and could be explained by any combination of processes: continued loss of β-cells; waning compliance (after initial losses, weight increased over time in the Lifestyle group and medication compliance decreased in both placebo and metformin DPP groups); possible loss of pharmacologic efficacy of metformin over time (consistent with reports by Kahn [66], Brown et al. [68], or Ekstrom et al. [81]).

The predicted rate of change of the β-cell population is relatively small. A thorough analysis of the likely values in-vivo of replication and apoptosis rates, hence of the rapidity of β-cell mass variation, was conducted previously [10]. Given the assessed parameters, the model predicts both a relatively slow increase in β-cell mass following hyperglycemic needs and a relatively slow decline of β-cell mass due to toxicity. While otherwise plausible, this behavior might seem inconsistent with the very rapid expansion of β-cell mass in pregnancy [82], and with its fast decline after delivery: different mechanisms than those here considered may be at play during pregnancy.

We note that the current model can be used to represent the effects of different therapeutic approaches: not only oral hypoglycemic agents of different kinds (metformin, sulfonylureas, SGLT2 inhibitors, thiazolidinediones) and insulin administration, but also bariatric surgery [83–87] or β-cell anti-inflammatory protectors such as Interleukin-1-receptor antagonists [15].

Parameter calibration

Model parameters were calibrated in order for the forecasts to adapt to DPP data: no statistical parameter estimation was conducted. In the present work the goal was not to fit as closely as possible endpoint averages from a given patient sample, obtaining statistical identification of the model for that specific population, but rather to show how reasonable assumptions on mechanisms of action translate into plausible effect parameters, which in turn translate into realistic matches to observed average behavior. While it can be readily seen that the model structure is consistent with available observations (because there exist parameter values that allow the model to capture observed behavior), no diagnostics of goodness-of-fit are therefore calculated.

It is theoretically possible that apparent physiologic insights from the modeling could be highly specific to the current calibration rather than represent generalizable outcomes. It cannot be claimed that the mechanisms underlying the different responses of DPP subjects to the interventions have been identified. What can be said is that the model structure and the chosen parameter values, consistent with reasonable hypotheses on the underlying physiology, produce outputs with are consistent with the actual observations.

It should be appreciated that, in spite of the large number of free parameters, the tight relationship among the state variables, produced by the mechanistic structure of the model, makes it impossible to obtain whatever behavior is desired for each one of several variables simultaneously, particularly when constraining parameters to physiologically acceptable ranges. Indeed, the very fact that the model is actually able to reproduce observed average behavior of several variables simultaneously on the basis of plausible parameter values is indicative of the robustness of the hypotheses made.

It should also be appreciated that some observations have intrinsically large variability: for example, different patients may have widely varying values of insulinemia at 30’ during OGTT since both the absolute peak values as well as the timing of the insulin peaks vary, and the variability of both is reflected in the large variability of the single observation recorded at 30 minutes. It would be pointless to require that the model agrees with observations more closely than the observations agree among themselves.

Model features

The model here proposed is rather complex, and choices had to be made for a reasonably parsimonious representation of the many mechanisms involved. These choices attempted to balance adherence to known physiology and mathematical simplicity. For example, gastrointestinal absorption of glucose was represented with a two-compartment, first-order deterministic process (Eqs 17, 18 and 19). Gastric emptying is in fact irregularly intermittent and pulsatile rather than linear and continuous [88]. The traditional simpler single linear elimination (see the exhaustive discussion in Yokrattanasak et al. [88]) was deemed sufficient and was adopted here as it speeds up simulations considerably. A constant rate of gastric emptying, such as is posited in Eq 16, represents only a crude first approximation to the likely dynamics of the entry of glucose-rich nutrients into the absorbing bowel. As early as Lehmann and Deutsch [89], the rate of gastric emptying was described as trapezoidal for sufficiently large (> 10 g) carbohydrate intake, with ascending, plateau and descending emptying rate phases. More recently, Li et al. [90] investigated functional forms similar respectively to a lognormal distribution or to the right-half of a normal distribution to describe the rate of appearance of glucose (or FFAs) from a mixed meal. On the other hand, Goel et al. [91] suggested that while the linear approach used in Eq 16 is appropriate for liquid meals, it may also be used for mixed meals, albeit with a different (slower) dynamics compared to liquid meals. This issue is worth addressing here because, in general, a poor description of the rate of glucose appearance may confound the effects of other glucose transport, secretion or production mechanisms. In the present case, however, precise identification of the glucose absorption dynamics was not the objective, and a simple, qualitatively plausible dynamics manages to produce intra-day or post-OGTT glycemic variations consistent with the observed trends of 30’ and 120’ post-prandial glycemias over several years. In particular, the simple gastric emptying model should be adequate to simulate delivery of nutrients to the small intestine after the OGTT (small volume, liquid, free of fat and protein), hence to reproduce the indices (post-OGTT glycemias at 30’ and 120’, insulinemia at 30’) tracked over slow time by the overall model.

It is of methodological interest to underscore that our model does not prescribe any explicit set-point, either of glycemia (such as is contemplated for instance in the so-called “minimal model” for the IVGTT) or in the fasting glycemia that beta-cells population dynamics may be assumed to target (as in Ribbing et al. [11]). Instead, what equilibria exist in the system are determined by the free interplay of the dynamics of the different determinants involved. The apparent set-points are mere emergent features of the complex underlying dynamics, and are apt to change, possibly dramatically so, when this dynamics evolves. The parameter values of the model’s underlying dynamics determine the actual observed values of the apparent set-points, such as the equilibrium level of fasting glycemia in the healthy young adult, as well as the values of other observed features, such as the age of development of diabetes in relation with the degree of insulin resistance.

The interplay between progressively developing insulin resistance and eventually failing compensatory pancreatic insulin hypersecretion is widely considered the hallmark of T2DM, but there are different interpretations (possibly corresponding to actual differences in pathophysiologic mechanisms between patient sub-populations) as to the causal chain leading to the eventual decompensation. It has been hypothesized [92] that some acute event (such as a surgical procedure, or a severe infectious episode), determining a sudden increase in insulin resistance, is responsible for the shift from compensation to decompensation, and some mathematical models of the development of T2DM indeed incorporate such an explicit shift [11, 13, 92]. However, over the clinical course of most T2DM patients it is not possible to identify such a triggering event. The present model makes no recourse to external triggering events and does not need the introduction of an explicit regime shift (a sudden change in parameter values, the introduction of a new external forcing function) in order to reproduce a rapid worsening of the clinical conditions at some point in the life of the subject. Instead, our model formalizes the concept that a persistent hyperglycemic insult, determined by long-standing, possibly progressive degrees of insulin resistance, brings about a progressive decline of insulin sensitivity, and that the even mild glucose toxicity connected with persistent insulin resistance eventually damages pancreatic replication reserve, determines an eventual decline of β-cell mass and an eventual failure of compensating insulin hypersecretion, resulting finally in rapid acceleration of hyperglycemia and in the overt clinical picture of frank T2DM.

The DPM15 model does not appear to support the hypothesis that primary insulin hypersecretion might be the causal factor of the development of T2DM. This hypothesis was discussed by Corkey [93] in hypothetical pathophysiological terms and by Goel [94] in mathematical terms, through an adaptation of the original Topp model to include a direct effect of insulin on β-cell dynamics. Simulations (not shown) have been conducted with our model assuming no primary insulin resistance (neither peripheral nor hepatic) and assuming conversely that the long-time behavior of glucose-driven insulin secretion by β-cells, instead of mildly decreasing with age, actually doubles. In this way the model expresses a progressive insulin hypersecretion, given the same prevailing glucose concentrations. All such simulations fail to reproduce a progressively increasing glycemia up to diabetic levels. The present model has not been modified to include a direct effect of insulin on β-cells, also because we agree with the statement [94] that the required causal signal linking insulin resistance and insulin hypersecretion may well be glycemia itself. In the present model the apparent Corkey paradox is replicated (glycemia not immediately rising upon worsening of insulin resistance), due to the controlling effects of increased insulin secretion up to the point where relative endocrine pancreatic insufficiency develops. Insulin hypersecretion could in fact be the mechanism by which glucose toxicity exerts its detrimental effects on β-cell turnover. Under this hypothesis, while hyperinsulinemia would not by itself be directly toxic, the glucose-induced, prolonged hypersecretory state would be damaging to the β-cell. This interpretation is consistent with the role of ‘ER stress’ that has been proposed for β-cell death: the overly taxed secretory apparatus of the cell results in accumulation of misfolded proteins which, in turn, trigger inflammatory and apoptotic mechanisms. In all these cases, from a modelling point of view we may take hyperglycemia as a fair indicator of the toxic situation for the β-cell population.

Limitations of the current work

One potential limitation is that the current model does not take into account the possible effect of glucose toxicity on insulin secretory function by existing β-cells. The effect of sustained, moderate hyperglycemia is here exerted exclusively on β-cell replication, hence on the maintenance of β-cell mass: should there be reasons to assume that a toxic effect is also exerted on insulin secretory mechanisms, this ought to be incorporated in the model. This is a moot point however, because current literature reports both increased and decreased insulin secretion with acute hyperglycemia [95].

Another limitation consists in not taking into account lipid metabolism, variations in fat mass, body size etc. as measurable indicators of peripheral insulin resistance and as possible contributing factor to low-key, systemic, continuous inflammation adding its toxic effect on β-cell replication.

A third limitation of the current model is its simplistic depiction of (mono-exponential) gastric emptying. Even so, the whole aggregated model captures relevant glycemic oscillations throughout the day, but future work will consider the introduction of a stochastic gastric emptying sub-model as well as the use of stochastically variable meal composition and size.

Finally, model validation against independent sets of observations is clearly desirable and will need to be addressed in the future, similarly to what was done for the previous version of the model in comparing its output with the CANOE study results. [18, 96]

Conclusion

A new mathematical model of the long-term development of Type 2 Diabetes Mellitus is consistent with available literature, is able to reproduce experimentally observed effects of therapeutic interventions on several endpoints, including fasting and post-prandial glycemias and insulinemias, and can simulate the evolution of common clinical and experimental indices in cohorts of virtual patients. Such an elaborate model will need to be further validated: as it stands now, it incorporates plausible physiology, agrees with available observations, and allows the investigator to formulate quantitative, testable questions for relevant patient populations.