## Abstract

Insulin administration during insulin-modified intravenous glucose tolerance test (IM-IVGTT) can induce transient hypoglycemia in healthy insulin-sensitive subjects. This triggers counterregulatory reflex (CRR) responses, which influence the kinetics of glucose and nonesterified fatty acids (NEFA), and undermines the accuracy of mathematical modeling methods that do not explicitly account for CRR. The aim of this study is to evaluate mathematical models of glucose and NEFA kinetics against experimental data in the presence or absence of CRR. Thirteen healthy nondiabetic subjects underwent a standard IM-IVGTT and a modified test (GC-IM-IVGTT) with a variable glucose infusion preventing hypoglycemia. While model predictions fit very well with glucose and NEFA data from GC-IM-IVGTT, they lagged behind observations from IM-IVGTT during recovery from hypoglycemia, independently of insulinemia, which did not differ significantly between protocols. A modification to the glucose minimal model, using the glucose concentration below a threshold as a signal for CRR, improves model predictions for both glucose and NEFA. The associated increase in endogenous glucose production correlates, among various CRR hormones, mainly with the dynamics of glucagon concentration. The modified minimal models introduce new parameters that quantify strength and duration of CRR following hypoglycemia. Although CRR represents an unwanted side-effect in IM-IVGTT occurring only in insulin-sensitive subjects, this study provides new insights leading to improved procedures for estimating insulin sensitivity from IM-IVGTT, which may also allow for assessing the individual capacity of recovery from hypoglycemic events in patients treated with insulin or insulin-releasing drugs.

- insulin sensitivity
- hypoglycemia
- mathematical modeling
- nonesterified fatty acids

type 2 diabetes results from a disruption of the closed loop regulation of insulin secretion and action occurring on the basis of alterations of insulin signaling in peripheral tissues and insulin release by pancreatic β-cells. This leads to progressive derangement of lipid and carbohydrate metabolism, starting with impaired glucose tolerance and fasting glucose concentrations. Chronically increased circulating glucose and lipids are major causes of diabetes complications by inducing mechanisms summarized as glucolipotoxicity (8).

Ideally, individual diagnosis of glucose intolerance and type 2 diabetes is based on glycosylated hemoglobin and on the oral glucose tolerance test (OGTT) with empirical rules applied to timed plasma glucose measurements (1). For clinical studies, glucose and insulin measurements at baseline or during OGTT also allow assessing insulin sensitivity of glucose disappearance using simple equations (10, 11, 12). Aside from the hyperinsulinemic-euglycemic clamp performed under steady-state conditions, which is considered the gold standard for assessing insulin sensitivity, other sophisticated techniques involving fitting of mathematical models to dynamic glucose and insulin concentration profiles have been developed to obtain estimates of insulin sensitivity and other metabolically meaningful parameters. In this regard, the minimal model of glucose disappearance (MINMOD; Ref. 2) is the most widely used mathematical tool to evaluate insulin sensitivity under the dynamic conditions of the insulin-modified intravenous glucose tolerance test (IM-IVGTT). The test consists in a glucose bolus, normalized with respect to body weight, of 0.3 g/kg given at time *t* = 0, an additional intravenous insulin injection of 0.03 IU/kg at *t* = 20 min, and frequent blood sampling for glucose and insulin concentration measurements. The administration of exogenous insulin aims at extending the range of glucose dynamics in patients with impaired insulin sensitivity and/or secretion. In normal insulin-sensitive subjects, however, administration of the standardized insulin dose usually provokes hypoglycemia that triggers metabolic counterregulatory reflex (CRR) responses. In such circumstances, the MINMOD analysis provides insulin sensitivity values that are underestimated if compared with insulin sensitivity measured with the euglycemic-hyperinsulinemic glucose clamp (6). This estimation bias can be corrected by administering also exogenous glucose during the test to avoid plasma glucose from falling below basal levels. This modified test will be indicated with GC-IM-IVGTT, where GC stands for glucose clamp. The additional glucose input can be easily taken into account in the MINMOD analysis (6).

With regard to lipotoxicity, excess lipolysis and ectopic fat deposition in nonadipose tissues have led to the hypothesis that altered kinetics of nonesterified fatty acids (NEFA) may play a role in the pathogenesis of type 2 diabetes (19, 20). However, despite increasing knowledge on the cellular metabolic fate and signaling pathways of carbohydrates and lipids, the interplay between glucose and NEFA regulation in physiological and pathological conditions in obesity and type 2 diabetes is not fully understood. In the past decade, several mathematical models have been proposed for quantitatively assessing NEFA kinetics (4, 14, 18, 22) but are not widely used probably due to limitations. Thus a practical method to describe NEFA kinetics during rapidly changing metabolic conditions is still needed.

This study presents a secondary modeling analysis of previously published data (6). The data set consists of repeated IM-IVGTT experiments performed in normal subjects without and with additional glucose infusion, GC-IM-IVGTT, to avoid counterregulatory response to hypoglycemia. The relevance of these data for the present study resides in providing replicated observations of insulin, glucose, and NEFA profiles in the same individuals under different conditions. The aims of the present study are reassessing the MINMOD approach for quantifying glucose disappearance during IM-IVGTT and testing the applicability to IM-IVGTT of a model of NEFA dynamics previously proposed for OGTT (22). The goals are to underline weaknesses and strengths of a combined glucose-NEFA kinetics model and to propose suitable modifications to the model to arrive at a joint model representation of glucose and NEFA kinetics during IM-IVGTT.

## MATERIALS AND METHODS

### Subjects and Protocol

Thirteen nondiabetic volunteers [7 male and 6 female, age: 26 ± 1 yr, body mass index (BMI): 22.1 ± 0.7 kg/m^{2}] were studied in random order during standard IM-IVGTT: 0.3 g/kg glucose at time 0 and 0.03 IU/kg insulin at 20 min and during a modified test (GC-IM-IVGTT) with additional glucose infusion adjusted manually to prevent plasma glucose concentration from falling below 100 mg/dl. Insulin, glucose, and NEFA plasma concentrations were measured at frequent intervals, from 15 min before the beginning of the test and during the following 3 h. Plasma concentrations of C-peptide, glucagon, cortisol, growth hormone, epinephrine, and norepinephrine were also measured at timed intervals. The protocol was approved by the local institutional review board and written informed consent was obtained from all subjects (for more details see Ref. 6).

The design of this protocol was motivated by observations made in normal insulin-sensitive subjects during standard IM-IVGTT in which glucose concentrations dropped markedly below baseline levels following the insulin injection. This stimulated CRR that accelerated recovery of glucose and stimulated NEFA production provoking sometimes an overshoot of NEFA levels above baseline.

### Mathematical Models

Previously proposed models of glucose (2) and NEFA (22) kinetics were initially fitted to experimental data after minor changes made to the original model equations. In particular, the glucose minimal model was modified to account for exogenous glucose infusion during the GC-IM-IVGTT, and the NEFA kinetics model was reparameterized to provide a clearer distinction between insulin-dependent suppression and insulin-independent NEFA production.

Following this initial validation phase, more substantial changes were made to the model equations to fit the accelerated recovery of glucose concentrations from hypoglycemia, and the observed NEFA overshoot attributed to CRR during IM-IVGTT. The suggested model changes are such that, under the hypothesis of identical insulin concentration profiles, the model of glucose kinetics behaves identically during the initial phase of both IM-IVGTT or GC-IM-IVGTT, up to the time when either the CRR or the exogenous glucose infusion begin.

#### Models for preliminary validation.

The initially tested models of glucose and NEFA kinetics are closely related to the original model equations, as follows.

##### MODEL OF GLUCOSE KINETICS.

The minimal model of glucose disappearance, adapted from Ref. 2, is described by the differential equations
(1)
(2)
where G(*t*) (mg/dl) and I(*t*) (mU/l), are glucose and insulin concentrations, respectively; *X*(*t*) (min^{−1}) represents insulin action in a remote compartment relative to deviations from measured basal insulin I_{b}, and therefore with initial condition *X*(0) = 0; the initial glucose concentration after the bolus is G(0^{+}) = G_{0} + D_{G}/V_{G}, where G_{0} is the measured basal, preinjection glucose concentration; in contrast G_{b} is the estimated basal reference glucose concentration; D_{G} is the glucose dose of 300 mg/kg and V_{G} (dl/kg) is the apparent glucose distribution volume normalized with respect to body weight; *u*_{G}(*t*) (mg·kg^{−1}·min^{−1}) is the exogenous glucose infusion, which is either zero during IM-IVGTT, or imposed manually during GC-IM-IVGTT; *S*_{G} (min^{−1}) represents glucose effectiveness (fractional clearance) at basal insulin; *S*_{I} (l·mU^{−1}·min^{−1}) is the insulin sensitivity index that quantifies the change in fractional glucose clearance with respect to unit changes in insulin concentration at steady state; and *P*_{2} (min^{−1}) determines the time lag between deviations of insulin from basal and changes in insulin action.

Parameter estimates of {*S*_{G}, *S*_{I}, V_{G}, G_{b}, *P*_{2}} were obtained by fitting the model predictions G(*t*), calculated at the blood sampling times, to glucose concentration measurements. Basal concentrations G_{0} and I_{b} were measured; the forcing input I(*t*) was obtained by linear interpolation of insulin concentration measurements; and the known glucose input *u*_{G}(*t*) was represented as piecewise constant infusions.

##### MODEL OF NEFA KINETICS.

The kinetics of NEFA concentration in plasma is described by the following differential equations
(3)
(4)
where N(*t*) (μmol/l) represents plasma NEFA concentration with fractional turnover rate *K*_{N} (min^{−1}); *Y*(*t*) (dimensionless) represents the signal for NEFA suppression in a remote compartment due to insulin deviations from basal I_{b}, with time lag characterized by *P*_{2}^{N} (min^{−1}); initial conditions are *Y*(0) = 0 and N(0) = N_{Ii} + N_{Id}, where N_{Ii} and N_{Id} (μmol/l) are insulin-independent and -dependent NEFA concentrations at basal insulin, respectively, i.e., N_{Ii} is the residual NEFA concentration at maximum suppression by insulin and N_{Id} is the maximal concentration gap; NEFA suppression by insulin is modeled in *Eq. 3*, as [1 − *Y* ]^{+}, where [·]^{+} represents the positive part of the argument, i.e., [*x*]^{+} = *x* if *x* ≥ 0 and [*x*]^{+} = 0 if *x* < 0; the steepness of the suppression line is determined by *S*_{I}^{N} (l/mU; *Eq. 4*), which represents thus the sensitivity of NEFA suppression with respect to unitary increase of insulin concentration around basal. Moreover, maximal NEFA suppression is achieved for *Y*(*t*) ≥ 1 and therefore for prolonged deviations of insulin from baseline equal to or greater than (*S*_{I}^{N})^{−1}. On the other side, since insulin effects are modeled around basal, a theoretically possible decline of insulin concentrations below I_{b} determines, according to the model, an increase in NEFA concentration. This increase is however, limited by the fact that I(*t*) ≥ 0 and thus *Y*(*t*) ≥ −*S*_{I}^{N}I_{b}. The ratio between maximum increase of NEFA at zero insulin and NEFA at basal insulin is given by (1 + *S*_{I}^{N}I_{b}), if N_{Ii} = 0.

As for the glucose minimal model, parameter estimates of {*K*_{N}, *S*_{I}^{N}, N_{Ii}, N_{Id}, *P*_{2}^{N}} were obtained by fitting the model predictions N(*t*) to available NEFA concentration measurements at given sampling times, with measured basal concentration I_{b}, and with forcing input I(*t*), obtained by linear interpolation of insulin concentration measurements.

Noteworthy are the similarities between the glucose minimal model and the NEFA kinetics model, in particular, that the mathematical descriptions of insulin action dynamics, *Eqs. 2* and *4*, are structurally identical. This will allow assessing the association of insulin action on glucose disappearance with that on NEFA suppression. Moreover, the model structure *Eq. 2* will be assumed as a general model of hormonal action mechanism.

#### Combined models in the presence of CRR.

The preliminary fitting of models of glucose and NEFA kinetics to experimental data showed better predictions of GC-IM-IVGTT measurements than of IM-IVGTT. During the latter protocol, counterregulatory hormones increased (6), whereas insulin concentration profiles did not change significantly compared with the GC-IM-IVGTT. In addition, model-derived glucose and NEFA concentrations lagged behind measurements. Taken together, these findings support the hypothesis of the presence of CRR and that its effects are significant for both glucose and NEFA. Attempts to incorporate into the kinetic models the profiles of counterregulatory hormones did not give satisfactory results, possibly because of unexplainable variability of hormone measurements and sparse sampling. For this reason, the role and significance in the present study of these hormones were assessed independently as described later.

A practical and effective representation of the CRR was obtained by describing insulin-independent glucose and NEFA production proportional to the response variable *R*(*t*), calculated from the following differential equation
(5)
where *R*(*t*) ≥ 0 (mg/dl) represents the memory of the CRR, with initial condition *R*(0) = 0; with time lag determined by *K*_{R}, (min^{−1}) and value proportional to the amplitude and duration of below-basal glucose excursions, quantified by [G(*t*) − G_{b}]^{−}, where [·]^{−} represents the negative part of the argument, i.e., [*x*]^{−} = −*x* if *x* ≤ 0 and [*x*]^{−} = 0 if *x* > 0.

The physiological assumption made in *Eq. 5* is that the CRR is initiated by glucose falling below the basal glucose threshold G_{b}. If this assumption were true, it would imply the existence of a biological sensor that measures absolute glucose concentrations around a threshold rather than, for instance, rates of changes in glucose concentrations that would allow for a more rapid anticipatory CRR. Such feedforward mechanisms exist with regard to increases of glucose concentration following oral intake, i.e., gut factors.

Glucose and NEFA response to CRR following hypoglycemia, was represented as the stimulated endogenous glucose production rate: EGP(*t*) = *P*_{G}*R*(*t*) (mg·dl^{−1}·min^{−1}) and the NEFA concentration: *P*_{N}*R*(*t*) (μmol/l), respectively. The difference in measurement units between the two response terms depends on the way glucose and NEFA kinetics have been defined. In particular, *R*(*t*) is expressed in units of mg/dl of glucose and *P*_{G} as a fractional rate constant with units of min^{−1}; while *P*_{N}*R*(*t*) represents a set point in the NEFA kinetic model, with *P*_{N} being a sensitivity with units of μmol·l^{−1}·dl·mg^{−1}.

The combined models of glucose and NEFA kinetics can be summarised as
(6)
(7)
(8)
(9)
(10)
with initial conditions G(0) = G_{0} + D_{G}/V_{G}; N(0) = N_{Ii} + N_{Id}, and *X*(0) = *R*(0) = *Y* (0) = 0. Estimated model parameters are {*S*_{G}, *S*_{I}, *P*_{2}, *V*_{G}, *G*_{b}, *K*_{N}, *P*_{2}^{N}, *S*_{I}^{N}, N_{Ii}, N_{Id}, *K*_{R}, *P*_{G}, *P*_{N}} while G_{0} and I_{b} are imposed to measured values. Model parameters were estimated by fitting the combined model, simultaneously, to glucose and NEFA concentrations measurements with insulin, I(*t*), as a model input obtained by linear interpolation of concentration measurements, and the external glucose infusion given as piecewise constant input, or zero.

#### Assessment of the role of CRR hormones.

The effects of CRR hormones on glucose and NEFA kinetics were initially investigated using the standard hormone action model described by first order linear dynamics forced by deviations from basal of the generic hormone *H*(*t*), i.e., analogously to the insulin action models *Eqs. 7* and *10*. Referring to stimulated endogenous glucose production, EGP(*t*) = *P*_{G}*R*(*t*), such hormone action model becomes
(11)
where *H*_{b} is the basal reference concentration; *K*_{H} is the rate constant of action; and *S*_{H} is the sensitivity of EGP with respect to unit variations of hormone *H*(*t*) from basal. Irrespective of whether such a single-hormone representation is valid or not, a drawback is that it requires the time course of *H*(*t*), which is determined, typically, by linear interpolation of concentration measurements. This may give only a coarse approximation of the true profile *H*(*t*) if sampling is sparse. The additional problem arising with multiple CRR hormones is that the effects of each single hormone may be difficult to disentangle from the others.

It may be for these reasons that the empirical description of CRR based on *Eq. 8*, combined with the two parameters *P*_{G} and *P*_{N}, performed much better in finding a solution of the model identification problem than using multiple single-hormone effects described according to *Eq. 11*. In addition to the trivial advantage of *Eq. 8* over *Eq. 11* that it does not require frequent CRR hormone concentration measurements to reconstruct their time profile, *H*(*t*), it does not rely on implicitly expected dynamic behaviors, such as constant, basal values during the initial phase of the experiment to exclude changes not attributable to hypoglycemia. Nevertheless, without accounting explicitly for a direct effect of CRR hormones in model *Eqs. 6*-*10*, it is not possible to make any inference as to which hormone fits within the proposed modeling paradigm best or to support interpretations of the modeling results.

To overcome these limitations, we consider a single-hormone effect on EGP(*t*) according to *Eq. 11* and express *H*(*t*) as a function of model parameters and the time courses of EGP(*t*) and dEGP(*t*)/d*t*, which are regarded as known from *Eqs. 6*–*10*. The resulting inverse problem is obtained by rearranging *Eq. 11* into the linear regression model
(12)
where subscript *i* represents the *i*-th experiment, *t*_{k} ∈, {*t*_{1}, . . ., *t*_{ni}} is the sampling schedule, *H*_{i}(*t*_{k}) is the measured response variable (hormone), and EGP_{i}(*t*_{k}) and dEGP_{i}(*t*_{k})/d*t* are the regressors calculated using individual parameter estimates in *Eqs. 6*–*10*.

The relationship between linear regression coefficients {α, β, γ} and model parameters {*H*_{b}, *S*_{H}, *K*_{H}} in *Eq. 11* is evident. However, positivity constraints on the regression coefficients pose less obvious requirements on the solution such as causality, i.e., that hormone effects lag behind its variations from basal. This condition was used to evaluate the plausibility of the solution of the linear regression problem, together with common statistical criteria such as precision of parameter estimates and accuracy of model predictions.

### Statistical Methods

The statistical analysis of glucose and NEFA data by means of the kinetic models was carried out within a nonlinear mixed-effects (NLME) population parameter identification framework. The NLME approach provides point estimates of so-called fixed effects that include average population values of kinetic parameters and estimates of the effects of covariates or individual characteristics, as well as second-order statistics of so-called random effects that characterize within and between-subject variability of parameters. For this purpose, mixed-effects models are fitted simultaneously to the available data collected in all subjects. In short, the population modeling approach describes experimental data collected in the *i*-th subject as
(13)
where **y**_{i} is the vector of all measurements, **f**(**t**_{i}, θ_{i}) are the model predictions at the sampling times **t**_{i} for given parameter values θ_{i}, and **e**_{i} represents measurement noise. Individual parameter estimates, θ̂_{i}, and population summaries, e.g., mean and standard deviation, can be obtained, in principle, using the traditional two-stage approach in which individual parameter estimates θ̂_{i} are obtained first by nonlinear weighted least squares fitting, and then analyzed as a new set of data. A more dependable approach consists in formulating a hierarchical model that describes individual responses according to *Eq. 13*, as well as the effects of covariates on selected components of the parameter vector θ_{i}, typically as a linear equation as
(14)
where *k* refers to the *k*-th component of the individual parameter vector θ_{i}; θ_{0k} represents the mean population value or the value characterizing a reference group; *X*_{ik} are the explicative variables (covariates) to predict, deterministically, interindividual variability by means of the *fixed effects*, α_{k}; and δ_{ik} are the *random effects* that represent unpredictable inter-individual variations (16). Instead of θ_{ik} the actual estimated parameters become thus θ_{0k} and α_{k}. The selection of the explicative variables in *Eq. 14* is carried out through the typical iterative procedure adopted in multiple linear regression modeling, based on the significance level attached to regression coefficient α_{k} with null hypothesis of α_{k} being equal to zero, within a tradeoff strategy between accuracy of model predictions and model complexity.

It is worth noting that θ_{0k} may also refer to the reference protocol, e.g., IM-IVGTT, with *X*_{ik} becoming an indicator variable, i.e., *X*_{ik} = 0 for IM-IVGTT and *X*_{ik} = 1 for GC-IM-IVGTT, such that α_{k} quantifies between-protocol differences. Moreover, the linear relationship, *Eq. 14*, is not a too strong limitation as it can be applied together with nonlinear transformations of parameters, e.g., θ_{ik} may represent the logarithm of a kinetic model parameter, such that the additive model, *Eq. 14*, becomes a multiplicative relationship for the back-transformed parameter.

Glucose and NEFA data were used either separately for identifying the mixed-effects parameters of the single-output glucose or NEFA kinetic models, respectively, or jointly for identifying the combined glucose-NEFA kinetic model. In particular, the glucose, *Eqs. 1* and *2*, and NEFA, *Eqs. 3* and *4*, kinetic models were fitted separately to all available glucose and NEFA data, respectively, excluding from the analysis glucose concentrations measured before 8 min to allow for complete mixing and NEFA measurements after 90 min in four subjects who showed a NEFA overshoot above basal. On the contrary, the complete model, *Eqs. 6*-*10*, was fitted to all available glucose and NEFA data simultaneously by excluding only the glucose measurements before 8 min. The actual estimated parameters and the associated statistical figures were obtained for the logarithms of the original model parameters, yet the results are presented, back-transformed, for the original parameterization.

Measurement errors for glucose and NEFA concentrations were assumed heteroschedastic following a power variance function with corresponding coefficients estimated along with the main model parameters. In the estimation procedure, weights of prediction errors, i.e., (observed − fitted), were thus calculated as the inverse of a power of the prediction itself. When used together with the combined model, NEFA and glucose data had two distinct power variance function coefficients. Population parameter estimates were obtained using the NLME package of the statistical software R (16, 17). The various glucose and NEFA kinetics models were defined and simulated using the modeling tool Pansym (21).

The inverse problem of reconstructing CRR hormone profiles from sparse hormone measurements was solved, using individually estimated time courses of EGP and dEGP/d*t*, within the population modeling framework of linear mixed-effects parameter estimation using the NLME package (16, 17). Basal hormone concentrations measured before administration of the glucose bolus were excluded from the analysis to ignore the immediate effect of a glucose bolus on hormone concentrations. Population mean parameter estimates (fixed effects) of {α, β, γ} (*Eq. 12*) and random effects variances were obtained for each measured CRR hormone, but only statistically significant results are reported.

The nonparametric, two-sample two-sided Kolmogorov-Smirnov (KS) test was used to assess the effect of protocol, i.e., IM-IVGTT vs. GC-IM-IVGTT, on the cumulative distribution of measurements of various substances.

The null hypothesis of statistical tests was rejected with cut-off levels of *P* = 0.05.

## RESULTS

### Exploratory Data Analysis

Average glucose and NEFA concentration profiles during the two protocols are depicted in Fig. 1. The glucose undershoot below baseline during IM-IVGTT and the efficacy of exogenous glucose infusion during GC-IM-IVGTT are clearly seen in Fig. 1, *A* and *B*, respectively. Average levels of NEFA concentrations during IM-IVGTT were, at baseline, lower than during GC-IM-IVGTT but reached at the end of the test higher values than at baseline as well as throughout the GC-IM-IVGTT. This observation suggested expressing short-term regulation of NEFA concentrations in relative terms with respect to baseline rather than in absolute terms.

Increased average end-test NEFA concentrations during IM-IVGTT were primarily due to NEFA overshoot observed in four subjects. For the initial modeling analysis, NEFA data beyond 90 min were therefore discarded for these four subjects. Instead, the final model that accounts for CRR was fitted against all available NEFA concentrations. On the contrary, in all modeling analyses involving glucose, measurements were considered starting from the 8-min sample. Steady state glucose concentrations were slightly higher during GC-IM-IVGTT compared with IM-IVGTT (Fig. 1, *A* vs. *B*). This justifies a higher endogenous insulin secretion during GC-IM-IVGTT vs. IM-IVGTT, as demonstrated indirectly by the C-peptide levels shown in Fig. 2, *D* and *C*, respectively. Nevertheless, insulin concentrations were equally influenced by the first phase insulin secretion and more importantly by the exogenous insulin injection (Fig. 2, *A* and *B*).

The empirical distributions of pooled insulin measurements collected in all subjects during IM-IVGTT and GC-IM-IVGTT were compared with the KS test, which did not evidence significant deviations from the identity (results not shown). This indicates that insulin concentration profiles were only marginally affected by the protocol if at all. Any between-protocol difference in experimental or modeling results, except for the highly specific C-peptide measure, can therefore hardly be attributed to changes in insulin concentrations.

The KS test applied to glucose and NEFA concentrations indicated, as expected, statistically significant differences in the measurement distributions between the two protocols. In particular, glucose was significantly higher during GC-IM-IVGTT than during IM-IVGTT at concentrations <100 mg/dl, whereas NEFA concentrations were higher during IM-IVGTT in the range between 70 and 300 μmol/l. This strengthened the qualitative observation that NEFA concentrations were higher during IM-IVGTT compared with GC-IM-IVGTT and, if taken together with the finding that distributions of insulin concentrations are comparable, this supports the hypothesis of intervention of insulin-independent CRR mechanisms for NEFA production during IM-IVGTT.

### Preliminary Modeling Analysis

The minimal model of glucose disappearance, *Eqs. 1* and *2*, and the model of NEFA kinetics, *Eqs. 3* and *4*, were first fitted, individually, to glucose and NEFA concentrations, respectively. This preliminary analysis evidenced discrepancies of behavior between the two protocols, endorsing the hypothesis of a significant effect of CRR. Limitations of the original models in describing the respective substrate dynamics were in fact not confined to poor fitting quality but showed important between-protocol variations of parameter estimates that partly compensated the different underlying kinetics. Preliminary results obtained with the glucose and NEFA kinetic models are presented separately in the following paragraphs.

#### Glucose minimal model.

Population kinetic modeling of glucose dynamics gave, in terms of average quality of fit, better results with data from GC-IM-IVGTT than from IM-IVGTT. Model prediction errors (Δ = observed − fitted) obtained with model *Eqs. 1* and *2* are connected in Fig. 3, *A* and *B*, with dashed lines and compared with the ±SE intervals of glucose measurements during protocols IM-IVGTT (*A*) and GC-IM-IVGTT (*B*), respectively. According to the 68% (SE) confidence intervals of mean concentration values, Fig. 3 shows that the model fit is on average acceptable as regards glucose data from the GC-IM-IVGTT (Fig. 3*B*), while model predictions are biased for the IM-IVGTT in the central part of the experiment with overestimation (Δ < 0) ∼60 min and underestimation (Δ > 0) ∼100 min.

Numerical values of parameter estimates are summarized in Table 1. Although the minimal model performs better with data from GC-IM-IVGTT, numerical results are presented with main reference to IM-IVGTT. Given that parameter variability was assessed using log-linear models, Δ*S*_{I} in Table 1 represents the multiplicative factor of *S*_{I} for GC-IM-IVGTT, compared with IM-IVGTT, irrespective of gender. This means that insulin sensitivity estimated in a same subject with GC-IM-IVGTT is expected to be 2.6 times (95% range [2, 3.4]) the value estimated with IM-IVGTT. Moreover, the apparent glucose distribution volume V_{G} (dl/kg) was found significantly correlated with BMI. In the specific case, the nominal value of V_{G} was calculated for a reference BMI of 24 kg/m^{2}, such that V_{G} varies according to the formula: V_{G}(BMI) = V_{G}·(∇V_{G})^{BMI−24}.

#### NEFA kinetics model.

Population kinetic modeling of NEFA dynamics gave better results, in terms of average quality of fit and similar to glucose, with data from GC-IM-IVGTT than from IM-IVGTT. Model prediction errors obtained with model *Eqs. 3* and *4* are connected in Fig. 3, *C* and *D*, with dashed lines and compared with the ±SE intervals of NEFA measurements during protocols IM-IVGTT (*C*) and GC-IM-IVGTT (*D*), respectively. It can be observed that during IM-IVGTT the model overestimates (Δ < 0) the NEFA nadir ∼50 min. The corresponding fixed-effects estimates are reported in Table 2. All between-subject variability of kinetic parameters was described by random effects only, as no covariate was found to significantly affect the fixed effects. The estimated value of insulin-independent NEFA concentration, N_{Ii}, was markedly lower, but not negligible, compared with the insulin-dependent NEFA concentration, N_{Id}, the former accounting for ∼13.7% of total NEFA at baseline. The estimated NEFA turnover rate, *K*_{N}, was in physiological range with a half-life of 4.2 min. The reciprocal of the insulin sensitivity index, *S*_{I}^{N}, provides an estimate of the least deviation of insulin from basal that switches insulin-dependent NEFA production completely off. This resulted equal to 17.7 mU/l.

### Combined Glucose-NEFA Model in the Presence of CRR

The combined minimal model of glucose and NEFA kinetics that accounts for CRR following insulin administration (*Eqs. 6*–*10*) was fitted to all available glucose and NEFA data, i.e., including the four cases with NEFA overshoot that were partially excluded in the preliminary analysis. For a clearer evaluation of the improvements introduced by the combined model, Fig. 3 shows the sequences of average prediction errors obtained with the original models (dashed lines) and with the combined model (solid lines). As expected, the two modeling approaches fit similarly well both mean glucose and NEFA concentrations during GC-IM-IVGTT (Fig. 3, *B* and *D*) and provide a significant refinement of the model predictions of glucose and NEFA during standard IM-IVGTT (Fig. 3, *A* and *C*). The main adjustment in the model equations is the introduction of the variable *R*(*t*), which determines, simultaneously, the profile of endogenous glucose and NEFA production stimulated by CRR.

The fixed-effects estimates of the combined model are reported in Table 3. It can be observed that nominal values, as well as dependencies of fixed effects on covariates, have changed with respect to the preliminary results (Tables 1 and 2). In particular, basal reference glucose concentration, G_{b}, resulted smaller in the combined model and significantly different between males and females; the remote insulin action rate constant *P*_{2} essentially confirmed the value estimated with the original model for the GC-IM-IVGTT group, as no significant difference was found between the two protocols; glucose effectiveness, *S*_{G}, behaved similarly between the two modeling approaches; and insulin sensitivity of glucose dissappearance, *S*_{I}, exhibits the most important differences between the two models as regards dependency on gender and test type (Table 1). It is worth noting that a unique value of the fixed-effect *S*_{I} in Table 3 does not exclude between- and within-subject variability of this parameter, which was effectively modeled by random effects, but only that there are no significant predictors of between-subject and between-protocol variability for the data set at hand. Finally, the apparent glucose distribution volume and its dependency on BMI were similar for the two models.

With regard to parameters that characterize NEFA kinetics, insulin-dependent and -independent NEFA concentrations, N_{Id} and N_{Ii}, further diverged and the fraction of insulin-independent NEFA concentration at basal declined to 7.8%, which is a consequence of improved fit of NEFA nadir during IM-IVGTT; insulin sensitivity of NEFA suppression, *S*_{I}^{N}, remained essentially unchanged; and for the combined model, the rate constants of insulin action, *P*_{2}^{N}, was found significantly different between IM-IVGTT and GC-IM-IVGTT.

Regarding the issue as to whether the rate constants of insulin action on glucose disappearance or on NEFA suppression are the same or not, from Table 3 it follows that the 95% confidence interval of *P*_{2} overlaps that of *P*_{2}^{N} for the GC-IM-IVGTT but not for the estimate of *P*_{2}^{N} determined from IM-IVGTT. With the caveat of a possible fortuitous coincidence, this leads to the conclusion that, in the absence of CRR, control of glucose disappearance and NEFA suppression share the same direct or indirect insulin signaling pathway. On the other hand, the significant difference of *P*_{2}^{N} between the two protocols suggests that the differences in insulin action on NEFA kinetics are not completely resolved by adding the CRR model, *R*(*t*) in *Eq. 9*.

#### Solutions of inverse problems for CRR hormone profiles.

Figure 4, *A* and *B*, shows the average model-reconstructed time courses of stimulated endogenous glucose production, EGP(*t*), during the two protocols. A marginal CRR activity was estimated also during GC-IM-IVGTT. Glucagon concentrations measured during the two protocols and the corresponding inverse problem solutions are summarized in Fig. 4, *C* and *D*. Glucagon is the only CRR hormone for which all three regression coefficients of the inverse problem resulted nonnegative and highly significantly different from zero. Numerical results of fixed and random effects and figures of precision are the following: α = 58.66 ± 3.06 ng/l (fixed effect ± SE), *P* < 0.0001 (SD = 15.3; SD of random effect); β = 7.882 ± 1.84 ng/l/(mg·dl^{−1}·min^{−1}), *P* < 0.0001 (SD = 6.72); γ = 88.36 ± 25.9, ng/l/(mg·dl^{−1}·min^{−2}), *P* < 0.001 (SD = 73.7).

With regard to other measured CRR hormones, i.e., cortisol, growth hormone, norepinephrine, and epinephrine, results exhibited a high variability of concentration profiles and a poor prediction of the data (not shown).

In summary, among the various measured hormones only glucagon supports the hypothesis that EGP(*t*) may be expressed interchangeably through *Eq. 8* as a function of glucose concentrations below basal or through *Eq. 11* as a causal, delayed effect of suprabasal hormone variations, yet with the appropriate time courses summarized in Fig. 4, *C* and *D*. However, if compared with the model reconstructed profile, all hormone measurements, including glucagon, present a large between-subject variability (±SE vs. mean).

## DISCUSSION

Various mathematical models have been proposed in the past for quantitative assessment of NEFA kinetics (4, 14, 18, 22). These differ mainly in the mathematical description of the factors regulating the rate of NEFA appearance; in particular, the control mechanisms for lipolysis have been described in terms of *1*) insulin-action alone (14, 22), which is in line with current knowledge on antilipolytic effects of insulin even at basal glucose (3); *2*) glucose-action alone (4), which is theoretically barely defendable but is supported in practice by successful prediction of NEFA concentration time curves (5); or *3*) the action of both insulin and glucose (18), however, with glucose assumed to stimulate lipolysis rather than potentiate the anti-lipolytic effect of insulin observed experimentally (15).

Several competing models have been compared in Periwal et al. (14), but the statistical methods used therein for parameter estimation and model comparison led to conclusions and recommendations that are likely biased towards a too complex model of NEFA kinetics. In fact, to overcome practical identifiability problems, the Bayesian parameter estimation scheme adopted in Periwal et al. relied on “prior knowledge” about estimated model parameters, which is unavailable during early stages of model developments. Moreover, Periwal et al. represented the pharmacodynamic effect of insulin on NEFA production by means of a Hill-type equation, which exhibits, depending on the model parameter values, either a zero or an infinite slope of the effect on NEFA of insulin variations around zero or basal. This aspect does not allow the definition of a physiologically useful sensitivity index of NEFA production with regard to insulin variations around basal insulin.

This study provides a new integrated perspective on glucose and NEFA kinetics and their control by insulin. The main results concern: *1*) a modification to the minimal model of glucose disappearance to account for CRR during IM-IVGTT; *2*) the validation of a previously proposed model of NEFA kinetics; and *3*) the proposal of a new combined model of glucose and NEFA kinetics, which fixes, on one side, a prediction bias for both glucose and NEFA during IM-IVGTT, and provides, on the other, a tool for the assessment of CRR to insulin-induced hypoglycemia. This latter may find application as diagnostic tool in insulin treatment.

The preliminary, exploratory analysis of IM-IVGTT vs. GC-IM-IVGTT data has confirmed the expectation that glucose profiles are significantly different but has provided evidence of a significant difference between NEFA concentration profiles. On the contrary, the distribution of insulin concentrations was not significantly different between the two protocols, despite the fact that endogenous insulin secretion, assessed indirectly through circulating C-peptide levels, was significantly higher during the late phase of the GC-IM-IVGTT.

The modeling analysis of insulin and glucose concentrations profiles, from the two protocols, with the minimal model of glucose disappearance (2), produced on average an acceptable fit, yet with appreciable prediction bias for the IM-IVGTT data (Fig. 3). Major inconsistencies emerged, however, rather from the structure and numerical values of the fixed-effects parameter estimates (Table 1). In particular, the rate constant of insulin action, *P*_{2}, was significantly higher during IM-IVGTT, and, more importantly, the insulin sensitivity index *S*_{I} result was 2.6 times higher during GC-IM-IVGTT. In summary, differences in glucose kinetics, with equivalent insulin profiles and accounting for exogenous glucose infusion during GC-IM-IVGTT, were compensated by differences in parameter values that can be interpreted in terms of protracted (smaller *P*_{2}) and stronger (increased *S*_{I}) insulin action. These results are inconsistent because the two protocols have to be considered identical during the initial phase until glucose infusion is started during GC-IM-IVGTT.

With regard to the analysis of insulin and NEFA concentration profiles by means of the model originally proposed for OGTT data (22), best fit model predictions were acceptable for both protocols, again with a slight prediction bias during the nadir phase of the IM-IVGTT. Fixed effects were limited to the main model parameters, i.e., NEFA kinetics was not significantly affected by any considered covariate, and all between- and within-subject variability was attributed to random effects. These results may be due to large variability of NEFA responses between different experiments. Nevertheless, the model of NEFA kinetics provided physiologically relevant and plausible parameter values such as the NEFA turnover rate, *K*_{N}; the rate constant of insulin action on NEFA suppression, *P*_{2}^{N}; and the insulin-dependent (suppressible) and -independent NEFA concentrations at baseline.

From these preliminary results, it became clear that the main conceptual problems were in the description of glucose kinetics by the minimal model of glucose disappearance, rather than in the mathematical description of NEFA kinetics, probably due to a much higher reproducibility of glucose time courses compared with NEFA. Among various possible competing alternatives that have been tested, the model *Eqs. 6*–*8* represent the final refinement proposed for the minimal model of glucose disappearance to account for CRR in the presence of hypoglycemia induced by exogenous insulin administration in normal insulin-sensitive subjects. These three differential equations are able to improve the average fit of glucose concentration during IM-IVGTT as shown in Fig. 3*A*, whereas no modification to the minimal model appears necessary in the absence of CRR (Fig. 3*B*). Thus *Eqs. 6* and *7*, with *R* = 0, should remain valid for analyzing regular IVGTT without exogenous insulin administration. Nevertheless, the proposed modification not only serves as a correction factor for the estimation bias in insulin sensitivity but represents a conceptual innovation in the description of blood glucose regulation during IM-IVGTT that may become useful to assess in clinical practice the counterregulatory response to insulin-induced hypoglycemia, e.g., in patients treated with insulin or insulin-releasing drugs. This issue has been already investigated previously (7) but with a much more demanding experimental setup than the standard IM-IVGTT considered in this study.

Once the CRR-correction term for glucose kinetics had been established, the same mechanism was postulated for the NEFA kinetics, showing a significant improvement especially at the nadir but also appreciable reduction of the average prediction errors during the whole test (Fig. 3*C*), despite the fact that the fit was based on all NEFA concentrations, including those that had been excluded from the preliminary analysis. Fixed-effects estimates of the combined model (Table 3) showed, besides the effect of gender and BMI on a few parameters, a significant effect of the protocol only on the rate constant of insulin action on NEFA suppression, *P*_{2}^{N}. The observation that *P*_{2}^{N} during GC-IM-IVGTT is statistically indistinguishable from the companion *P*_{2} of the minimal model means that the rate constant of onset and termination of insulin action (in the absence of CRR) may be the same for glucose and NEFA and these may thus share macroscopically the same signaling pathway. The difference in *P*_{2}^{N} between the two protocols indicates that there might be some unmodeled mechanism that reduces the time lag of NEFA suppression by insulin, most likely due to unmodeled effects of counterregulatory hormones. As for the minimal model described previously, the model described by *Eqs. 9* and *10*, with *R* = 0, may be used in the absence of hypoglycemia, e.g., for OGTT, which has been the first application of the model.

The presented results support the validity of the combined model in the presence of CRR, as well as the validity of the NEFA model proposed in Ref. 22 in the absence of CRR. The latter model provides parameter estimates of physiologically relevant characteristics of NEFA kinetics and its control by insulin. However, the validity of the NEFA kinetic model has been questioned in a previous study (Ref. 14, on the basis of a Bayesian inference approach). The latter requires prior knowledge on parameter distributions that may strongly influence the outcome of the inference problem represented by the so-called posterior distribution of the parameters. Such prior knowledge on parameter distributions is in general missing during the early stages of model development, as the aim is to select the best model among alternatives that have not been tested before. Fixing arbitrarily lower and upper bounds on prior parameter distributions, as done in Ref. 14, can lead to biased estimates because the posterior parameter distribution will necessarily lie within the prior distribution interval. To impose ordinary constraints on parameters, nonlinear transformations are useful, e.g., positivity can be imposed with log transformation, as done in this study. As an alternative to Bayesian estimation with interval bounds on prior distribution, a more dependable statistical approach would have been the population modeling method based on hierarchical Bayesian models in which the between-subject distribution of the model parameters is inferred from the data along with the posterior distribution of the parameters of each single individual (9, 13). The alternative approach adopted in this study is the nonlinear mixed-effects modeling approach, which does not rely on prior information on parameter distributions (16).

Irrespective of the statistical approach used to identify a particular model and to compare it to competing alternatives, mathematical models should not be evaluated only on the basis of goodness of fit and number of model parameters but also in terms of physiological plausibility of underlying hypotheses and functionality of the produced results in characterizing different metabolic conditions in different cohorts or in single subjects. In this regard, the counterregulatory hormones of insulin had been measured with the intent to assess their effect on glucose and NEFA kinetics. However, they could not be included easily into the model probably because of sparse sampling, measurement noise, or large between-subject variability. Nevertheless, they provided a means for validating the postulated effect of CRR to hypoglycemia described in terms of lagged action of below-basal glucose concentrations. In particular, by calculating the model predictions of stimulated EGP (and its time derivative) and by assuming a typical hormone action mechanism with first order dynamics and linearized effect around basal, it was possible to formulate an inverse problem that was solved for each hormone by fitting linear mixed-effects models to concentration measurements, as if it were the single cause of CRR. The results indicate that glucagon fits within this modeling framework well, while nonconclusive results are obtained for the other CRR hormones due to lack of precision and accuracy in model predictions (not shown).

Another important evaluation criterion is model complexity, which influences the possibility to identify the model parameters. The minimalistic modeling approach adopted in this and other studies is therefore justified by the intent to individualize model parameters through estimation from data experiments. For this purpose, population modeling approaches, such as NLME and Bayesian methods, overcome many problems of practical identification by fitting to experimental data broader models that comprise the kinetic model as well as some explicit or implicit assumptions about statistical distributions of parameters. In this context, a “minimal” model may be not as minimal as required for parameter estimation from individual experiments data, and its validity may be restricted to particular types of experiments. Further studies are therefore required to assess the range of suitability of the presented minimal models of glucose and NEFA kinetics under various conditions of clinical experiments.

### Perspectives and Significance

This study addresses an underemphasized or ignored problem related to the standardized IM-IVGTT protocol combined with the MINMOD analysis that leads to underestimates of insulin sensitivity (*S*_{I}) in normal and hyperresponsive insulin-sensitive subjects. While classification or grading of *S*_{I}, e.g., as normal, increased, or reduced, remains unaffected if cutoff levels are chosen as quantiles, underestimation of normal or increased *S*_{I} compresses the range of higher values and compromises the clinical utility that *S*_{I} estimates could have in insulin dose prescription for avoiding hypoglycemic episodes. To overcome this limitation, a correction term in the MINMOD equations is proposed to account for CRR to hypoglycemia, which allows unbiased estimation of *S*_{I} if compared with estimates obtained from the more demanding GC-IM-IVGTT protocol. Although reduction or individualized adjustment of the insulin dose administered during IM-IVGTT could eliminate the problems related to CRR all together, the standard IM-IVGTT with modified minimal model equations could become a tool for estimating both *S*_{I} as well as the patient's capacity to recover from transient hypoglycemia. Within a broader view of short-term regulation of fuel substrates, NEFA kinetics during IM-IVGTT has been described with a simple mathematical model that has proven successful in describing accurately NEFA concentration profiles during IM-IVGTT and GC-IM-IVGTT, showing however, with the former protocol similar problems in the presence of CRR as MINMOD. A combined glucose and NEFA kinetic model has then been proposed that relies on the same CRR mechanism that acts seemingly in an insulin-independent way on both glucose and NEFA. The combined model may therefore become useful for simultaneous analysis of glucose and NEFA profiles during IM-IVGTT for investigating short-term regulation of the two major fuel substrates.

## GRANTS

The work of M. Roden is supported in part by a grant from the Austrian Science Foundation (FWF, P15656), the European Foundation for the Study of Diabetes (ESFD), and the Schmutzler-Stiftung. The work of K. Thomaseth and G. Pacini is supported in part by Regione Veneto (Biotech DGR 2702/10-09-04).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: K.T., A.B., G.P., and M.R. conception and design of research; K.T. and A.P. analyzed data; K.T., A.B., A.P., G.P., and M.R. interpreted results of experiments; K.T. prepared figures; K.T. and A.P. drafted manuscript; K.T., A.B., G.P., and M.R. edited and revised manuscript; K.T., A.B., G.P., and M.R. approved final version of manuscript; A.B. performed experiments.

## ACKNOWLEDGMENTS

We thank Lorenzo Finesso for helpful discussion on system theoretical implications of numerical solutions of the inverse problem.

Present address of K. Thomaseth: Institute of Electronics, Computer and Telecommunication Engineering (IEIIT), National Research Council (CNR), Padua, Italy.

## Glossary

- D
_{G} - Standard IVGTT glucose dose (300 mg/kg)
- EGP(
*t*) - Counterregulatory endogenous glucose production (mg·dl
^{−1}·min^{−1}) - G(
*t*) - Instantaneous glucose concentration (mg/dl)
- G
_{b} - Steady-state glucose concentration (mg/dl)
- G
_{0} - Initial glucose concentration (mg/dl)
- I(
*t*) - Instantaneous insulin concentration (mU/l)
- I
_{b} - Basal insulin concentration (mU/l)
*K*_{N}- Fractional NEFA turnover rate (min
^{−1}) *K*_{R}- Rate constant of CRR signal (min
^{−1}) - N(
*t*) - Instantaneous NEFA concentration (μmol/l)
- N
_{Id} - Insulin-dependent component of basal NEFA (μmol/l)
- N
_{Ii} - Insulin-independent component of basal NEFA (μmol/l)
*P*_{2}- Turnover rate constant of insulin action
*X*(*t*) (min^{−1}) *P*_{2}^{N}- Turnover rate constant of insulin action
*Y*(*t*) (min^{−1}) *P*_{G}- Glucose recovery rate from hypoglycemia (min
^{−1}) *P*_{N}- Sensitivity of NEFA increase following hypoglycemia (μmol·l
^{−1}·dl·mg^{−1}) *R*(*t*)- CRR signal associated with hypoglycemia (mg/dl)
*S*_{G}- Glucose effectiveness (fractional clearance) at basal insulin (min
^{−1}) *S*_{I}- Insulin sensitivity of glucose disappearance (l·mU
^{−1}·min^{−1}) *S*_{I}^{N}- Insulin sensitivity of inhibition of NEFA production (l/mU)
- V
_{G} - Apparent glucose distribution volume (dl/kg)
*X*(*t*)- Insulin action on glucose disposal (min
^{−1}) *Y*(*t*)- Insulin action on suppression of NEFA production (1)

- Copyright © 2014 the American Physiological Society