A Novel Mathematical Model Detecting Early Individual Changes of Insulin Resistance

Abstract

Background:

Insulin resistance (IR) and hyperinsulinemia as well as obesity play a key role in the metabolic syndrome (MetS), type 2 diabetes (T2D), and associated cardiovascular disease. Unfortunately, IR and hyperinsulinemia are often diagnosed late (i.e., when the MetS is already clinically evident). An earlier diagnosis of IR would be desirable to reduce its clinical consequences, in particular in view of the increasing prevalence of obesity and diabetes conditions. For this purpose, we developed a mathematical model capable of detecting early onset of IR through small variations of insulin sensitivity, glucose effectiveness, and first- or second-phase responses.

Materials and Methods:

Murine models provide controlled conditions to study various stages of IR. Various degrees of hypercholesterolemia, obesity, IR, and atherosclerosis were induced in low-density lipoprotein receptor-deficient mice by feeding them cholesterol- or sucrose-rich diets. IR was assessed by oral glucose tolerance tests. Controls included animals fed exclusively, or switched back to, regular chow. A nonlinear mathematical model of the order of 5 was developed by refining Bergman's “Minimal Model” and then applied to experimental data.

Results:

Different metabolic constellations consistently corresponded to specific and close-meshed changes in model parameters. Reduced second-phase glucose sensitivity characterized an early impaired glucose tolerance. Later stages showed an increased first-phase glucose sensitivity compensating for decreased insulin sensitivity. Finally, T2D was associated with both first- and second-phase sensitivities close to zero.

Conclusions:

The new mathematical model detected various insulin-sensitive or -resistant metabolic stages of IR. It can therefore be implemented for quantitative metabolic risk assessment and may be of therapeutic value by anticipating the start of therapeutic interventions.

Introduction

I nsulin resistance (IR), hyperinsulinemia, dyslipidemia, and obesity play key pathophysiological roles in the metabolic syndrome (MetS) and type 2 diabetes (T2D), which in turn are associated with a markedly increased risk of cardiovascular diseases (CVD). Unfortunately, IR and hyperinsulinemia are often diagnosed late (i.e., when the MetS is already clinically advanced and accompanied by metabolic co-morbidities and CVD). An earlier diagnosis of IR could reduce clinical consequences by prompting earlier treatment but would require a diagnostic method more sensitive to early stages. For this purpose, we developed a mathematical model capable of detecting early onset of IR through small variations of insulin sensitivity, glucose effectiveness, and first- or second-phase responses.

Depending on ethnicity, genetic, and lifestyle risks, 30–40% of the population of developed countries develops MetS.¹ Although the importance of glucose intolerance, dyslipidemia, obesity, and hypertension is well recognized, various definitions of the MetS exist.^1

–6 For example, the World Health Organization and the European Group for the Study of Insulin Resistance include IR among defining central criteria; in contrast, the International Diabetes Foundation uses central obesity as a diagnostic criterion.^1

–6

IR and hyperinsulinemia are strongly linked to an early onset of T2D, obesity, dyslipidemia, and endothelial dysfunction, as well as CVD.^3,7 Initially, an insufficient effect of insulin, resulting from either insufficient insulin secretion or insufficient insulin effect at the cellular level, is met by a compensatory increase in insulin secretion, and plasma glucose levels remain in the euglycemic range (Figure 1a). Prolonged excessive demand for insulin, associated, for example, with obesity, often leads to decreasing function of pancreatic β-cells, increasing hyperglycemia, and eventually manifest T2D. The diagnostic criteria for the final stages of IR and T2D in humans are well established.^5,6

FIG. 1.

Development of insulin resistance, overview of murine studies, and optimization of the mathematical model. (a) Stages of insulin resistance and type 2 diabetes (adapted from Fonseca³). Increasing pancreatic insulin production initially compensates for increasing insulin resistance. Once the capacity for compensation is exhausted or β-cell damage decreases insulin secretion, glucose levels gradually increase until manifest diabetes is reached. (b) Experimental design of murine studies. Groups of low-density lipoprotein receptor^−/− mice were subjected to an oral glucose tolerance test (OGTT) before (Group A), during (Group B), and after (Group C) being fed an obesogenic diet and compared with controls (Groups D and E). (c) Optimization of the mathematical model by parameter identification based on in vivo data. The deviation of the experimental G _exp and I _exp and the model outputs G and I is assessed by the objective function J. An optimization algorithm minimizes J with respect to the parameter vector p of the mathematical model. Color graphics available online at www.liebertonline.com/dia

Dyslipidemic conditions, especially increased levels of free fatty acids, also contribute to IR. Free fatty acids result mainly from adipose tissues because of the lipolysis of triglyceride-enriched lipoproteins.^1

–4,7 Conversely, increasing IR generates more free fatty acids, mainly by increasing lipolysis, and thus worsens dyslipidemic conditions.^1

–4,7 Therefore, we used an in vivo model susceptible to both diet-induced IR and dyslipidemia.

Because IR is initially clinically silent, it may remain undetected while already exerting pathogenic effects relevant for diabetes and CVD. The main difficulty in establishing a mathematical model capable of detecting the onset of IR lies in the lack of well-defined physiological correlates of early-stage IR that might be used for a noninvasive diagnosis.

In both humans and mice, insulin glucose homeostasis is usually first assessed by measurements of fasting or fed plasma glucose and insulin respectively C-peptide. A widely used diagnostic follow-up procedure is the oral glucose tolerance test (OGTT), which measures glucose and insulin responses to a defined glucose bolus administered orally. Although OGTTs can detect reduced insulin secretion, conclusive evidence for IR (as opposed to pancreatic β-cell damage) requires hyperinsulinemic euglycemic clamps.^8,9 Other assays include insulin tolerance tests and insulin suppression tests. Many of these are invasive. Mathematical models permitting extraction of additional information on the onset of IR from OGTTs would therefore be of value, in particular if they predict later T2D.

We here report the development of a novel mathematical model of insulin glucose homeostasis based on Bergman's “Minimal Model,” which provides measures of glucose effectiveness 1/T _G and insulin sensitivity k_X . To improve modeling of the physiological first-phase insulin response, we introduced the respective first- and second-phase parameters k _G1 and k _G2 into the equation and carried out data analysis in two sequential steps. The new model was first validated with the help of a human reference dataset based on intravenous glucose tolerance tests (IVGTTs).¹⁰ It was then applied to data obtained by performing OGTTs on five groups of low-density lipoprotein receptor (LDLR)-deficient (LDLR ^−/−) mice, an established model of hypercholesterolemia, insulin resistance, diet-induced MetS, and atherosclerosis.^{11

–20} In these murine groups, various degrees of IR and stages of the MetS were induced by feeding them either regular rodent chow or an obesogenic high-fat, high-cholesterol diet for various periods. We also switched one group back to regular chow, to mimic the return to a healthy diet after onset of the MetS. We thus induced various metabolic constellations resulting from the underlying genetic defect of LDLR ^−/− mice or the various dietary regimens. We then characterized these constellations by conventional clinical parameters and by using our mathematical model, in the hope that the latter would indicate early stages of the MetS and thus permit early diagnosis of the pathogenic process and earlier therapy. Preliminary data from this study have been presented previously in abstract form.²¹

Materials and Methods

Murine model of the MetS

In order to mimic different stages of human IR, an established murine model of atherosclerosis and T2D was used.^{11

–20} Five groups of age-matched male LDLR ^−/− mice (n=43) were fed normal rodent chow for 74 days, at which time an OGTT was performed on Group A. Group B was fed an established obesogenic and atherogenic diet containing 1.25% cholesterol and 21% fat (Teklad TD96121; Harlan, Indianapolis, IN) for 109 days, and an OGTT was performed after 78 days. Group C followed the same regimen but was then switched back to regular chow to detect the benefit of a low-risk diet after the onset of IR. Groups D and E were controls fed regular chow. A schematic overview of experimental groups is provided in Figure 1b.

OGTTs were performed after 16 h of fasting. Twenty-five microliters of blood was collected from the retroorbital plexus of isofluorane-anesthetized mice prior to as well as 15, 30, 60, 90, and 120 min after gavaging with a sterile 10% glucose solution (1.5 mL/kg of body weight). Plasma glucose levels were determined using glucose test strips (Devine Medical Supplies, Santa Fe Springs, CA) and a glucometer (MediSense; Abbott, Chicago, IL) four times a day. Plasma insulin level was determined by an ultrasensitive mouse enzyme-linked immunosorbent assay (catalog number 10-1150-01; Mercodia, Uppsala, Sweden).

To assess metabolic changes by the criteria used in humans, the homeostasis assessment model of IR (HOMA-IR) index was calculated from fasting insulin I ₀ and fasting glucose G ₀: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*}{ \rm HOMA - IR } = \frac { I_0 \cdot G_0 } {405 \frac {\mu \rm U} {\rm mL } \frac {\rm mg } { \rm dL }} \end{align*} \end{document}

Mathematical model of glucose–insulin homeostasis

Any mathematical model of the glucose–insulin homeostasis must reflect the complex metabolic changes observed in vivo. One of the first models of the glucose insulin homeostasis was the “VT-model,”²² also known as the “Minimal Model.”²³ It provides established parameters with a physiological meaning, such as glucose efficiency and insulin sensitivity. We previously used the Minimal Model^24
–26 and now expanded it as follows.

The Minimal Model considers three compartments: the plasma glucose G(t), the plasma insulin I(t), and the interstitial insulin X(t), which summarizes insulin in non-plasma tissues and is also called the remote insulin compartment according to Bergman et al.²³ The change of plasma glucose G(t) is driven by an internal effectiveness (with time constant T _G) and a nonlinear expression that is controlled by the current interstitial insulin X(t): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \begin{split} \frac {dG ( t )} {dt} = - \frac {1} {T_G} \Delta G ( t ) - k_X X ( t ) G ( t ) , \\ \Delta G ( t ) = G ( t ) - G_{{ \rm basal}} , \quad G ( 0 ) = G_0\end{split} \tag{1}\end{align*} \end{document}

The reciprocal time constant 1/T _G is often interpreted as “glucose effectiveness,” and the gain k_X is known as “insulin sensitivity” of the glucose system.²⁷ To simulate the dynamic situation of an OGTT, we started the simulation at t=0, when exogenously administered glucose reaches the plasma. Hence, in our model an initial plasma glucose concentration G ₀ is assumed that is already elevated compared with its basal level: G ₀>G _basal. Note that this is not identical to the glucose measurement at “T0” of an OGTT, which denotes the basal equilibrium before glucose administration.

The change of plasma insulin I(t) is modeled in a similar way. An internal clearance (with time constant T _I) and an external input v(t) are distinguished: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \frac { dI ( t ) } { dt } = - \frac { 1 } { T_I } \Delta I ( t ) + \frac { 1 } { T_I } v ( t ) , \quad \Delta I ( t ) = I ( t ) - I_ { { \rm basal } } , \quad I ( 0 ) = I_0 \tag { 2 } \end{align*} \end{document}

Here, the basal insulin concentration is used as initial condition I ₀=I _basal because no external insulin excitation is applied, and the insulin responses are presumed to lag behind the initial glucose increase. The external input v(t) is considered to be linear. The input v(t) must be controlled by the plasma glucose G(t). Bergman et al.^22,23 suggested the time-variant function v(t)=max[0, γt(G(t)−G _basal)]. From our point of view, this approach has two disadvantages: an early plasma insulin peak (for small times t) cannot be modeled, but it is observed in OGTT as a consequence of an increased plasma glucose level. Bergman et al.^22,23 solve this by assuming an initial insulin concentration I ₀>I _basal, but this does not correspond to the test conditions because no initial insulin is applied to the subject. Second, there is no physiological motivation for an increasing insulin reaction just because the time t increases.

To make the model more consistent with actual physiological events, we suggest differentiating between a first- and second-phase insulin reaction. The first-phase reaction is controlled by the time derivative d(ΔG(t))/dt of plasma glucose in combination with a lag element (small time constant T ₁). The second-phase reaction is controlled by ΔG(t) itself again in combination with a lag element (larger time constant T ₂>T ₁). Both parts are weighed with gains k _G1 and k _G2 and summed up in v(t). The frequency domain transfer function is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \begin{split}& V ( s ) = \left( \frac {k_{{\rm G}1^S}} {1 + T_1s} + \frac{k_{{\rm G}2}} {1 + T_2s} \right) \cdot \\ & \Delta G ( s ) = \frac{k_{{\rm G}1}} {T_1} \cdot \Delta G ( s ) - \underbrace{\frac{\frac{k_{{\rm G}1}} {T_1}} {1 + T_1s} \cdot \Delta G ( s ) }_{V_1 ( s ) } \underbrace { + \frac {k_{{\rm G}2}} {1 + T_2s} \cdot \Delta G ( s )}_{V_2 ( s )} \end{split} \tag{3}\end{align*} \end{document}

Hence, we introduce k _G1 and k _G2 as “first-phase” and “second-phase glucose sensitivity” of the insulin system, respectively. This simplifies the model approach of Grodsky,²⁸ who distinguished a labile (fast) and a stable (slow) insulin compartment with a threshold for insulin release. It can also be considered to be a model of the controller that regulates the insulin secretion as a function of the current plasma glucose. The use of a proportional (P) and a time derivative (D) part find its equivalent in technical PD controllers.²⁹

The second-order transfer function in Eq. 3 is transformed into two first-order differential equations (Eqs. 4 and 5) in the time domain. Among the infinite number of state–space representations we obtain the Jordan Canonical Form by introducing the states v ₁(t) and v ₂(t) according to Eq. (3). The output equation (Eq. 6) includes a feed-through term: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \frac { dv_1 ( t ) } { dt } = - \frac { 1 } { T_1 } v_1 ( t ) - \frac { k_ { { \rm G } 1 } } { T_1^2 } \Delta G ( t ) , \quad v_1 ( 0 ) = 0 \tag { 4 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \frac { dv_2 ( t ) } { dt } = - \frac { 1 } { T_2 } v_2 ( t ) + \frac { k_ { { \rm G } 2 } } { T_2 } \Delta G ( t ) , \quad v_2 ( 0 ) = 0 \tag { 5 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*}v ( t ) = v_1 ( t ) + v_2 ( t ) + \frac { k_ { { \rm G } 1 } } { T_1 } \Delta G ( t ) \tag { 6 } \end{align*} \end{document}

The third compartment for the interstitial insulin X(t) is modeled in line with Bergman et al.^22,23 As in Eq. 2, there is an internal clearance (time constant T_X ) and a linear excitation that is controlled by the insulin I(t): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \frac { dX ( t ) } { dt } = - \frac { 1 } { T_X } X ( t ) + \frac { 1 } { T_X } \Delta I ( t ) , \quad \Delta I ( t ) = I ( t ) - I_ { { \rm basal } } , \quad X ( 0 ) = 0 \tag { 7 } \end{align*} \end{document}

Because we assign the insulin sensitivity k_X to Eq. 1, the state X(t) has the dimension of an insulin concentration, which seems to be more appropriate than the dimension of min⁻¹ used by Bergman et al.^22,23

In order to obtain a state–space model²⁹ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \frac { d { \bf x } ( t ) } { dt } = { \bf f ( x ) } \tag { 8 } \end{align*} \end{document}

Equation 6 is inserted into Eq. 2, and Eqs. 1, 2, 4, 5, and 7 are summarized to a vector differential equation of the order n=5 with the state vector x(t)=[G(t), I(t), v ₁ (t), v ₂ (t), X (t)]^T.

Parameter identification

For parameter identification, the mathematical model is typically evaluated over the time horizon of an OGTT, which is usually 120 min. As the initial condition, a glucose excitation G ₀>G _basal and a basal insulin concentration I ₀=I _basal are assumed. As illustrated in Figure 1c, the model outputs for plasma glucose level G(k) and plasma insulin level I(k) are compared with the corresponding in vivo measurements G _exp(k) and I _exp(k). The deviations are weighed with w _G and w _I and assessed within an objective function that sums up the square errors over the measurements \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$k = 1 , \ldots , N$$ \end{document} : \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} J(\textbf{p}) = \begin{cases} \sqrt{ \frac {1} {N} \sum\limits_{k= 1}^N [W_G (G (k , \textbf{p}) - G_{\exp} (k))^2 + w_I (I (k , \textbf{p}) - I_{\exp} (k))^2]} \quad {\rm if} \ \textbf{p} \ {\rm is} \ {\rm valid} \\ \infty \qquad {\rm if} \ \textbf{p} \ {\rm is \ not \ valid} \end{cases} \tag {9} \end{align*} \end{document}

We try to reach a balanced weighting and a dimensionless objective function J by choosing w _G=(1 mg/dL)⁻² and w _I=(0.667 μU/mL)⁻². If no measurement is available, the corresponding w _G or w _I is set to zero. In order to consider parameter constraints, J is set to infinity in case of forbidden parameter combinations.

An optimization algorithm has to find an optimal set of model parameters p that minimizes the objective function J(p) (see Fig. 1c). For that purpose we apply the heuristic simplex search algorithm from Nelder and Mead.³⁰ The iterative algorithm stops when the initial simplex volume could be reduced by the factor 10⁻⁶ or when the maximum number of 150 iteration steps has been reached. Model and optimization algorithm are implemented in MATLAB (MathWorks, Natick, MA).

Validation of the model with human data

As a first step to quantify the model, we followed the program of Pacini and Bergman.¹⁰ Their values determined from a frequently sampled IVGTT were converted to define the time constants T _G, T _I, and T_X and the gain k_X . As basal values for G _basal and I _basal, the last measurements were taken (Table 1). In a second step, the respective time constants T ₁ and T ₂ for the first- and second-phase insulin response were defined in line with the step response in Grodsky,²⁸ and an optimization was carried out in order to identify the corresponding gains k _G1 and k _G2 (Table 1). In vivo measurements and in silico results of the IVGTT are shown in Figure 2. The evaluation of the objective function from Eq. 9 resulted in J=24.76, which corresponds to a maximum root mean squared error of 24.76 mg/dL for plasma glucose G(t) or of 16.51 μU/mL for plasma insulin I(t). The model obtained serves as a human reference model in the following analyses.

FIG. 2.

Comparison of in vivo measurements of a frequently sampled intravenous gluocse tolerance test from Pacini and Bergman¹⁰ and its corresponding in silico simulation of our model. Color graphics available online at www.liebertonline.com/dia

Table 1.

Parameters of the Human Reference Model

Parameter	From (reference)	Value	Units
T _G	1/p ₁=(Pacini and Bergman¹⁰)	32.45	min
G _basal	G(t _final)=(Pacini and Bergman¹⁰)	90.0	mg/dL
k_X	p ₃/p ₂=(Pacini and Bergman¹⁰)	507×10⁻⁶	ml/μU mg/dL/min
T _I	1/n=(Pacini and Bergman¹⁰)	3.33	min
I _basal	I(t _final)=(Pacini and Bergman¹⁰)	7.0	μU/mL
T ₁	Similar to Grodsky²⁸	4.0	min
k _G1	Optimization	3.601	μU/mL dL/mg min
T ₂	Similar to Grodsky²⁸	12.0	min
k _G2	Optimization	0.356	μU/mL dL/mg
T_X	1/p ₂=(Pacini and Bergman¹⁰)	47.78	min

Results

Murine model

Body weights, plasma cholesterol, and plasma triglyceride levels before, during, and after the obesogenic diet are summarized in Figure 3a–c. As expected, endogenous and dietary risk factors could be linked to unbalanced metabolic conditions, such as altered plasma cholesterol and triglyceride levels (compare Group B vs. Groups A and D in Fig. 3b and c). The change back to normal chow resulted in declining plasma lipids in Groups C and E but no significant changes in body weights (Fig. 3a–c).

FIG. 3.

(a–f) Metabolic characterization of the murine model and (g–i) results of the in silico experiments. Significances were determined by two-sided t test for unpaired groups. Groups of low-density lipoprotein receptor^−/− mice were subjected to an oral glucose tolerance test (OGTT) before (Group A), during (Group B), and after (Group C) being fed an obesogenic diet and compared with controls (Groups D and E). HOMA-IR, homeostasis assessment model of insulin resistance. Color graphics available online at www.liebertonline.com/dia

Fasting measurements are displayed in Figure 3d and e. OGTT measurements are shown in Figure 4. Several phenotypes with different degrees of IR were evident. Compared with Groups A and D, Group B showed significantly increased fasting insulin levels but no increase in fasting glucose levels. This indicates that insulin levels required to maintain euglycemia were higher than those on regular chow and that Group B had IR but did not suffer from a manifest diabetes. With the return of insulin-resistant obese mice to regular chow in Group C as well as an increasing age and persisting weight gain, insulin levels were reduced significantly compared with Group B, but glucose levels continued to increase. Thus, β-cell damage progressed irreversibly over time, and a persistent degree of IR and an onset of T2D characterize Group C.

FIG. 4.

Comparison between the actual oral glucose tolerance test data of murine Groups A–E and the predictions by the mathematical model. Plasma glucose G(t) and insulin I(t) were measured before glucose administration and 15, 30, 60, 90, and 120 min thereafter. Data of individual subjects are shown in gray, and group means with SD error bars in black. Additionally, the corresponding in silico predictions are plotted. The solid line denotes the group mean; dashed lines denote the corresponding SD. Color graphics available online at www.liebertonline.com/dia

Note, however, that data of the control groups (A, D, and E) do not represent physiological levels in wild-type mice but strongly reflect the underlying genetic defect resulting in spontaneous hyperlipidemia. Furthermore, animal and human glucose responses cannot be compared directly. The same is true for the threshold for manifest T2D. By human standards specified in the guidelines of the World Health Organization,⁶ Group C comes closest to the human threshold for T2D. Increased HOMA-IR also indicated IR in Groups B and C³ but could not distinguish the different underlying mechanisms. For this purpose, dynamic reactions such as those retrieved from an OGTT have to be considered, which can be done by a dynamic model.

Application of the novel mathematical model to the murine data

We then applied our mathematical model to see whether it was sensitive to early stages of IR. From our point of view the time constant T _G of glucose effectiveness and the insulin sensitivity k_X , as well as the glucose sensitivities k _G1 and k _G2 are of particular diagnostic interest. Furthermore, the initial plasma glucose level G ₀ is identified in order to allow the model to adapt to the different form of initial glucose administration of an OGTT.

In the in silico simulations based on the murine data we observed that T _G and k_X were closely associated. Identification of both constants in one model may therefore lead to inconsistent results. For this reason, we first determined the parameter set {T _G,k _G1,k _G2,G ₀} in an initial analysis. As Table 2 shows, T _G is growing to large values for all control groups (A, D, and E) that are only limited by the maximum number of iterations of the optimization algorithm. Therefore, T _G can be assigned to an arbitrary large value. We chose T _G=1,000 min for simplicity and optimized {k_X ,k _G1,k _G2,G ₀} in the main analysis. Table 3 indicates the objective functions obtained for all groups in the initial and the main analysis. As expected, the objective function J performed better in the main analysis (i.e., yielded smaller values than in the initial analysis for all groups). It is interesting that the in silico experiments achieved smaller values of J for the groups not subjected to the extreme high-fat diet (A, D, and E), which means that the model characterization is better in these cases.

Table 2.

Results of the Initial Analysis

Group	Parameter T_G (in min)	Parameter k_G1 (in μU/mL dL/mg min)	Parameter k_G2 (in μU/mL dL/mg)
As in Table 1	32.4	3.563	0.348
A	7,298.2±11,274.3^a	2.807±1.953	0.045±0.060
D	8,908.1±16,915.8^a	2.874±1.920	0.052±0.043
E	9,972.3±10,879.0^a	4.594±1.131	0.103±0.089

Data are mean±SD values for model parameters T _G, k _G1, and k _G2 for Groups A, D, and E. Differences from the human reference model are observed for T _G and k _G2, whereas k _G1 tends to be in the same order.

Simplex optimization does not converge and is stopped after reaching the maximum of 150 iteration steps.

Table 3.

Resulting Objective Function J of the Initial and Main Analyses

Group	Initial analysis	Main analysis
A	30.80±18.72	17.92±11.44
B	85.64±35.66	46.12±20.41
C	56.54±29.72	38.49±18.45
D	25.64±14.09	14.21±6.66
E	24.74±11.73	14.80±6.96

Data are mean±SD values. A small value J represents a good model representation (see Eq. 9 in Materials and Methods).

In order to assess the mathematical model, its in silico predictions of the murine OGTTs were compared with the actual in vivo measurements (Fig. 4). A good agreement was seen for Groups A, D, and E. The markedly impaired glucose tolerance (IGT) of Groups B and C resulted in high insulin levels for t<20 min and larger SDs.

Model-based differentiation among different stages of the MetS

IGT associated with LDLR deficiency

We then determined whether metabolic changes meeting the criteria for IGT can be detected in the experimental groups not exposed to the high-fat diet (i.e., Groups A, D, and E). The HOMA-IR index was greater than 5 in Group A, which suggests IR, but well below 5 in Group D. The mathematical model provided more consistent information. Table 2 indicates the parameters identified in mice and compares them with those of the human reference model. Within the obvious limitations of such a comparison, the following conclusions are supported:

• The first-phase glucose sensitivity k _G1 is in the same order as in the human reference model. Hence, the first-phase insulin response seems to be similar.

• In contrast, the second-phase glucose sensitivity k _G2 is reduced by a factor of 3–8. This may indicate increased glucose resistance in the second phase.

• During parameter identification the estimate of T _G is growing continuously to very large values until the iteration is aborted. If T _G is striving to infinity, the glucose efficiency 1/T _G is going to zero. It seems to be a common characteristic of the LDLR ^−/− mice that their ability to clear glucose (without insulin stimulation) is disrupted.

Together, initial analysis of data with the novel model indicates an early diabetes phenotype charactecterized by a reduced second-phase insulin response and a marginal self-contained glucose clearance in LDLR ^−/− mice.

Combined effect of LDLR deficiency and obesogenic high-cholesterol high-fat diet

In Group B, fasting glucose levels were normal, but fasting insulin levels were markedly increased, most likely in order to compensate for the rapid onset of IR.

Again, the mathematical model may provide an explanation. Figure 3g–i shows the parameters identified by the main analysis of all five groups. Compared with the previously discussed Groups A, D, and E, the insulin sensitivity k_X is decreased in Group B. At the same time the glucose sensitivities k _G1 and k _G2 are increased, most likely in order to compensate for the loss of insulin sensitivity. For example, in relation to Group A the gain ratios are \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$k_X^{ ( B ) } / k_X^{ ( A ) } = 0.237$$ \end{document} on the one hand and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$k_{{\rm G}1}^{ ( B ) } / k_{{\rm G}1}^{ ( A ) } = 4.859$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$k_{{\rm G}2}^{ ( B ) } / k_{{\rm G}2}^{ ( A ) } = 6.602$$ \end{document} on the other hand. However, the total loop gain only slightly differs by 0.237×4.859=1.152 for the first-phase response and 0.237×6.602=1.565 for the second-phase response. The model clearly differentiates between the metabolic conditions induced by the genetic defect only and or in combination with diet.

T2D

The mean fasting glucose level of Group C in Figure 3d was 143.0 mg/dL. Thus, by human standards, these animals had a manifest T2D. The in silico model yields exactly the same diagnosis: the glucose sensitivities k _G1 and k _G2 are significantly lower than those of Group B, whereas the insulin sensitivity is not able to compensate this.

Compared with the IGT group (Group A), it is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$k_X^{ ( C ) } / k_X^{ ( A ) } = 0.381$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$k_{{\rm G}1}^{ ( C ) } / k_{{\rm G}1}^{ ( A ) } = 0.421$$ \end{document} , such that the total loop gain is 0.381×0.421=0.160. For parameter k _G2 it is 0.381×1.041=0.397, respectively. The total loop gain has dropped significantly in both cases below 1.0, which can be explained by β-cell dysfunction.

In summary, the new mathematical model allows detections of different stages of IR and is sensitive to its early progression.

Discussion

Insights from the mathematical model

We here introduce a novel model designed to be sensitive to early patterns of the MetS and thus to facilitate early diagnosis of IR. Compared with previous methods, reviewed for example by Matsuda,³¹ our model not only provides a measure of insulin sensitivity,³² but also determines additional parameters, which allow us to distinguish between different close-meshed patterns of IR.

The key requirement of our mathematical model was to accurately reproduce the glucose–insulin homeostasis observed in vivo. It was derived from the ordinary differential equation type VT (or Minimal) Model^22,23 by adding two compartments. Such additions have been successfully performed in order to consider external influences (e.g., diets³³), to reflect more complex mechanisms,³⁴ or to distinguish different organs,³⁵ but they usually come at the price of requiring additional measurements. Our model has the advantage of being time-invariant and of differentiating between first- and second-phase insulin responses. For a broad review of model-based methods in the field of diabetes see the review of Cobelli et al.²⁷; different modeling approaches are compared by Makroglou et al.³⁶ and Mari.³⁷

After successful validation of the model with IVGTT data from a human reference population, we applied the novel model to different age- and genetically matched groups of LDLR ^−/− mice, which were fed either with standard chow for various time periods (dyslipidemia reflecting only the genetic defect) or with a standardized high-fat/high-cholesterol diet (combined effect of genetic defect and high-risk lifestyle). As shown in Table 4, four different stages of IR could be distinguished with the help of the mathematical model. These are in line with the prevailing hypothesis of a gradual process beginning with IGT and ending in severe T2D. The model allowed us to separate early stages of IR characterized by increasing insulin production from later stages showing decreased insulin levels and rising glucose levels and from T2D. It therefore appears that the model parameters {T _G, k_X , k _G1, k _G2} provide the means for an early diagnosis of IR in conjunction with a simple OGTT.

Table 4.

Stages of Murine Insulin Resistance Detected In Vivo and by In Silico Simulation

Pattern	Mechanisms	In vivo	In silico
Normal (comparable to the human reference group)	Physiological response to glucose challenge	Physiological plasma glucose and insulin levels	Parameters {T _G, k_X , k _G1, k _G2} as shown in Table 1
IGT: genetic dyslipidemia, increasing age (Group A)	Increasing IR (decreasing glucose effectiveness); some compensatory insulin increase	Raised HOMA-IR index, delayed insulin response	Reduced glucose effectiveness 1/T _G, reduced second-phase glucose sensitivity k _G2
IGT plus high-risk lifestyle: high-fat diet (Group B)	Compensatory insulin increase; beginning insulin deficiency	Increased fasting plasma insulin or increased 2-h plasma glucose	Reduced insulin sensitivity k_X , increased glucose sensitivity k _G1, but normal or slightly increased loop gain (k_X ×k _G1)
T2D (Group C)	Insulin production insufficient to maintain euglycemia (IR and β-cell damage)	Increased fasting glucose	Glucose sensitivity k _G1 collapses, reduced loop gain (k_X ×k _G1)

HOMA-IR, homeostasis assessment model of insulin resistance; IGT, impaired glucose tolerance; IR, insulin resistance; T2D, type 2 diabetes.

Provided that further validation in human subjects supports these findings, the method may offer practical and cost advantages. More important is that the method may provide an additional level of certainty in diagnosing early stages of IR and thus be of value for the prevention of the clinical consequences of T2D, mainly CVD.

With regard to basic science applications in mice, the emphasis clearly lies on early detection. Currently, there is no clear consensus about what constitutes borderline hyperinsulinemia and hyperglycemia in the many genetically modified mouse models available. Nor is there an agreement on basic aspects of euglycemic clamps, such as the duration of fasting, the use of an initial insulin bolus, or the subsequent constant insulin infusion rate.⁸ Technical aspects further complicate the comparison of data obtained by various groups. For example, no uniformly recognized standards exist for the time interval between catheter implantation, anesthetic agents, and the blood sampling technique from arterial catheters or peripheral sites. Glucose clamp studies in mice therefore often represent far more stressful conditions than those encountered in humans. In contrast, OGTTs are carried out under fairly standardized conditions even in mice, except for the duration of fasting. Our model, being sensitive to early changes of IR, should therefore be of great value. However, the very absence of uniform criteria for IR in mice, as well as the absence of histologic or other verifiable changes associated with the onset of IR, makes it difficult to establish the physiological significance—and relevance—of the new parameters introduced by our model.

Novel diagnostic indices

The HOMA-IR index is widely used for human screening purposes but was not able to differentiate stages of murine IR, probably because it is based on fasting measurements. Indices based on responses triggered, for example, by oral glucose, like in the OGTT provide more information. Stumvoll et al.³⁸ investigated several indices for IR, termed MCR, ISI, 1^st PH, and 2^nd PH, that are calculated as linear combinations of the (nonfasting) measurements of an OGTT. All of these indices showed good agreement with clamp measurements. They are, however, not obtained from a dynamic model (defined by differential equations), and their model coefficients could not be interpreted in the scope of a physiological model. Furthermore, these indices also utilize the body mass index, which hinders comparison with animal models. In contrast, Pillonetto et al.³⁹ provided an analytical solution for a truly dynamic index of insulin sensitivity, based on Bergman's model. Compared with this, our approach requires an optimization procedure to obtain the solution but determines a set of parameters. For the classification of different patterns of IR proposed in Table 4, a mere scalar index will be insufficient. For this purpose, a multidimensional set of parameters has to be provided, based on a model such as the one presented here or that proposed by Stumvoll et al.³⁸ The combination of in vivo data and mathematical procedures offers far more detailed pathophysiological insights into the glucose–insulin homeostasis and therefore permits an early diagnosis and greater therapeutic options.

Today, the computational resources needed for model evaluation and optimization no longer impede its widespread application.^24,25 Altogether, our model offers a novel perspective to diagnose and target IGT and early evolution of IR.

Footnotes

Acknowledgments

Work was supported by grants HL-067792 and HL-089559 (to W.P.) from the National Institutes of Health. C.E. was supported by a German Cardiac Society Research Fellowship.

Author Disclosure Statement

No competing financial interests exist for any of the authors. C.E. designed and carried out the experiments in the murine model. C.E. and C.A. developed the mathematical model and adopted it to the measured data. All authors analyzed and discussed the results and contributed to writing the manuscript.

References

Day

. Metabolic syndrome, or what you will: definitions and epidemiology. Diabetes Vasc Dis Res, 2007; 4:32–38.

Eckel

, Grundy

, Zimmet

. The metabolic syndrome. Lancet, 2005; 365:1415–1428.

Fonseca

. Clinical Diabetes—Translating Research into Practice. Philadelphia: Saunders Elsevier, 2006.

Stone

, Bilek

, Rosenbaum

. Recent National Cholesterol Education Program Adult Treatment Panel III update: adjustments and options. Am J Cardiol, 2005; 96:53–59.

American Diabetes Association:. Consensus Development Conference on Insulin Resistance. Diabetes Care, 1998; 21:310–314.

World Health Organization. Definition and Diagnosis of Diabetes Mellitus and Intermediate HyperglycemiaReport of a WHO/IDF ConsultationGeneva: World Health Organization, 2006.

Renard

, Kramer

, Johansson

, Lamharzi

, Tannock

, von Herrath

, Chait

, Bornfeldt

. Diabetes and diabetes-associated lipid abnormalities have distinct effects on initiation and progression of atherosclerotic lesions. J Clin Invest, 2004; 114:659–668.

Ayala

, Bracy

, McGuinness

, Wasserman

. Considerations in the design of hyperinsulinemic-euglycemic clamps in the conscious mouse. Diabetes, 2006; 55:390–397.

Heikkinen

, Argmann

, Champy

, Auwerx

. Evaluation of glucose homeostasis. Curr Protoc Mol Biol, 2007; Chapter 29Unit 29B3.

10.

Pacini

, Bergman

. MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsively from the frequently sampled intravenous glucose tolerance test. Comput Methods Programs Biomed, 1986; 23:113–122.

11.

Yamashita

, Freigang

, Eberle

, Pattison

, Gupta

, Napoli

, Palinski

. Maternal immunization programs postnatal immune responses and reduces atherosclerosis in offspring. Circ Res, 2006; 99:e51–e64.

12.

Eberle

, Merki

, Yamashita

, Johnson

, Armando

, Quehenberger

, Napoli

, Palinski

. Maternal immunization affects in utero programming of insulin resistance and type 2 diabetes. PLoS One, 2012; 7:e45361.

13.

Hevener

, Olefsky

, Reichart

, Nguyen

, Bandyopadyhay

, Leung

, Watt

, Benner

, Febbraio

, Nguyen

, Folian

, Subramaniam

, Gonzalez

, Glass

, Ricote

. Macrophage PPARγ is required for normal skeletal muscle and hepatic insulin sensitivity and full antidiabetic effects of thiazolidinediones. J Clin Invest, 2007; 117:1658–1669.

14.

Ishibashi

, Goldstein

, Brown

, Herz

, Burns

. Massive xanthomatosis and atherosclerosis in cholesterol-fed low density lipoprotein receptor-negative mice. J Clin Invest, 1994; 93:1885–1893.

15.

, Binder

, Gutierrez

, Brown

, Plotkin

, Pattison

, Valledor

, Davis

, Willson

, Witztum

, Palinski

, Glass

. Differential inhibition of macrophage foam cell formation and atherosclerosis in LDL receptor-deficient mice by PPARα, β/δ and γ J Clin Invest, 2004; 114:1564–1576.

16.

Reaven

, Palinski

. The use of animal models to study diabetes and atherosclerosis and potential anti-atherosclerotic therapies. Johnstone

, Veves

. Contemporary Cardiology: Diabetes and Cardiovascular Disease, 2 ^nd. Totowa, NJ: Humana Press, 2005; 259–280.

17.

Tangirala

, Rubin

, Palinski

. Quantitation of atherosclerosis in murine models: correlation between lesions in the aortic origin and in the entire aorta, and differences in the extent of lesions between sexes in LDL receptor-deficient and apolipoprotein E-deficient mice. J Lipid Res, 1995; 36:2320–2328.

18.

Ishibashi

, Brown

, Goldstein

, Gerard

, Hammer

, Herz

. Hypercholesterolemia in low density lipoprotein receptor knockout mice and its reversal by adenovirus-mediated gene delivery. J Clin Invest, 1993; 92:883–893.

19.

Kennedy

, Ellacott

, King

, Hasty

. Mouse models of the metabolic syndrome. Dis Model Mech, 2010; 3:156–166.

20.

Merat

, Casanada

, Sutphin

, Palinski

, Reaven

. Western-type diets induce insulin resistance and hyperinsulinemia in LDL receptor-deficient mice but do not increase aortic atherosclerosis compared with normoinsulinemic mice in which similar plasma cholesterol levels are achieved by a fructose-rich diet. Arterioscler Thromb Vasc Biol, 1999; 19:1223–1230.

21.

Eberle

, Palinski

, Ament

. Identifikation und Modellierung komplexer Stoffwechselvorgänge am Beispiel der Insulinresistenz ermöglicht frühzeitige Diagnose, Verlaufsbeobachtung und Vorhersageanalyse des Diabetes mellitus Typ II [abstract P-36] Diabetol Stoffwech, 2012; 7:S36.

22.

Bergman

, Ider

, Bowden

, Cobelli

. Quantitative estimation of insulin sensitivity. Am J Physiol, 1979; 236:E667–E677.

23.

Bergman

, Phillips

, Cobelli

. Physiologic evaluation of factors controlling glucose tolerance in man: measurement of insulin sensitivity and beta-cell glucose sensitivity from the response to intravenous glucose. J Clin Invest, 1981; 68:1456–1467.

24.

Eberle

, Ament

. The Unscented Kalman Filter estimates the plasma insulin from glucose measurement. Biosystems, 2011; 103:67–72.

25.

Eberle

, Ament

. Identifiability and online estimation of diagnostic parameters with in the glucose insulin homeostasis. Biosystems, 2012; 107:135–141.

26.

Eberle

, Ament

. Real-time state estimation and long-term model adaptation: a two-sided approach toward personalized diagnosis of glucose and insulin levels. J Diabetes Sci Technol, 2012; 6:1148–1158.

27.

Cobelli

, Dalla Man

, Sparacino

, Magni

, De Nicolao

, Kovatchev

. Diabetes: models, signals, and control. IEEE Rev Biomed Eng, 2009; 2:54–96.

28.

Grodsky

. A threshold distribution hypothesis for packet storage of insulin and its mathematical modeling. J Clin Invest, 1972; 51:2047–2059.

29.

Nise

. Control Systems Engineering, 7 ^th. Hoboken, NJ: Wiley Text Books, 2013.

30.

Nelder

, Mead

. A simplex method for function minimization. Comput J, 1965; 7:308–312.

31.

Matsuda

. Measuring and estimating insulin resistance in clinical and research settings. Nutr Metab Cardiovasc Dis, 2010; 20:79–86.

32.

Ferrannini

, Mari

. How to measure insulin sensitivity. J Hypertens, 1998; 16:895–906.

33.

Kildegaard

, Randløv

, Poulsen

, Hejlesen

. The impact of non-model-related variability on blood glucose prediction. Diabetes Technol Ther, 2007; 9:363–371.

34.

Roy

, Parker

. Dynamic modeling of free fatty acid, glucose, and insulin: an extended "minimal model." Diabetes Technol Ther 2006;8:617–626. Erratum in: Diabetes Technol Ther, 2007; 9:133.

35.

Tolic

, Mosekilde

, Sturis

. Modeling the insulin-glucose feedback system: the significance of pulsatile insulin secretion. J Theor Biol, 2000; 207:361–375.

36.

Makroglou

, Li

, Kuang

. Mathematical models and software tools for the glucose-insulin regulatory system and diabetes: an overview. Appl Numeric Math, 2006; 56:559–573.

37.

Mari

. Mathematical modeling in glucose metabolism and insulin secretion. Curr Opin Clin Nutr Metab Care, 2002; 5:495–501.

38.

Stumvoll

, Mitrakou

, Pimenta

, Jenssen

, Yki-Järvinen

, Van Haeften

, Renn

, Gerich

. Use of the oral glucose tolerance test to assess insulin release and insulin sensitivity. Diabetes Care, 2000; 23:295–301.

39.

Pillonetto

, Caumo

, Sparacino

, Cobelli

. A new dynamic index of insulin sensitivity. IEEE Trans Biomed Eng, 2006; 53:369–379.