- Email: [email protected]
Dynamic Control of Insulin on Glucose Kinetics: Tracer Experiment Design and Time-Varying Interpretative Models
CO)l\Tig ht © IF.-\C IlIth Tri e llllial \\'" rld C "" g ITSS , \(ullit' h , FRC , I ~ 1 8 i
DYNAMIC CONTROL OF INSULIN ON GLUCOSE KINETICS: TRACER EXPERIMENT DESIGN AND TIME-VARYING INTERPRETATIVE MODELS C. Cobelli*, A. Mari*, G. Toffolo*, A. D. Cherrington** and O. P. McGuinness** *Depa)'tlllellt of Elert)'oll ies allr/ IlIfo ),lIIatifs, U ll il'l'I'sity of Par/o1'(/ , all r/ C \'R IlI stitlltl' jill' 5.1',1/ 1' 111 DY/(/lIIin allt! Biul'IIgilll'l')'illg, PIlt/m 'a, Ita l\, ** Dl'pa )'tllll'll t of Physiulo/",,)', \ 'al/(ll'I'bilt L'lI h'l'l'Sity, S ash"il/I' , TI'III1I',ISI:I' , [ 'SA
Abstract. In this paper a linear time-varying two-compartment model has been developed for studying glucose kinetics from tracer experiments in the non-steady state. Parameters have been estimated using data from a pilot experiment performed in a dog. The model has been used to compute the time-course of glucose production, which is not directly measurable. More traditional non-steady-state methods were also comparatively evaluated. The results suggest that the effect of insulin on glucose production is slower and less marked than currently thought. The inadequacy of the standard simplified approaches is confirmed, in agreement with a recently published theoretical work. Keywords. Medical systems; Modelling; System analysis; Identification; Parameter estimation; Tracer experiment .
traditional limitations.
INTRODUCTION The use of radioactive tracers has made possible basic advances in the knowledge of the kinetics of glucose and of the mechanisms by which its concentration is regulated in the body. The effects of insulin per se on glucose kinetics can be studied with tracers by clamping glucose concentration at its basal value and elevating insulin levels. (DeFronzo, Tobin and Andres, 1979). In the steady state, i.e . both at the basal and at the hyperinsulinemic state, the problem has been recently studied by Cobelli, Toffolo and Ferrannini (1984) and Ferrannini and others (1985). In the non-steady state, i.e. during the transition from the basal to the normoglycemic hyperinsulinemic steady state, the evaluation of the effect of insulin on glucose production and utilization is much more critical from a modeling point of view, as discussed by Cobelli, Mari and Ferrannini (1987). In that paper it has been shown that widely used simplified models, such as Steele's model (Steele, 1959), can be not sufficently accurate. A linear time-varying compartmental model, which incorporates the knowledge of the system in the initial and final steady state, has been proposed for a more accurate non-steady-state analysis. The present work improves that analysis, which was performed on average data, by including in the experimental protocol the study of the individual steady-state kinetics. Furthermore, the protocol make possible to evaluate experimentally the role of the specific activity, which has been shown theoretically (Cobelli, Mari and Ferrannini, 1987) to be a determinant of the error of any model. The new approach has been also compared with more
Time-varying
systems;
evidencing
ones,
their
EXPERIMENTAL PROTOCOL A pilot experiment has been performed on a 21 kg healthy dog, using radioactively 3
labeled glucose ([3- HJglucOse) as a tracer of the native substrate (the tracee). At the start of the experiment (time 0), an impulsive dose of tracer (94.55.10 6 dpm) was injected and followed by a constant infusion (781456 dpm/min). Tracer concentration was then monitored for 150 min to assess basal glucose kinetics. At time 150 min the basal steady state was perturbed by infusing insulin at a constant rate (5 mU·min-1.kg- 1 ) and keeping glucose concentration constant through a variable glucose infusion controlled by a feedback loop (glucose clamp technique, see e.g. DeFronzo, Tobin and Andres, 1979). Throughout the rest of the experiment insulin was infused at the same rate and glucose was clamped at its basal values. Tracer infusion rate was left unchanged until time 330 min, and in this period tracer concentration falls due to the effect of insulin. At time 330 min, labeled glucose was infused at a rate which exactly parallels that of the glucose infused for the clamp. The ratio of tracer infusion rate to cold glucose infusion rate was in fact kept constant to 21015 dpm/mg. From time 330 min to 480 min tracer concentration rises again due to the increased tracer infusion rate. At time 480 min tracer infusion was stopped, and the washout concentration curve monitored until 540 min to assess glucose kinetics in the new normoglycemic, hyperinsulinemic steady state.
25
c . Cohelli
1'1
1/1.
MODEL DESCRIPTION Tracer and tracee dynamics are described by the same linear time-varying system. Outputs are tracer and tracee concentrations and inputs are tracer infusion rate and tracee rate of appearance, i.e. exogenous glucose infusion rate plus endogenous production. The system structure in our experimental conditions must reflect the fact that the protocol realizes a transition between two steady states. The system matrix is thus a matrix whose entries are constant before the starting of the perturbation (time to) and then vary with time to reach asymptotically a new constant value. Thus, the system can be described by A(t)x(t) + Bu(t) y( t)
Cx( t)
A(t)
[a . . ( t ) ]
(1)
Fig. 1.
The time-varying two-compartment model for tracer and tracee kinetics. Inputs and outputs are specified in the text.
=
We have assumed P01(t) 0 because kOl(t), which represents the insulin-independent glucose utilization, is constant, as discussed in Cobelli, Toffolo and Ferrannini (1984). CALCULATION OF GLUCOSE PRODUCTION
1)
with o
lim t-++oo
a ..
1)
We will assume that the system has a compartmental structure (Carson, Cobelli and Finkelstein, 1983) and inputs and outputs take place in the first compartment. The system matrices Band C are thus B = [10 .. 0]', C = [cO .. O], where c is the reciprocal of the first compartment volume. The entries of the n by n matrix A(t) are
The calculation of glucose production from concentration measurements is an input estimation problem for the system of equation 1. In our experimental conditions, the input is the sum of the glucose flux produced endogenously and that infused exogenously to clamp glucose concentration. Since we are interested in endogenous glucose production, g(t), we have evaluated it by describing its time-course with a parametrized function of time and by parameter estimation. The chosen function is g(t)
kij(t)
( 2)
i*j
and it is of employed for parameters .
i=n a . . (t) JJ
(4)
- E kij(t) i=o
1*)
the the
same kind of that time-varying system
MODEL IDENTIFICATION
with ( k~ .
Identifiability Analysis
t
o i 1) 00 0 ~ k~ . + (kij-kij)Pij(t-to) 1)
t)t -
0
where Pij(t) are functions which range from zero, when t=O, to one, when t approaces infinity. The choice of the system order stems from the analysis of the system in steady state. A two-compartment model has been shown adequate for glucose kinetics in previous studies (Cobelli, Toffolo and Ferrannini, 1984; Radziuk, Norwich and Vranich, 1978), and has been adopted here (Fig. 1). The functions Pij(t) represent the control of insulin on glucose kinetics. We have chosen a simple expression for these functions, i.e. -\I . .
1 - e
1)
t ij
21,12,02
(3 )
Identifiability of system parameters, in the deterministic sense, can be tested in various ways if the system is linear and time-invariant (see e.g. Carson, Cobelli and Finkelstein, 1983, chapter 7) . When the system is time-varying, the methods based on the transfer function fail, while those based on the canonical representation of the system matrices can still be employed under rather general conditions. In the time-varying case, the time-course of the parameters of some canonical form (e.g. the controllability canonical form) is uniquely determinable from the knowledge of the impulse response. Thus, if the original system is transformed by equivalence transformation (Silverman, 1971) into a canonical form, it is possible to evaluate identifiability by examining if the parameters of the original system, as they appear in the canonical form, are (locally or globally) determinable from the knowledge of the time-varying parameters of the canonical form . If, for instance, the mapping of the original parameters into the canonical ones is one-to-one, the system is uniquely
27
Dyna mic Contro l o f Insulin o n Glu cose Kin e tics
identifiable. Although this approach is in general more complex in the time-varying case because the equivalence tranformation equations are no longer algebraic but differential, it can still work for system of low dimension, as in this study . Parameter Estimation
!i
_, ~i 614
~~
I~
....
2 0 ~----------------~
!i
Parameter estimation has been performed by least squares. It is assumed that the output of the system of equation 1 is corrupted by a white noise of constant variance. The minimization routine employed is based on the Marquardt algorithm (VBOIA, Harwell subroutine library) and requires the derivatives of the system output with respect to the parameters which have been computed by integration of the sensitivity equations. The approximation of the Hessian obtained from the gradient has been used to compute the asymptotic covariance matrix of the parameters. The noise variance has been estimated from the residuals . The parameters of the model of equations 1-3 have been estimated from tracer data, and, with these estimates, the parameters of the glucose production function (equation 4) have been subsequently estimated from tracee data.
~~ec 300 200
tl a,100 ~.
I
0
..
~
0.8
COMPARISON WITH TRADITIONAL ANALYSES Glucose production has been also estimated by means of the standard tools of steady-state analysis and by Steele equation in non-steady state (Steele 1959; see also Radziuk, Norwich and Vranich, 1978, equation 1). Steady-state analysis is performed directly on the observed data averaged over the intervals 120-150 min, 300-330 min, and 450-480 min. Data smoothing to compute the derivatives of tracer and tracee concentration required in Steele equation has been obtained by the predictions of the time-varying model. The pool fraction used in Steele equation was 0.65 and the total distribution volume 250 ml/kg (cf. e.g. Radziuk, Norwich and Vranich, 1978).
The system parameters estimated by least squares are shown in Table 1, while Fig . 2 shows the model fit, and Fig. 3 the residuals, i.e. the differences between
0 '60 '120'180'240'300'360'420'480'540 IIIln
Fig. 2 .
~~~
~~
....
RESULTS Time-varying Model
_ _ _---J
....
Upper part : tracer and tracee infusion rates. Lower part: tracer and tracee concentrations . Circles are observed data and solid lines model predictions.
-10' 3 2 1 0 -1 -2 -3
~.~:~l
~'O
~~0.1 ~
,A" ....... TV\,.... ....
.A
"
."tf\_ 1\ VV'''v--y
Si -0. 2 , , , , , , , , , , , , , , , , , , , ....
0
60 120 180 240 300 360 420 480 540 IIIln
TABLE l
Fig. 3. Tracer and tracee residuals.
Model Parameters
Parameter
k2l
k12
k02
0 initial value (k ) 0.0656 0.0492 0.00292 ij 00 (k ) 0.111 0.00618 0.0497 final value ij time constant (\I . . ) 0.637 0.0571 0.00269
1)
Symbols are defined in equations 2-3.
All 1
the parameters are given in minA 1 value of 0.01 min- has been assumed for k . The value of c in the C vector Ol (equation 1) is 0.000387 ml- 1 •
predicted and observed values. The standard deviation of the error estimated 3
from the residuals was 1.18.10 dpm/ml. Some of the coefficients of variation were above 100 %, but it was not possible to improve precision by grouping some of the time constants \/ .. without a loss of fit. 1)
On the other hand, several trials with a more complex structure of the p .. (t) did 1) not lead to a significant decrease of the residual sum of squares nor yielded a more satisfactory pattern of the residuals.
C. CoiJclli 1"1 (/1.
TABLE 2
Glucose Production The parameters gl and ~ of the glucose production function (equation 4) could not be estimated with an acceptable coefficient of variation, and were remarkably correlated. The function g(t) was thus simplified by assuming that the asymptotic production is zero, i.e. go = gl' and the parameters could be reliably estimated without a loss of fit. The parameters were in this case go = 39.7 mg/min, ~ = 0.00114 min- , and the coefficients of variation were below 40 %. Fig. 2 shows the model fit and Fig. 3 the residuals. The estimated standard deviation of the error was 0.063 mg/ml. The time course of glucose production is shown in Fig. 4. It can be seen that glucose production, although affected by insulin action, is suppressed only to a modest extent.
Glucose Production 450-480
120-150
300-330
Steady-state analysis
44.4
-8.2
30.2
Steele equation
42.2
12.7
29.7
Time-varying model
39.7
29.7
27.7
Interval
1
i Q
c
E
:~ 1 20
~, 0 ~ 0. :g E -20
. . /\
i.. ------------~----~:~..-~...~:~:-------------,. .:....... ..
(.)
~
-40 0""-.-6'T'0......... ' t""2'""0""'-tT'e-0'''-2"'T;'-0T''3-0"-'0...,'....3... 6-0,...'4""2-0'T'-4T'eo""''--'s;'0 mln
Fig. 4.
Glucose production computed with the time-varying model (solid line) and with Steele equation (dotted line).
Steele Equation The time course of glucose production obtained by Steele equation (Fig. 4) is dramatically different from that of the time-varying model. It takes on physically unrealizable negative values for long periods and shows also unlikely rises above the baseline. It must be pointed out that the predicted production shows a jump (down to less than -400 mg/min and up to about 63 mg/min) when the tracer infusion is drastically increased (time 330 min, Fig. 2). The reasons of this behaviour are explained by the sudden change of the derivative of the specific activity and have been treated in detail in the paper by Cobelli, Mari and Ferrannini (1987). The irregularities of the profile are due to small discontinuities in the first derivative of the tracee concentration, which is that predicted by the time-varying model (Fig. 2). The profile computed assuming a constant value of 1.09 mg/ml for the tracee concentation, i.e. the mean value during the clamp, was by far more erratic with an initial sudden increase to a rate of over 160 mg/min and several negative nadirs of less than -40 mg/min (data not shown) . Steady-state Analysis The results of steady-state analysis are shown in Table 2. The steady-state production values are compared with the average values over corresponding intervals obtained by Steele equation and the time-varying model. It clearly
Values (mg/min) are averages corresponding 30 min interval.
over
the
appears that while in the periods 120-150 min and 450-480 min the computed values are similar, in the period 300-330 min, i.e. at the end of the 3rd hr of the clamp, they differ markedly. Steady-state production is in fact remarkably negative, and Steele's production is almost one third of that of the time varying model. DISCUSSION System identification for linear time-varying systems is a difficult problem when the structure of the time-varying parametrs is not known in advance. Assumptions are in fact needed to reconstruct the time-course of these parameters, the correctness of which can rarely be tested. We have assumed a very simple structure of k .. (t), i.e. an l.) exponential transition between two constant values, bearing in mind that insulin exerts its control on the parameters of the glucose system from compartments remote from plasma (e.g. Bergman and others, 1979; Insel and others, 1975). This simple assumption allowed us to estimate the whole set of parameters with a reasonable precision, and to obtain a plausible description of the data (Fig. 2 and Fig. 3). The variance of the error estimated from the residuals of the tracer concentration was however higher than that of the pure assay error. Although it is expected that the total measurement error is larger than the pure assay error, it is likely that some model error is still pr2sent in our analysis. However, we were not able from this single case to develop a more satisfactory model, in spite of several attempts with different and more complex structures of the time-varying parameters. The calculation of glucose production is an input estimation problem which, being highly ill-conditioned, requires some regularization technique. The choice of a prescribed function of time for the production seemed to us a straightforward choice, because it is natural to assume that the control exerted by insulin on glucose production is similar to that exerted on the system parameters. The goodness of this assumption can be evaluated from the ability of the whole model to predict tracee concentratl.on (Fig. 2 and Fig. 3). Since the variance of the residuals appears to be remarkably higher than the assay error, it can be concluded that the estimated time-course of glucose rroduction (Fig. 4) cannot
lhll;tlllil ( :ollln,1 of IllslIlill Oil (;llIlosl' hilll'lilS
account for all the changes in glucose concentration. The scatter of the points around the model predicted concentration suggests that the production may not have a smooth t\me-course. Alternatively, the presence in the system of a third, small and fast equilibrating compartment may be an other explanation of the rapid concentration variations observed in response to the frequent changes in exogenous glucose infusion rates. The deep undershoot observed between 150 and 180 min remains also unexplained. In this period, to achieve a better fit the production should have increased considerably above the basal level. Although this unlikely event cannot be excluded, it seems more reasonable to explain this phenomenon with a defect in the time-varying model, rather than in the production function. Further investigation and experiments are needed to clarify this issue. The traditional non-steady-state analysis with Steele equation appears to be completely inadequate. Apart from the problem of ill-conditioning of the approach (cf. Cobelli and Mari, 1987), the method, which is based on a monocompartmental approximation, is by far too simplistic in this experimental situation. Physically unrealizable negative values are obtained between 180 and 270 min, and an unphysiological jump is observed at 330 min. The latter phenomenon demonstrate a clear dependence of the result on the time-course of the specific activity, in full accord with the theory developed by Cobelli, Mari and Ferrannini (1987). The results obtained by the steady-state analysis clearly show that a false assumption of steady state may lead to erroneous results. While in fact in the basal period (120-150 min) and in the last period of the clamp (450-480 min) all the calculations are in agreement, they differ markedly between 3CO and 330 min (Table 2). Our analysis thus indicates that the negative values of glucose production reported in recent publications which deal with similar experiments (Bell, Firth and Rizza, 1986; Prager, Wallace and Olefsky, 1986), can be entirely due to a model error. Finally, it is worthwhile to emphasize that the modeling approach here employed, led to a although perfectible, has conclusion which challenge the current physiological thinking, i.e. a rapid and complete inhibitory effect of insulin on glucose production. Even if based on a single case, our results suggest that glucose production is slowly and possibly not completely suppressed by insulin. CONCLUSION Cobelli, Mari and Ferrannini (1987) have already emphasized the importance of proper modeling in the non-steady state. They proposed a two-compartment time-varying model for data analysis which includes the most recent physiological knowledge and overcomes the limitations of previously published models. This paper improves that analysis, which was performed on average data, showing that the model can be identified from
experimental tracer uata in a single individual. Furthermore, glucose production has been calculated as an input estimation problem by taking into proper account the small glucose fluctuations. The results of our analysis indicate that the effect of insulin on glucose production is less rapid and marked than currently thought. The study also indicates, in full agreement with the theory (Cobelli, Mari and Ferrannini, 1987), that oversimplified approaches, such as Steele equation or improper steady-state analysis, are not reliable for computing glucose production. REFERENCES Bell, P. M., R. G. Firth, and R. A. Rizza (1986). Assessment of insulin action in insulin-dependent diabetes mellitus
using
[6
14
3
CJglucose, [3 HJglucose,
3
and [2 HJglucose. Differences in the apparent pattern of insulin resistance depending on the isotope used. ~ Clin. Invest. 78:1479-1486. Bergman, R. N., Y. z. Ider, C. R. (1979). C. Cobelli Bowden, and estimation of insulin Quantitative Physiol. sensitivity. Am. ~ 236 :E66 7-E6 77. R., C. Cobelli, and L. Carson, E. (1983). The Mathematical Finkelstein Modeling of Metabolic and Endocrine Systems. John Wiley, New York. Cobelli, C., G. Toffolo, and E. Ferrannini (1984). A model of glucose kinetics and their control by insulin, compartmental and noncompartmental approaches. Math. Biosci. 72:291-315. Cobelli, C., and A. Mari (1987). Models for the estimation of glucose fluxes in non-steady state from concentration measurements. Preprints 10th IFAC World Congress. Cobelli, C., A. Mari, and E. Ferrannini (1987). Non-steady state: error analysis of Steele's model and developments for glucose kinetics. Am. ~ Physiol._ (in press) . Tobin, and R. DeFronzo, R. A., J. D. Andres (1979). Glucose clamp technique: a method for quantifying insulin secretion ~ Physiol. and resistance. Am. 237 :E214-E223. Ferrannini, E., J. D. Smith, C. Cobelli, G. Toffolo, A. Pilo, and R. A. DeFronzo (1985). Effect of insulin on the distribution and disposition of glucose in man. ~ Clin. Invest. 76:357-364. Insel, P. A., J. E. Liljenquist, J. D. Tobin, R. S. Sherwin, P. Watkins, R. Andres, and M. Berman (1975). Insulin control of glucose metabolism in man. A new kinetic analysis. ~ Clin. Invest. 55:1057-1066. Prager, R., P. Wallace, and J. M. Olefsky (1986) . In vivo kinetics of insulin action on peripheral glucose disposal and hepatic glucose output in normal and obese subjects. ~ Clin. Invest. 78:472-481.
30
C. Cobelli
Radziuk, J., K. H. Norwich, and H. Vranic (1978). Experimental validation of measurements of glucose turnover in nonsteady state. Am . ~ Physiol. 234:E84-E93. (1971). Silverman, L. H. Realization of linear dynamical systems. IEEE Trans. Automat. Contr. 16:554-567 .
Steele, R (1959). Influences of glucose loading and of injected insulin on hepatic glucose output. Ann. NY Acad. ScL 82:420-430.
PI
al.
DYNAMIC CONTROL OF INSULIN ON GLUCOSE KINETICS: TRACER EXPERIMENT DESIGN AND TIME-VARYING INTERPRETATIVE MODELS C. Cobelli*, A. Mari*, G. Toffolo*, A. D. Cherrington** and O. P. McGuinness** *Depa)'tlllellt of Elert)'oll ies allr/ IlIfo ),lIIatifs, U ll il'l'I'sity of Par/o1'(/ , all r/ C \'R IlI stitlltl' jill' 5.1',1/ 1' 111 DY/(/lIIin allt! Biul'IIgilll'l')'illg, PIlt/m 'a, Ita l\, ** Dl'pa )'tllll'll t of Physiulo/",,)', \ 'al/(ll'I'bilt L'lI h'l'l'Sity, S ash"il/I' , TI'III1I',ISI:I' , [ 'SA
Abstract. In this paper a linear time-varying two-compartment model has been developed for studying glucose kinetics from tracer experiments in the non-steady state. Parameters have been estimated using data from a pilot experiment performed in a dog. The model has been used to compute the time-course of glucose production, which is not directly measurable. More traditional non-steady-state methods were also comparatively evaluated. The results suggest that the effect of insulin on glucose production is slower and less marked than currently thought. The inadequacy of the standard simplified approaches is confirmed, in agreement with a recently published theoretical work. Keywords. Medical systems; Modelling; System analysis; Identification; Parameter estimation; Tracer experiment .
traditional limitations.
INTRODUCTION The use of radioactive tracers has made possible basic advances in the knowledge of the kinetics of glucose and of the mechanisms by which its concentration is regulated in the body. The effects of insulin per se on glucose kinetics can be studied with tracers by clamping glucose concentration at its basal value and elevating insulin levels. (DeFronzo, Tobin and Andres, 1979). In the steady state, i.e . both at the basal and at the hyperinsulinemic state, the problem has been recently studied by Cobelli, Toffolo and Ferrannini (1984) and Ferrannini and others (1985). In the non-steady state, i.e. during the transition from the basal to the normoglycemic hyperinsulinemic steady state, the evaluation of the effect of insulin on glucose production and utilization is much more critical from a modeling point of view, as discussed by Cobelli, Mari and Ferrannini (1987). In that paper it has been shown that widely used simplified models, such as Steele's model (Steele, 1959), can be not sufficently accurate. A linear time-varying compartmental model, which incorporates the knowledge of the system in the initial and final steady state, has been proposed for a more accurate non-steady-state analysis. The present work improves that analysis, which was performed on average data, by including in the experimental protocol the study of the individual steady-state kinetics. Furthermore, the protocol make possible to evaluate experimentally the role of the specific activity, which has been shown theoretically (Cobelli, Mari and Ferrannini, 1987) to be a determinant of the error of any model. The new approach has been also compared with more
Time-varying
systems;
evidencing
ones,
their
EXPERIMENTAL PROTOCOL A pilot experiment has been performed on a 21 kg healthy dog, using radioactively 3
labeled glucose ([3- HJglucOse) as a tracer of the native substrate (the tracee). At the start of the experiment (time 0), an impulsive dose of tracer (94.55.10 6 dpm) was injected and followed by a constant infusion (781456 dpm/min). Tracer concentration was then monitored for 150 min to assess basal glucose kinetics. At time 150 min the basal steady state was perturbed by infusing insulin at a constant rate (5 mU·min-1.kg- 1 ) and keeping glucose concentration constant through a variable glucose infusion controlled by a feedback loop (glucose clamp technique, see e.g. DeFronzo, Tobin and Andres, 1979). Throughout the rest of the experiment insulin was infused at the same rate and glucose was clamped at its basal values. Tracer infusion rate was left unchanged until time 330 min, and in this period tracer concentration falls due to the effect of insulin. At time 330 min, labeled glucose was infused at a rate which exactly parallels that of the glucose infused for the clamp. The ratio of tracer infusion rate to cold glucose infusion rate was in fact kept constant to 21015 dpm/mg. From time 330 min to 480 min tracer concentration rises again due to the increased tracer infusion rate. At time 480 min tracer infusion was stopped, and the washout concentration curve monitored until 540 min to assess glucose kinetics in the new normoglycemic, hyperinsulinemic steady state.
25
c . Cohelli
1'1
1/1.
MODEL DESCRIPTION Tracer and tracee dynamics are described by the same linear time-varying system. Outputs are tracer and tracee concentrations and inputs are tracer infusion rate and tracee rate of appearance, i.e. exogenous glucose infusion rate plus endogenous production. The system structure in our experimental conditions must reflect the fact that the protocol realizes a transition between two steady states. The system matrix is thus a matrix whose entries are constant before the starting of the perturbation (time to) and then vary with time to reach asymptotically a new constant value. Thus, the system can be described by A(t)x(t) + Bu(t) y( t)
Cx( t)
A(t)
[a . . ( t ) ]
(1)
Fig. 1.
The time-varying two-compartment model for tracer and tracee kinetics. Inputs and outputs are specified in the text.
=
We have assumed P01(t) 0 because kOl(t), which represents the insulin-independent glucose utilization, is constant, as discussed in Cobelli, Toffolo and Ferrannini (1984). CALCULATION OF GLUCOSE PRODUCTION
1)
with o
lim t-++oo
a ..
1)
We will assume that the system has a compartmental structure (Carson, Cobelli and Finkelstein, 1983) and inputs and outputs take place in the first compartment. The system matrices Band C are thus B = [10 .. 0]', C = [cO .. O], where c is the reciprocal of the first compartment volume. The entries of the n by n matrix A(t) are
The calculation of glucose production from concentration measurements is an input estimation problem for the system of equation 1. In our experimental conditions, the input is the sum of the glucose flux produced endogenously and that infused exogenously to clamp glucose concentration. Since we are interested in endogenous glucose production, g(t), we have evaluated it by describing its time-course with a parametrized function of time and by parameter estimation. The chosen function is g(t)
kij(t)
( 2)
i*j
and it is of employed for parameters .
i=n a . . (t) JJ
(4)
- E kij(t) i=o
1*)
the the
same kind of that time-varying system
MODEL IDENTIFICATION
with ( k~ .
Identifiability Analysis
t
o i 1) 00 0 ~ k~ . + (kij-kij)Pij(t-to) 1)
t)t -
0
where Pij(t) are functions which range from zero, when t=O, to one, when t approaces infinity. The choice of the system order stems from the analysis of the system in steady state. A two-compartment model has been shown adequate for glucose kinetics in previous studies (Cobelli, Toffolo and Ferrannini, 1984; Radziuk, Norwich and Vranich, 1978), and has been adopted here (Fig. 1). The functions Pij(t) represent the control of insulin on glucose kinetics. We have chosen a simple expression for these functions, i.e. -\I . .
1 - e
1)
t ij
21,12,02
(3 )
Identifiability of system parameters, in the deterministic sense, can be tested in various ways if the system is linear and time-invariant (see e.g. Carson, Cobelli and Finkelstein, 1983, chapter 7) . When the system is time-varying, the methods based on the transfer function fail, while those based on the canonical representation of the system matrices can still be employed under rather general conditions. In the time-varying case, the time-course of the parameters of some canonical form (e.g. the controllability canonical form) is uniquely determinable from the knowledge of the impulse response. Thus, if the original system is transformed by equivalence transformation (Silverman, 1971) into a canonical form, it is possible to evaluate identifiability by examining if the parameters of the original system, as they appear in the canonical form, are (locally or globally) determinable from the knowledge of the time-varying parameters of the canonical form . If, for instance, the mapping of the original parameters into the canonical ones is one-to-one, the system is uniquely
27
Dyna mic Contro l o f Insulin o n Glu cose Kin e tics
identifiable. Although this approach is in general more complex in the time-varying case because the equivalence tranformation equations are no longer algebraic but differential, it can still work for system of low dimension, as in this study . Parameter Estimation
!i
_, ~i 614
~~
I~
....
2 0 ~----------------~
!i
Parameter estimation has been performed by least squares. It is assumed that the output of the system of equation 1 is corrupted by a white noise of constant variance. The minimization routine employed is based on the Marquardt algorithm (VBOIA, Harwell subroutine library) and requires the derivatives of the system output with respect to the parameters which have been computed by integration of the sensitivity equations. The approximation of the Hessian obtained from the gradient has been used to compute the asymptotic covariance matrix of the parameters. The noise variance has been estimated from the residuals . The parameters of the model of equations 1-3 have been estimated from tracer data, and, with these estimates, the parameters of the glucose production function (equation 4) have been subsequently estimated from tracee data.
~~ec 300 200
tl a,100 ~.
I
0
..
~
0.8
COMPARISON WITH TRADITIONAL ANALYSES Glucose production has been also estimated by means of the standard tools of steady-state analysis and by Steele equation in non-steady state (Steele 1959; see also Radziuk, Norwich and Vranich, 1978, equation 1). Steady-state analysis is performed directly on the observed data averaged over the intervals 120-150 min, 300-330 min, and 450-480 min. Data smoothing to compute the derivatives of tracer and tracee concentration required in Steele equation has been obtained by the predictions of the time-varying model. The pool fraction used in Steele equation was 0.65 and the total distribution volume 250 ml/kg (cf. e.g. Radziuk, Norwich and Vranich, 1978).
The system parameters estimated by least squares are shown in Table 1, while Fig . 2 shows the model fit, and Fig. 3 the residuals, i.e. the differences between
0 '60 '120'180'240'300'360'420'480'540 IIIln
Fig. 2 .
~~~
~~
....
RESULTS Time-varying Model
_ _ _---J
....
Upper part : tracer and tracee infusion rates. Lower part: tracer and tracee concentrations . Circles are observed data and solid lines model predictions.
-10' 3 2 1 0 -1 -2 -3
~.~:~l
~'O
~~0.1 ~
,A" ....... TV\,.... ....
.A
"
."tf\_ 1\ VV'''v--y
Si -0. 2 , , , , , , , , , , , , , , , , , , , ....
0
60 120 180 240 300 360 420 480 540 IIIln
TABLE l
Fig. 3. Tracer and tracee residuals.
Model Parameters
Parameter
k2l
k12
k02
0 initial value (k ) 0.0656 0.0492 0.00292 ij 00 (k ) 0.111 0.00618 0.0497 final value ij time constant (\I . . ) 0.637 0.0571 0.00269
1)
Symbols are defined in equations 2-3.
All 1
the parameters are given in minA 1 value of 0.01 min- has been assumed for k . The value of c in the C vector Ol (equation 1) is 0.000387 ml- 1 •
predicted and observed values. The standard deviation of the error estimated 3
from the residuals was 1.18.10 dpm/ml. Some of the coefficients of variation were above 100 %, but it was not possible to improve precision by grouping some of the time constants \/ .. without a loss of fit. 1)
On the other hand, several trials with a more complex structure of the p .. (t) did 1) not lead to a significant decrease of the residual sum of squares nor yielded a more satisfactory pattern of the residuals.
C. CoiJclli 1"1 (/1.
TABLE 2
Glucose Production The parameters gl and ~ of the glucose production function (equation 4) could not be estimated with an acceptable coefficient of variation, and were remarkably correlated. The function g(t) was thus simplified by assuming that the asymptotic production is zero, i.e. go = gl' and the parameters could be reliably estimated without a loss of fit. The parameters were in this case go = 39.7 mg/min, ~ = 0.00114 min- , and the coefficients of variation were below 40 %. Fig. 2 shows the model fit and Fig. 3 the residuals. The estimated standard deviation of the error was 0.063 mg/ml. The time course of glucose production is shown in Fig. 4. It can be seen that glucose production, although affected by insulin action, is suppressed only to a modest extent.
Glucose Production 450-480
120-150
300-330
Steady-state analysis
44.4
-8.2
30.2
Steele equation
42.2
12.7
29.7
Time-varying model
39.7
29.7
27.7
Interval
1
i Q
c
E
:~ 1 20
~, 0 ~ 0. :g E -20
. . /\
i.. ------------~----~:~..-~...~:~:-------------,. .:....... ..
(.)
~
-40 0""-.-6'T'0......... ' t""2'""0""'-tT'e-0'''-2"'T;'-0T''3-0"-'0...,'....3... 6-0,...'4""2-0'T'-4T'eo""''--'s;'0 mln
Fig. 4.
Glucose production computed with the time-varying model (solid line) and with Steele equation (dotted line).
Steele Equation The time course of glucose production obtained by Steele equation (Fig. 4) is dramatically different from that of the time-varying model. It takes on physically unrealizable negative values for long periods and shows also unlikely rises above the baseline. It must be pointed out that the predicted production shows a jump (down to less than -400 mg/min and up to about 63 mg/min) when the tracer infusion is drastically increased (time 330 min, Fig. 2). The reasons of this behaviour are explained by the sudden change of the derivative of the specific activity and have been treated in detail in the paper by Cobelli, Mari and Ferrannini (1987). The irregularities of the profile are due to small discontinuities in the first derivative of the tracee concentration, which is that predicted by the time-varying model (Fig. 2). The profile computed assuming a constant value of 1.09 mg/ml for the tracee concentation, i.e. the mean value during the clamp, was by far more erratic with an initial sudden increase to a rate of over 160 mg/min and several negative nadirs of less than -40 mg/min (data not shown) . Steady-state Analysis The results of steady-state analysis are shown in Table 2. The steady-state production values are compared with the average values over corresponding intervals obtained by Steele equation and the time-varying model. It clearly
Values (mg/min) are averages corresponding 30 min interval.
over
the
appears that while in the periods 120-150 min and 450-480 min the computed values are similar, in the period 300-330 min, i.e. at the end of the 3rd hr of the clamp, they differ markedly. Steady-state production is in fact remarkably negative, and Steele's production is almost one third of that of the time varying model. DISCUSSION System identification for linear time-varying systems is a difficult problem when the structure of the time-varying parametrs is not known in advance. Assumptions are in fact needed to reconstruct the time-course of these parameters, the correctness of which can rarely be tested. We have assumed a very simple structure of k .. (t), i.e. an l.) exponential transition between two constant values, bearing in mind that insulin exerts its control on the parameters of the glucose system from compartments remote from plasma (e.g. Bergman and others, 1979; Insel and others, 1975). This simple assumption allowed us to estimate the whole set of parameters with a reasonable precision, and to obtain a plausible description of the data (Fig. 2 and Fig. 3). The variance of the error estimated from the residuals of the tracer concentration was however higher than that of the pure assay error. Although it is expected that the total measurement error is larger than the pure assay error, it is likely that some model error is still pr2sent in our analysis. However, we were not able from this single case to develop a more satisfactory model, in spite of several attempts with different and more complex structures of the time-varying parameters. The calculation of glucose production is an input estimation problem which, being highly ill-conditioned, requires some regularization technique. The choice of a prescribed function of time for the production seemed to us a straightforward choice, because it is natural to assume that the control exerted by insulin on glucose production is similar to that exerted on the system parameters. The goodness of this assumption can be evaluated from the ability of the whole model to predict tracee concentratl.on (Fig. 2 and Fig. 3). Since the variance of the residuals appears to be remarkably higher than the assay error, it can be concluded that the estimated time-course of glucose rroduction (Fig. 4) cannot
lhll;tlllil ( :ollln,1 of IllslIlill Oil (;llIlosl' hilll'lilS
account for all the changes in glucose concentration. The scatter of the points around the model predicted concentration suggests that the production may not have a smooth t\me-course. Alternatively, the presence in the system of a third, small and fast equilibrating compartment may be an other explanation of the rapid concentration variations observed in response to the frequent changes in exogenous glucose infusion rates. The deep undershoot observed between 150 and 180 min remains also unexplained. In this period, to achieve a better fit the production should have increased considerably above the basal level. Although this unlikely event cannot be excluded, it seems more reasonable to explain this phenomenon with a defect in the time-varying model, rather than in the production function. Further investigation and experiments are needed to clarify this issue. The traditional non-steady-state analysis with Steele equation appears to be completely inadequate. Apart from the problem of ill-conditioning of the approach (cf. Cobelli and Mari, 1987), the method, which is based on a monocompartmental approximation, is by far too simplistic in this experimental situation. Physically unrealizable negative values are obtained between 180 and 270 min, and an unphysiological jump is observed at 330 min. The latter phenomenon demonstrate a clear dependence of the result on the time-course of the specific activity, in full accord with the theory developed by Cobelli, Mari and Ferrannini (1987). The results obtained by the steady-state analysis clearly show that a false assumption of steady state may lead to erroneous results. While in fact in the basal period (120-150 min) and in the last period of the clamp (450-480 min) all the calculations are in agreement, they differ markedly between 3CO and 330 min (Table 2). Our analysis thus indicates that the negative values of glucose production reported in recent publications which deal with similar experiments (Bell, Firth and Rizza, 1986; Prager, Wallace and Olefsky, 1986), can be entirely due to a model error. Finally, it is worthwhile to emphasize that the modeling approach here employed, led to a although perfectible, has conclusion which challenge the current physiological thinking, i.e. a rapid and complete inhibitory effect of insulin on glucose production. Even if based on a single case, our results suggest that glucose production is slowly and possibly not completely suppressed by insulin. CONCLUSION Cobelli, Mari and Ferrannini (1987) have already emphasized the importance of proper modeling in the non-steady state. They proposed a two-compartment time-varying model for data analysis which includes the most recent physiological knowledge and overcomes the limitations of previously published models. This paper improves that analysis, which was performed on average data, showing that the model can be identified from
experimental tracer uata in a single individual. Furthermore, glucose production has been calculated as an input estimation problem by taking into proper account the small glucose fluctuations. The results of our analysis indicate that the effect of insulin on glucose production is less rapid and marked than currently thought. The study also indicates, in full agreement with the theory (Cobelli, Mari and Ferrannini, 1987), that oversimplified approaches, such as Steele equation or improper steady-state analysis, are not reliable for computing glucose production. REFERENCES Bell, P. M., R. G. Firth, and R. A. Rizza (1986). Assessment of insulin action in insulin-dependent diabetes mellitus
using
[6
14
3
CJglucose, [3 HJglucose,
3
and [2 HJglucose. Differences in the apparent pattern of insulin resistance depending on the isotope used. ~ Clin. Invest. 78:1479-1486. Bergman, R. N., Y. z. Ider, C. R. (1979). C. Cobelli Bowden, and estimation of insulin Quantitative Physiol. sensitivity. Am. ~ 236 :E66 7-E6 77. R., C. Cobelli, and L. Carson, E. (1983). The Mathematical Finkelstein Modeling of Metabolic and Endocrine Systems. John Wiley, New York. Cobelli, C., G. Toffolo, and E. Ferrannini (1984). A model of glucose kinetics and their control by insulin, compartmental and noncompartmental approaches. Math. Biosci. 72:291-315. Cobelli, C., and A. Mari (1987). Models for the estimation of glucose fluxes in non-steady state from concentration measurements. Preprints 10th IFAC World Congress. Cobelli, C., A. Mari, and E. Ferrannini (1987). Non-steady state: error analysis of Steele's model and developments for glucose kinetics. Am. ~ Physiol._ (in press) . Tobin, and R. DeFronzo, R. A., J. D. Andres (1979). Glucose clamp technique: a method for quantifying insulin secretion ~ Physiol. and resistance. Am. 237 :E214-E223. Ferrannini, E., J. D. Smith, C. Cobelli, G. Toffolo, A. Pilo, and R. A. DeFronzo (1985). Effect of insulin on the distribution and disposition of glucose in man. ~ Clin. Invest. 76:357-364. Insel, P. A., J. E. Liljenquist, J. D. Tobin, R. S. Sherwin, P. Watkins, R. Andres, and M. Berman (1975). Insulin control of glucose metabolism in man. A new kinetic analysis. ~ Clin. Invest. 55:1057-1066. Prager, R., P. Wallace, and J. M. Olefsky (1986) . In vivo kinetics of insulin action on peripheral glucose disposal and hepatic glucose output in normal and obese subjects. ~ Clin. Invest. 78:472-481.
30
C. Cobelli
Radziuk, J., K. H. Norwich, and H. Vranic (1978). Experimental validation of measurements of glucose turnover in nonsteady state. Am . ~ Physiol. 234:E84-E93. (1971). Silverman, L. H. Realization of linear dynamical systems. IEEE Trans. Automat. Contr. 16:554-567 .
Steele, R (1959). Influences of glucose loading and of injected insulin on hepatic glucose output. Ann. NY Acad. ScL 82:420-430.
PI
al.