Modeling Insulin-Glucose Dynamics During Insulin Induced Hypoglycemia. Evaluation of Glucose Counterregulation

A dynamical network model of insulin-glucose interactions in subjects with Type I Diabetes was developed and applied to data sets for 40 subjects. Each data set contained the amount of dextrose + insulin infused and blood glucose (BG) determinations, sampled every 5 minutes during a one-hour standardized euglycemic hyperinsulinemic clamp and a subsequent one-hour BG reduction to moderate hypoglycemic levels. The model approximated the temporal pattern of BG and on that basis predicted the counterregulatory response of each subject. The nonlinear fits explained more than 95% of the variance of subjects' BG fluctuations, with a median coefficient of determination 97.7%. For all subjects the model-predicted counterregulatory responses correlated with measured plasma epinephrine concentrations. The observed nadirs of BG during the tests correlated negatively with the model-predicted insulin utilization coefficient (r = -0.51, p < 0.001) and counterregulation rates (r = -0.63, p < 0.001). Subjects with a history of multiple severe hypoglycemic episodes demonstrated slower onset of counterregulation compared to subjects with no such history (p < 0.03).

A dynamical network model of insulin-glucose interactions in subjects with Type I Diabetes was developed and applied to data sets for 40 subjects. Each data set contained the amount of dextrose + insulin infused and blood glucose (BG) determinations, sampled every 5 minutes during a one-hour standardized euglycemic hyperinsulinemic clamp and a subsequent one-hour BG reduction to moderate hypoglycemic levels. The model approximated the temporal pattern of BG and on that basis predicted the counterregulatory response of each subject. The nonlinear fits explained more than 95% of the variance of subjects' BG fluctuations, with a median coefficient of determination 97.7%. For all subjects the model-predicted counterregulatory responses correlated with measured plasma epinephrine concentrations. The observed nadirs of BG during the tests correlated negatively with the model-predicted insulin utilization coefficient (r = -0.51, p < 0.001) and counterregulation rates (r = -0.63, p < 0.001). Subjects with a history of multiple severe hypoglycemic episodes demonstrated slower onset of counterregulation compared to subjects with no such history (p < 0.03).
Kqwvrds: Dynamic models. hypoglycemia, counterregulation, hyperinsulinemic clamp Insulin Dependent Diabetes Mellitus (IDDM) occurs when the pancreas produces insufficient insulin to prevent hyperglycemia, necessitating administration of exogenous insulin by injection. Excessive insulin, relative to metabolic needs, leads to low blood glucose (BG) or hypoglycemia, BG <3.9 mM as defined by the DCCT Study Group (1993). In most individuals with IDDM, low BG triggers the release of counterregulatory hormones. This in turn prompts the release of stored glucose into the bloodstream to restore euglycemia. However, in IDDM subjects, the ability to counterregulate is frequently impaired by factors such as long-standing diabetes, autonomic neuropathy, and intensive therapy (Amiel et nl., 1987(Amiel et nl., , 1988Cryer et nl., 1985). Insufficient or absent counterregulatory responses allow BG to "Corresponding Author: Uni\wsi~y of Virginia Health Sciences Center. Bux 137, Charlottes\illc. VA 22908. USA: Tel: (804)-921-8656: Fax: (801)-924-8652: E-mail: bpk2u@Virginia.edu fall further until stupor, unconsciousness or seizure occurs. This condition, referred to as severe hypoglycemia (SH), is responsible for four percent of the deaths among patients with IDDM (DCCT Study Group, 1991). The risk for SH in IDDM has been attributed not only to relative insulin excess, but also to impaired glucose counterregulation (Cryer et al., 1985;Gerich, 1988;White et al., 1983).
Glucose metabolism has been studicd with isotopic tracer methods in animal and human studies (Brier et al., 1977;Hetenyi, Ninomiya and Wrenshall, 1966) and described by several mathematical models that included network modeling of glucose metabolism in normal ideal man (Guyton et al., 1978) and diabetic dogs (Yamasaki, Tiran and Albisser, 19841, and multicompartment models (Insel et al., 1975;Steele, Rostami and Altszuler, 1974). The minimal model, suggested as an alternative of hyperinsulinemic euglycemic clamping for measuring insulin sensitivity in vivo (Bergman et al., 1979), received considerable attention, support and critiques (Bergman et al., 1987;Cobelli et al., 1986;Cobelli, Brier and Ferrannini, 1990;Mari, 1997;Quon et al., 1994). While some of these models focused specifically upon the role of the liver (Carson and Cramp, 1976), investigation of counterregulation through insulin infusion has typically been pursued in two ways: 1) during the induction of hypoglycemia various hormones are sampled to determine whether their levels increase as BG falls, 2) whether or not BG spontaneously rises or platcaus despite the continual infusion of regular insulin (Bolli et al., 1984). These basic approaches to describing glucose counterregulation included quantifying plasma hormonal concentrations, but they did not yield a precise mathematical model of glucose counterregulation.
In this manuscript we present a mathematical model of insulin-glucose dynamics during a euglycemic hyperinsulinemic clamp and subsequent reduction of BG to hypoglycemic levels. Using this model, we approximated the temporal pattern of each subject's BG fluctuations and evaluated parameters of insulin and glucose sensitivity. On that basis we computed dynamic estimates for the onset and rate of these subject's counterregulatory responses. These results were verified by correlating for each subject the model-estimated counterregulatory dynamics with subsequently analyzed plasma epinephrine concentrations and applied to the study of relationships and group effects pertinent to SH.

Subjects
Forty subjects were recruited through newsletters, notices posted in diabetes clinics and direct physician referral. All subjects had to have diabetes for at least two ycars and have taken insulin since the time of diagnosis. There were 16 males and 24 females, with mean age 35.5 yr (SEM = 1.3), mean duration of disease 16.9 yr (SEM = 1.51, mean insulin units/kilogram per day 0.59 (SEM = 0.031, and mean glycosylated hemoglobin 8.6% (SEM = 0.3). The non-diabetic range for the glycosylated hemoglobin assay was 4.4 to 6.9%. Twenty-five subjects reported a history of multiple SH episodes in the past year while 15 subjects had no such history.

Procedure
All subjects attended orientation meetings and signed consent forms. To ensure that wbjccts' BGs were not in a low rangc (e3.9 mM) for 72 hrs prior to the study, their insulin dose was reduced by 10% and long acting insulin was discontinued 36 hrs prior to the study. Subjects were instructed to eat prophylactically 10 g of glucose whenever BG was (5.6 mM and were required to test their BG five times a day (1 hr before each meal, at bedtime and 4 hrs into their sleep). If low BG occurred, the study was rescheduled. Subjects were admitted to the University of Virginia General Clinical Research Center. Upon admission, subjects were given a physical exam, including an assessment for autonomic neuropathy. BG was maintained overnight between 5.6-8.3 mM with intravenous regular human insulin as per a previously published insulin infusion protocol (Bolli et al., 1984). Subjects were given dinner and a bedtime snack the evening before the study, but remained fasting on the morning of the study. No caffeinated beverages were consumed after hospital admission.
On the morning of the study, IV's were placed in the nondominant forearm. Insulin was continuously infused at a constant rate of 1.0 mU/kg/min and a 20% dextrose solution was infused at a variable rate to maintain BG at 6 mM. Figure 1 presents the design of the study.
During Phase 1 (euglycemia) BG was maintained between 5.6 and 8.3 mM. During Phase 2 (BG reduction), BG was steadily lowered to a target level of 2.2 mM by varying the dextrose infusion. Adjustments in dextrose infusion were made every 5 min. The insulin infusion was discontinued during Phase 3 (recovery). The protocol was discontinued if manifestations of SH occurred (e.g. severe lethargy, disorientation, confusion or inappropriate behavior). Arterialized blood (achieved by warming the hand in a heated glove to 50 "C) was sampled for glucose concentration every 5 rnin and for plasma epinephrine concentration every 10 min.

Mathematical Model of Insulin-Glucose Dynamics
A dynamical network model was developed to describe BG dynamics as a function of two principal temporal variables, insulin and dextrose infusion. Since the subjects did not eat and their physical activity was negligible during the study, neither of these parameters was assumed to influence BG. In addition to the two principal variables, a counterregulatory response was anticipated at lower BG levels. The network of functional interactions is shown in Figure 2.
The system was defined in terms of three timedependent state variables: 1) BG level, 2) insulin infusion, and 3) dextrose infusion. A network of processes provides the functional regulatory interactions responsible for BG control: BG level was positively affected by dextrose infusion (process D), as well as by the potential for replenishment when BG is low from available liver stores (process CR). The negative effect on BG level by insulin injections is denoted by process I. A regulatory loop (Reg) between liver stores and BG is implemented by way

Insulin Infusion
Liver Stores FIGURE 2 Network of functional interactions during hyperinsulinemic clamp: RG level\ is positively affected by the dextrose inl'usion (D) and by the pokntial for replenishment when BG is low (CR), and ncgativcly affected by the insulin infusion (I). A regulatory loop (Reg) between liver stores and BG inhibits the process CR at elcvatcd BG Icvcls.

Dextrose Infusion
of inhibition of process CR by elevated BG levels (i.e., a release from inhibition of process CR below some threshold BG level, thus providing for counterregulatory recovery from low BG levels by recruitment of available liver stores). Time rates of change of BG (i.e., d[BG]/ d r) were described by a nonlinear ordinary differential equation. The counterregulatory response was modeled as a release of glucose from a multi-compartment storage pool. It was assumed that during the first phase no counterregulatory response occurred. This permitted estimation of each individual's insulin-glucose dynamics parameters. During the second phase, an adaptive stepwise procedure was used to determine the onset and the rate of counterregulation.

Phase 1: Maintained euglycemia
It was assumed that: 1) the dextrose infusion influenced BG positively through an unknown dextrose conversion parameter a , and 2) BG decay rate was inversely proportional to the BG level, through an unknown insulin utilization parameter b. This led to the following nonlinear differential equation for the BG time rate of change: where D ( t ) is the variable dextrose infusion rate (mglkglmin) and I is the constant insulin infusion.
where E is a small constant, was used instead of standard l / B G ( t ) for a better computational stability.

Phase 2: BG reduction
During the second phase of the test counterregulation was anticipated and equation [A] was expanded by an additional term: We allowed the counterregulation term C R ( t ) to be a uni-or bi-modal function, corresponding to one-or two-compartment modeling. We would not impose a specific analytical form on the counterregulation function. In general it needs to be a positive pulsatile function that, depending on the subject's data, has one, two (or possibly more) additive components: For this particular application in order to be able to approximate our data, we suggest each counterregulation component to be defined by: This way the function CRl (t), would be zero when t < T I and would increase at time T I . Thus, the parameter T I would interpreted as time of onset of counterregulation, while the product C l .rf would be the counterregulation slope at onset (the derivative at time t = T I ) . The same would be valid for the second component CR2(t). The analytical form of CR(t) was carefully selected to allow for this physiological interpretation of its parameters, however it is not restricted and other functions that meet certain mathematical requirements would be suitable descriptors of counterregulatory responses.

Parameter Estimation
The input data for the model were each subject's dextrose infusion records and corresponding BG levels, each of which were recorded every 5 minutes.
An automated algorithm for analysis of these data was developed as follows: Prior to identification by the algorithm of onset of counterregulation, each subject's parameters a and b were estimated, along with a maximumlikelihood estimate of their initial BG level. This was accomplished by a modified Gauss-Newton nonlinear least squares parameter estimation algorithm (Johnson and Frasier, 1985;Straume, Frasier-Cadoret and Johnson, 1991) in which the differential equation [A] was integrated numerically for BG(t) by a fourth-order Runge-Kutta method (Press et al., 1989). Applied to each subjects' data set, this procedure successfully evaluated individually for each subject these characteristics of the dynamics of dextrose utilization and BG elimination during euglycemia and descent into hypoglycemia prior to onset of counterregulation.
The algorithm was initialized to consider first only those time points in which there was clearly no potential for counterregulation i.e., (those points comprising the euglycemic phase of the experiment; Phase 1 of Figure 1). The parameters a and b and the initial BG level were nonlinear least squares estimated to this subset of data, followed by evaluation of the standard deviation of fit to the BG data. At this point, the difference between the observed BG level for the next time point and the BG level predicted by the model in the absence of counterregulation was evaluated. If this difference was less than two standard deviations of fit, this next observed BG level was considered to be prior to onset of counterregulation and was included as an additional point for estimation of the parameters a and b and an initial BG level by equation [A]. This process was repeated iteratively until a BG level was identified that differed from the predicted value by more than two standard deviations of fit, indicating onset of counterregulation. From this point onward, the parameters a and b and the estimated initial BG level were fixed, onset of counterregulation (parameter T I ) was defined as the time of the previous time point, and the model began fitting to differential equation [B]. The process again proceeded iteratively one point at a time until either the remainder of the data set was successfully considered or until the need for a second component of the counterregulatory response was identified (in the same manner as above).

RESULTS
The average BG level during the first euglycemic phase (Phase 1) of the study was 6.3 mM (SEM = 0.1). The average nadir of BG reached during the second phase of the study was 2.5 mM (SEM = 0.08). The average plasma epinephrine concentration during Phase 1 was 52 (SEM = 6.3). During Phase 2 the average epinephrine peak was 367 (SEM = 44).

Goodness-of-fit of the Model
The algorithm was applied to the data sets of each of the 40 subjects. The accuracy of the data fit was tested by the coefficient of determination, usually interpreted as the percentage of the total variation of the dependent variable around its mean that is explained by the fitted model (Kvalseth, 1985), and by the mean square error (MSE) per data point. Since the model is nonlinear and a standard ANOVA p-value cannot be computed, the goodness-of-fit of each model was evaluated by the closeness of its coefficient of determination to 100% (note that our model is intrinsically nonlinear, and therefore the usual F-statistic and its significance level cannot be used). The median coefficient of determination across subjects was 97.7947, with a maximum of 99.7% and a minimum of 86.5%. The median MSE of the model fits was 0.16 mM, with a range from 0.07 to 0.44 mM, i.e. for all subjects the model-predicted BG fluctuations were within 0.44 mM from the observed BG values, with half of the subjects within 0.16 mM. This indicates an extremely good model fit for all subjects. Table I presents goodness-of-fit data for all subjects ordered by their coefficients of determination. Seven subjects had coefficients of determination above 99%, while only 2 subjects had coefficients of determination below 90%.
In order to better illustrate our model Figure 3 presents a sample (not the best but above average) data fit for subject #9 whose coefficient of determination was 98.8%, MSE = 0.14 mM. (thin black linc). The dcxtrose infusion is plotted along the right y-axis as a stepwise p r q line.
The x-axis of Figure 3 is the elapsed time in minutes. The left y-axis (BG (mM)) refers to three variables: 1) The BG data plotted by circles.
2) The full model fit, based on equation [B] (thick black line).
3) The model-predicted BG decay if counterregulation did not occur, based on equation [A] (thin black line). The dextrose infusion that is plotted along the right y-axis, is represented by a stepwise grey line, is a constant within each 5-minute interval and is adjusted every 5 minutes. The figure includes this subject's coefficient of determination R~ = 98.8% and MSE = 0.14 mM.

Counterregulation
The rate of counterregulation of each subject was estimated in units equivalent to mgkgtmin dextrose infusion, on the basis of the difference between models [B] and [A], as explained above (Parameter Estimation). In Figure 3 the counterregulation would be equivalent to dextrose infusion needed to elevate the BG level fit from the thin to the thick black line, i.e. from model [A] to model [B]. Consequently the onset of counterregulation will be the point where these two lines split, i.e. shortly after minute 90 in Figure 3.
An external validation of the predicted counterregulation rate was done using its correlations with the corresponding epinephrine data for each subject. Table I1 presents the correlation coefficients between model-predicted counterregulation and logarithm of epinephrine concentration together with their significance level. For subjects who counterregulated the median correlation coefficient was 0.82, range from 0.47 to 0.97 with all correlations significant at p = 0.05. The model estimated that four subjects did not counterregulate and this was confirmed by their non-increasing plasma epinephrine concentrations (#17, #36, #37 and #39 in Table 11). Figure 4 illustrates the relationship between the model-predicted counterregulation and the plasma epinephrine concentrations, recorded every 10 minutes, for subject #9 whose data and model curves were plotted in Figure 3. As with Figure 3, the xaxis represents the elapsed time in minutes. The left y-axis refers to logarithm of epinephrine concentration that is measured every 10 minutes and presented by a grey stepwise line. The right y-axis represents the model-evaluated counterregulation in units equivalent to mglkglmin dextrose infusion. The counterregulatory response of this subject began shortly after minute 90 and then increased rapidly with a small setback at minute 105. The correlation between epinephrine and counterregulation was significant, R = 0.91, p < 0.001.

Parameters of the Model
As we discussed in the previous section, our model has four essential parameters. Two of them, glucose   conversion a and insulin utilization coefficient b as evidenced by the negative correlation of the were determined and fixed during the euglycemic insulin utilization coefficient with the nadir of BG, Phase 1 of the study. The other two, time of counter-R = -0.51(p < 0.001). Lower nadir of BG was regulation onset T I and the counterregulation slope associated with higher counterregulation rate as eviat onset C 1 . r:, were determined during Phase 2. denced by the negative correlations between nadir of In addition, based on the model we computed the BG and the average and maximum counterregulation average and maximal counterregulation rate for each rates, R = -0.63 ( p < 0.00 1 ) and R = -0.64(p < subject as well as his BG level at onset of counter-0.001) respectively. The average epinephrine resregulation. Table I11 presents descriptive character-ponse per subject correlated with the average counistics for these parameters.
terregulation rate, R = 0.4, p = 0.005, while the Faster insulin utilization was associated with maximum epinephrine response correlated with the lower observed nadir of BG during the study maximal counterregulation, R = 0.4, p = 0.006.

Group Effects
The two groups of subjects, with and without a history of SH, did not differ in terms of age, duration of diabetes, insulin units/lulogram per day or glycosylated hemoglobin. During the test subjects from both groups reached similar nadirs of BG and had similar average epinephrine responses. Subjects with no history of SH demonstrated higher maximal epinephrine response, 493 pglml (SEM = 72) vs. 289 pglml (SEM = 51), p = 0.025. Table IVA summarizes these observations. Based on the model, the two subject groups did not differ in terms of glucose conversion and insulin utilization coefficients, time and BG at onset of counterregulation, and average counterregulation rate. Subjects with no history of SH demonstrated greater countesregulation slope at onset, i.e. faster onset of counterregulation, 2.4 (SEM = 0.6) vs. 1.1 (SEM = 0.2), p = 0.026 and marginally higher maximal counterregulation rate, p = 0.13. Table IVB summarizes these results.

DISCUSSION
The deterministic differential equation model developed in the present study accounts, in a highly reliable manner, for the dynamics of both exogenous dextrose-infusion-rate-dependent changes in blood glucose levels as well as endogenous physiological countesregulation during euglycemic hyperinsulinemic clamping of patients with IDDM. During euglycemia in the absence of counterregulatory response, the model i~, mechanistically parameterized, individually and separately for each patient considered, in terms of two physiological processes: 1) the dextrose-to-blood glucose conversion efficiency and 2) the insulin utilization efficiency for elimination of blood glucose. The model was implemented during analysis to use an objective criterion for identifying the time of onset of physiological counterregulation at low blood glucose levels. Counterregulation was then parameterized, again, individually and separately for each patient considered, in terms of the rate and volume of counterregulatory response, as well as for the potential for bimodal counterregulation.
The objective analysis performed by this implementation of the model was successful in all 40 of the blood glucose-dextrose infusion data sets of IDDM patients considered in this study, typically accounting for greater than 95% 01 the observed variance in blood glucose time series. Additionally, the counterregulatory responses predicted by the model are consistent with observed plasma epinephrine concentrations, as indicated by the high correlation between the two. Interestingly, the modeling results suggest the previously unrecognized possibility that blood glucose counterregulation may be a multicomponent process (as a bimodal counterregulatory response was predicted for 22 of the 40 individuals examined).
Lower nadir of BG during the study was associated by our model with a faster insulin utilization. On the other hand, a lower nadir of BG prompted higher and more aggressive counterregulatory response, but did not result in clearly larger epinephrine response as evidenced by a non-significant nadir BG-epinephrine correlation R = -0.2, p = 0.1.
Finally, our data indicated that a history of SH was associated with less aggressive epinephrine response to low BG, while our model demonstrated that this effect is primarily due to a less aggressive onset of counterregulation and only partly due to lower maximal counterregulation response. This result refines research findings that associate risk for SH with "impaired glucose counterregulation" (Cryer and Gerich, 1985;Gerich, 1988;White et al., 1983) clarifying one dimension of this counterregulation impairment. As a result, our model may shed light not only on which subjects may be poor candidates for intensive insulin therapy because of vulnerability to SH; but also why they may be so (in functionally mechanistic terms).