- Email: [email protected]
A circulatory model for calculating non-steady-state glucose fluxes. Validation and comparison with compartmental models
Computer Methods and Programs in Biomedicine 71 (2003) 269 /281 www.elsevier.com/locate/cmpb
A circulatory model for calculating non-steady-state glucose fluxes. Validation and comparison with compartmental models Andrea Mari a,*, L. Stojanovska b, J. Proietto c, A.W. Thorburn c a
Institute of Systems Science and Biomedical Engineering, National Research Council, LADSEB-CNR, Corso Stati Uniti 4, 35127 Padova, Italy b Faculty of Engineering and Science, Victoria University of Technology, PO Box 14428, Melbourne City 8001, Australia c Department of Medicine, Royal Melbourne Hospital, Parkville, Vic. 3040, Australia Accepted 5 July 2002
Abstract This study presents a circulatory model of glucose kinetics for application to non-steady-state conditions, examines its ability to predict glucose appearance rates from a simulated oral glucose load, and compares its performance with compartmental models. A glucose tracer bolus was injected intravenously in rats to determine parameters of the circulatory and two-compartment models. A simulated oral glucose tolerance test was performed in another group of rats by infusing intravenously labeled glucose at variable rates. A primed continuous intravenous infusion of a second tracer was given to determine glucose clearance. The circulatory model gave the best estimate of glucose appearance, closely followed by the two-compartment model and a modified Steele one-compartment model with a larger total glucose volume. The standard one-compartment model provided the worst estimate. The average relative errors on the rate of glucose appearance were: circulatory, 10%; two-compartment, 13%; modified one-compartment, 11%; standard one-compartment, 16%. Recovery of the infused glucose dose was 939/2, 949/2, 929/2 and 859/2%, respectively. These results show that the circulatory model is an appropriate model for assessing glucose turnover during an oral glucose load. # 2002 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Glucose kinetics; Tracer method; Oral glucose test; Mathematical models
1. Introduction The study of glucose metabolism often requires the calculation of glucose fluxes in non-steadystate conditions, for instance during an oral
* Corresponding author. Tel.: /39-049-829-5753; fax: /39049-829-5763 E-mail address: [email protected] (A. Mari).
glucose load or a glucose clamp. This calculation is based on tracer methodology, and requires a mathematical model. Unless the specific activity of the tracer is kept constant, which is often difficult or impossible, the accuracy of the calculated glucose fluxes depends on the model used. The most commonly used models are compartmental models with one or two compartments [1,2]. These models have proved useful in many
0169-2607/02/$ - see front matter # 2002 Elsevier Science Ireland Ltd. All rights reserved. PII: S 0 1 6 9 - 2 6 0 7 ( 0 2 ) 0 0 0 9 7 - 4
270
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
situations. However, discrepancies between the glucose fluxes calculated with these models and with independent methods have been found. Steele’s one-compartment model [1] has been shown to underestimate glucose production during a standard euglycemic hyperinsulinemic glucose clamp [3]. In experiments which simulate an oral load with an exogenous glucose infusion, both Steele’s model and the two-compartment model by Radziuk and colleagues [2] predicted the rate of appearance with some distortion [2,4]. During an oral glucose load, the rate of appearance calculated using the two-compartment model and the arteriovenous difference method were also somewhat different [5]. The problems encountered with the compartmental models may not originate entirely from the inadequacy of the models (cf. [5]). Nevertheless, compartmental models do not represent the physiological system as it is [6,7], and for this reason may introduce errors [8]. To avoid this drawback, circulatory models have been developed [8,9], which more appropriately represent the physiological system and are based on the solid principles of the theory developed for organ kinetics by Zierler [10,11]. Circulatory models have been used in various artificial non-steady-state conditions (e.g. [12,13]), but a model for general use has not been proposed and tested. In this work, we present a circulatory model for the calculation of nonsteady-state glucose fluxes in the general case. We evaluate the model performance in experiments in anesthetized rats that simulate an oral glucose load through an exogenous glucose infusion, by comparing the model-calculated rate of appearance with the known glucose infusion. We also compare the circulatory model performance with that of the more traditional compartmental approaches.
2. Methods 2.1. Experimental procedures 2.1.1. Basal glucose kinetics To determine the parameters for basal glucose kinetics, six adult male Sprague /Dawley rats were anesthetized with pentobarbitone sodium (Nem-
butal, Boehringer Ingelheim, NSW, Australia) administered by intraperitoneal injection (60 mg kg1 body weight). A catheter was placed in the right jugular vein (for tracer infusion) and left carotid artery (for blood sampling) [14]. The rats were given a single injection of [6-3H]glucose (10 mCi) and blood samples were taken at 0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 25, 30, 40, 50 and 60 min for analysis of plasma tracer. 2.1.2. Simulated oral load Another group of five male Sprague /Dawley rats were anesthetized and implanted with catheters as for the basal experiments. At time 0 min, a primed (185 mmol kg1) continuous infusion (3.8 mmol min 1 kg1) of [6,6-D2]glucose (Tracer Technologies, Somerville, MA, USA) in 0.9% NaCl was commenced. Three blood samples were taken at 60, 65 and 70 min for the estimation of the tracer-to-tracee ratio of [6,6-D2]glucose and plasma glucose. At 70 min, an infusion of 15 mCi ml1 [6-3H]glucose in 25% glucose commenced in an algorithm designed to mimic the appearance of gut-derived glucose following the administration of an oral glucose load. From time 70 to 90 min the glucose infusion rate was increased every minute to a maximum infusion rate of /48 mmol min1. This infusion rate was maintained for 15 min (time 90/105 min) after which the infusion rate was decreased every minute until it was stopped at 125 min. Blood samples (300 ml) were taken at 5 min intervals between 70 and 125 min for measurement of plasma glucose, the tracer-to-tracee ratio of [6,6-D2]glucose, and [6-3H]glucose specific activity. At the end of the study, timed collections of the [6,6-D2]glucose infusate were taken for accurate measurement of the constant [6,6-D2]glucose infusion rate. Timed collections of the variable glucose infusate were also taken to determine accurate time courses of the cold glucose delivery rate in each experiment. The [6-3H]glucose specific activity of the variable glucose infusate was also determined. 2.2. Analytical methods Plasma and infusate glucose levels, and plasma [6-3H]glucose tracer radioactivity and [6,6-D2]glu-
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
cose enrichment were measured as described previously [14,15]. [6,6-D2]glucose enrichment was analyzed using a gas chromatograph mass spectrometer (Shimadzu model GCSM-QP2000, Shimadzu Corporation, Kyoto, Japan) on ion monitoring mode to determine the relative intensity of the 98 and 100 molecular weight fragments i.e. (peak area at 100)/(peak area at 100/peak area at 98). A standard curve was run in parallel to convert relative intensity into relative tracer abundance (A ). The tracer-to-tracee ratio was then calculated as A /(1/A ). The study protocol was approved by the Royal Melbourne Hospital’s Animal Ethics Committee. 2.3. Circulatory model 2.3.1. Model structure Glucose kinetics are described using the circulatory model of Fig. 1, the mathematical theory of which has been presented previously [8,9,16]. In the model, the heart chambers and the lungs are represented in the heart /lungs block (upper block), while the remaining tissues are lumped into the periphery block (lower block). Blood flow for both blocks is cardiac output. In the experimental configuration adopted here, glucose is injected at the venous side, while blood is sampled at the arterial side (as indicated in Fig. 1). The two tissue blocks shown in Fig. 1 can be regarded as single inlet /single outlet organs, coupled in a feedback arrangement. The basis of the mathematical description of a block (as for an organ) is its impulse response [17]. The impulse response is the tracer efflux (concentration times blood flow) observed at the organ outlet after a bolus injection of a unit tracer dose at the organ inlet (this definition assumes tracer does not
recirculate). Cardiac output and the impulse responses of the blocks fully determine the circulatory model, and are specified in this section. Cardiac output was expressed as flow of blood per kilogram of body weight. Cardiac output was fixed to the value of 236 ml min 1 kg1, which is appropriate for anesthetized rats [18]. Since glucose concentration was measured in plasma and not in blood, this cardiac output value was corrected using the ratio between glucose concentration in whole blood and in plasma (0.53 in our rats) to ensure that the product of concentration and flow gives the actual glucose mass flux. A constant cardiac output value during the test was assumed (see Section 4). The impulse response of the heart /lungs block [rlung(t)] was modeled as a two-exponential function, which is the minimal representation ensuring a time course of the response in agreement with the experimental data. The impulse response starts from 0, rapidly increases to a peak value, and returns to zero as a single-exponential function. To simplify the notation, it is convenient to represent rlung(t ) as the convolution of two single-exponential functions: rlung (t)bebt vevt
(1)
where the symbol / denotes the convolution operator. As the coefficients of the single-exponential functions equal their exponents, the integral from 0 to of each exponential is 1. This also implies that the integral from 0 to of rlung(t) is 1, i.e. the glucose fractional extraction in the heart / lungs system is 0 [12]. The exponent that characterizes the rising phase of rlung(t) was fixed (b / 15 min 1, cf. [12]). The exponent of the decaying phase (v , min1) was calculated using the property that the mean transit time of rlung(t) is the ratio of volume to blood flow [17]. From this property, the following equation was obtained (cf. [12]): v
Fig. 1. The circulatory model.
271
bCO bVlung CO
(2)
where cardiac output (CO) and the heart /lungs glucose volume (Vlung) were assumed known
272
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
(CO /0.53 /236 ml min 1 kg1, see above; Vlung /17 ml kg1, see [12]). The impulse response of the periphery block [rper(t)] in the basal steady-state period can be expressed as [17]: rper (t) (1E)p(t)
(3)
where E (dimensionless) is the glucose fractional extraction and p (t) is the glucose transit time density function. The transit time density function p (t ) was modeled as a four-exponential function, starting from zero, rapidly increasing to a peak value, and returning to 0 as a three-exponential function. Similarly to rlung(t), we have represented p(t) as a convolution of a single-exponential function with a three-exponential function: p(t) gegt [w1 a1 ea1 t w2 a2 ea2 t (1w1 w2 )a3 ea3 t ]
(4)
where g , a1, a2, a3 (min 1) and w1, w2 (dimensionless) are parameters. In Eq. (4), the integral from 0 to of the threeexponential function in square brackets is w1/ w2/(1/w1/w2), i.e. 1. This ensures that the integral from 0 to of p (t) is also 1, as required for a transit time density function. The parameters w1 and w2 determine the contribution in the total response of the three exponential terms. As b in rlung(t), g determines the initial rising phase of p (t ). The parameter g was fixed (g/10 min 1, [12]), while a1, a2, a3, w1, w2 and E were estimated from the basal kinetic experiments. Eq. (3) can be extended to the non-steady-state [9]. In non-steady-state, the glucose fractional extraction E varies with time under the action of insulin, and the transit time density function p(t) may also change. However, in the present model we have assumed that the transit time density function is not affected by insulin. This is supported by the finding that the main effect of insulin is exerted on the fractional extraction of the periphery, while the transit time density function does not change substantially [8,12]. Thus, Eq. (3) was assumed to be valid also in the non-steady-
state period, with the same p (t) obtained from the basal period, and E (t) time-varying. To simulate the model of Fig. 1, the multiexponential equations for the impulse responses of the heart /lungs and periphery blocks are first represented as ordinary differential equations. These differential equations are then combined according to the feedback arrangement of Fig. 1, and then solved numerically using MATLAB (see Appendix A for details). 2.3.2. Analysis of basal steady-state experiments For the present circulatory model global identifiability of the model parameters is ensured in steady-state [16]. Mean parameters of the periphery transit time density function p (t) and the basal glucose fractional extraction E (six parameters) were obtained by least-squares fit of the mean [3-3H]glucose concentration in the basal kinetic experiments. The estimated parameters were subsequently used in the analysis of the non-steadystate data (a single parameter set for all nonsteady-state experiments). Parameter estimation was performed using a Levenberg /Marquardt least-squares algorithm (MATLAB function leastsq). 2.3.3. Analysis of non-steady-state experiments In non-steady-state, as p(t) was assumed to be time-invariant, the only time-varying parameter to be determined was E (t). E (t ) was assumed constant in the basal period preceding the start of glucose infusion, and represented as a piecewise constant function of time over 2 min time intervals thereafter (28 elements, E1 /E28, were necessary to cover the 55 min non-steady-state period). The Ek values were estimated from [6,6-2H2]glucose concentration using the model and a method traditionally used for deconvolution, as previously described [12]. In brief, the Ek values were estimated by least-squares fit of [6,6-2H2]glucose concentration, including in the least-squares a penalty term dependent on the second derivative of E (t ), calculated by forward differences of the Ek values. This additional term is necessary to obtain a smooth E (t), which would otherwise exhibit spurious oscillations.
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
273
2.3.4. Exogenous glucose appearance By dividing the [3-3H]glucose concentration by the specific activity of the glucose infusate, the glucose concentration component due to the exogenous glucose infusion was calculated (‘exogenous’ glucose concentration). From the exogenous glucose concentration and the model (using the tracer-determined E (t)), the rate of appearance of exogenous glucose was calculated using a deconvolution method [12]. For this purpose, we have approximated the exogenous glucose appearance rate as a piecewise constant function of time over 2 min time intervals, as E (t). 2.3.5. Glucose production The glucose concentration component due to glucose production (‘endogenous’ glucose concentration) was calculated by subtracting the exogenous glucose concentration from the total (measured) glucose concentration. From the endogenous glucose concentration and the model, glucose production was calculated using a deconvolution method, as for exogenous glucose. Glucose production was approximated as a constant value before t/70 min, and a piecewise constant function of time over 2 min time intervals thereafter. The calculated values were corrected for the non-negligible mass flux of the stable-isotope tracer (/3.8 mmol min 1 kg1), as the endogenous glucose concentration component includes the contribution of the tracer.
Fig. 2. Compartmental models. Top: single-compartment model; bottom: two-compartment model. See Eqs. (5) and (6) for the symbols.
We have used two values for the glucose distribution volume VS. The first value represents a standard volume choice (pool fraction of 0.65 with a total glucose volume of 250 ml kg1, i.e. VS /0.65 /250 ml kg1, [19]). The second value is the total glucose volume as calculated with the circulatory model in the basal period (VS /191 ml kg1). The time course of tracer and glucose concentration (exogenous or endogenous), and the derivative of specific activity were obtained from the data fit of the circulatory model (see Appendix A). This provides the necessary data smoothing and ensures that the differences between Steele’s model and the circulatory model are not due to differences in data smoothing and interpolation.
2.4. Compartmental models 2.4.1. One-compartment model The glucose rate of appearance [Ra(t )] was calculated every 2 min (time instants t0 /0, t1 / 2,. . . min) using Steele’s single-compartment model (Fig. 2 top, [1]), according to the equation (see [19]): Ra (tk )
R(t a k) a(tk )
VS
G(tk )a(t ˙ k) a(tk )
(5)
where Ra (tk ) is the tracer infusion rate, G (tk ) is arterial glucose concentration, a (tk ) is the glucose specific activity (the dot represents the time derivative), and VS is the compartment volume.
2.4.2. Two-compartment model The rate of glucose appearance was calculated using the two-compartment model of Fig. 2, according to the equations developed in [19]: Ra (tk )
Ra(tk )
G(tk )a(t ˙ k) k k V1 V1 12 21 a(tx ) a(tk ) k12 g(tk ) g(tk ) a(tk )
g(tk1 )b1 g(tk )b2 G(tk )b3 G(tk1 ); g(t0 )G(t0 )
274
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
g(tk1 ) b1 g(tk )b2 G(tk )b3 G(tk1 ); g(t0 )G(t0 )
(6)
where the star denotes the tracer, V1 is the volume of the first compartment, k21 and k12 are the intercompartmental rate constants (cf. Fig. 2), g(t) and g *(t ) are auxiliary variables (delayed glucose concentration), and the bk ’s (cf. [19]) are constants calculated from k12 and the calculation interval tk1/tk , which was 2 min as for Steele’s and the circulatory model. The parameters of the two-compartment model were obtained from the basal steady-state experiments, as for the circulatory model. The twocompartment model did not fit the initial part of the tracer disappearance curves accurately, as expected [20]. We have thus discarded the first 2 min in the analysis of the basal tracer curves. As for Steele’s model, concentrations of tracer and glucose, and the derivative of specific activity were obtained from the data fit of the circulatory model.
2.5. Evaluation of the model accuracy The accuracy of the models was evaluated by comparing the model-calculated rate of appearance with the actual glucose infusion rate, the error being the difference of the two. Besides the direct comparison of the time course of the mean rate of appearance, we have calculated the following indices of model performance: (1) percent recovery of the infused glucose dose, i.e. the ratio between the integral of the model-calculated rate of appearance and the total quantity of infused glucose, expressed in percent; (2) integral absolute error, i.e. the integral of the absolute value of the error on the rate of appearance; (3) cumulative distribution of the percent error (absolute value of the error divided by the actual glucose infusion rate). The cumulative distribution function was calculated by pooling the error values in all rats and at all time points. The frequency of occurrence was divided into ten percentiles.
3. Results 3.1. Basal glucose kinetics 3.1.1. Circulatory model The parameters of the circulatory model estimated from the tracer bolus injection in the basal state are reported in Table 1. The coefficients of variation of the parameters provided by the leastsquares algorithm were below 23%. The circulatory model predicted the [3-3H]glucose concentration curve accurately. 3.1.2. Two-compartment model The two-compartment model parameters were: V1 /119 ml kg 1; k01 /0.033; k21 /0.073; k12 / 0.11 min1; clearance /3.9 ml min 1 kg1; total volume /197 ml kg1. The coefficients of variation were below 20%. In the first 2 min, which were excluded from the fit, the two-compartment model underestimated [3-3H]glucose concentration, while in the remaining period the fit was accurate. 3.2. Non-steady-state glucose kinetics 3.2.1. Glucose and tracer concentrations Fig. 3 shows the mean [6,6-2H2]glucose concentration, total glucose concentration, and its exogenous and endogenous components. The solid lines represent the circulatory model fit. 3.2.2. Rate of appearance of exogenous glucose The mean rate of appearance of exogenous glucose calculated with the circulatory and the two-compartment model are shown in Fig. 4. The results for Steele’s model, obtained using the two glucose distribution volumes, are shown in Fig. 5. The standard Steele calculation (Fig. 5 top) was the least accurate, while the performance of the other models was not markedly different on average. All the models underestimated the downslope of the glucose rate of infusion. The twocompartment model most accurately predicted the downslope, while the circulatory model most accurately predicted the upslope. The performance of the modified Steele’s model with the volume set to the total glucose volume was similar to that of the circulatory model.
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
275
Table 1 Model parameters in the basal state for the anesthetized rats CO ml min 1 kg 1
125
E
0.0293
Heart /lungs
Periphery
b (min1)
v (min 1)
g (min1)
a1 (min1)
a2 (min 1)
a3 (min 1)
w1
w2
15.0
14.4
10.0
3.69
0.534
0.0788
0.684
0.265
CO, cardiac output; E , fractional extraction. For the remaining symbols see Eqs. (1) /(4). In the basal state, the clearance was 3.7 ml min 1 kg1 and the total volume 191 ml kg1. The parameter v was calculated from cardiac output and the volume of the heart / lungs block (Vlung /17 ml kg 1) using Eq. (2). The cardiac output value reported in the table is corrected for the distribution of glucose between plasma and red cells (CO /0.53 /236 /12 ml min 1 kg1, see Section 2).
Fig. 3. Mean (9/S.E.M.) tracer concentration, total glucose concentration, exogenous glucose concentration and endogenous glucose concentration. The solid line represents the circulatory model fit.
Table 2 reports the recovery of the infused glucose and the average absolute integral error for all models. Depending on the model, the correlation coefficient between the integral of the calculated rate of appearance and the total quantity of infused glucose was /0.85 /0.88, with a borderline significance (P :/0.05 /0.07). These integral accuracy indexes basically confirmed the results of Figs. 4 and 5.
Fig. 4. Mean (9/S.E.M.) exogenous glucose rate of appearance calculated with the circulatory (top) and the two-compartment model (bottom). The solid lines represent the true mean glucose infusion rate (standard errors are omitted for clarity).
The distribution of the pooled relative error is shown in Fig. 6. The curves show the percentage of the rate of appearance values that are affected by an error which does not exceed a specified value. For instance, for the circulatory model 60% of the calculated rates of appearance (read off the ordinate) are affected by an error which is not greater than 12% (read off the abscissa). In this representation, the curves of the more accurate models lie above those of the less accurate models. According to this analysis, which shows the error
276
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
Fig. 6. Cumulative distribution of the pooled relative error. Results for all models are reported, as indicated in the figure. The ordinate shows the percentage of the calculated rates of appearance that do not have an error which exceeds the level shown on the abscissa. For clarity, only 90% of the calculated rates are included. The error value on the abscissa which corresponds to 100% of the calculated rates (maximal error) is in fact much higher: circulatory model, 145%; two-compartment model, 121%; Steele’s model, 107% (V/162.5 ml kg 1) and 134% (V /191 ml kg 1). Fig. 5. Mean (9/S.E.M.) exogenous glucose rate of appearance calculated with Steele’s model and the two glucose volume values indicated in the figure. The solid lines represent the true mean glucose infusion rate (standard errors are omitted for clarity).
in 90% of the calculated points for clarity, the overall performance of the circulatory model is the best and that of the standard Steele’s model the worst. However, with the exception of the standard Steele’s model, the differences are small. 3.2.3. Glucose production Fig. 7 shows mean glucose production calculated using the circulatory model. Glucose production decreased to 0 at 25 min after the start of glucose infusion, and slowly returned towards the
Fig. 7. Mean (9/S.E.M.) glucose production calculated with the circulatory model.
Table 2 Recovery and integral absolute error
Recovery (%) Integral absolute error (mmol kg 1) a
Circulatory model
Two-compartment model
Steele’s model V/162.5 ml kg 1
Steele’s model V/191 ml kg 1
939/2 7099/89
949/2 7179/72
859/2 8629/106a
929/2 7719/93
P B/0.05 or less vs. circulatory model (paired t -test). The recovery was significantly less than 100% for all models (P B/0.05 or less, one-sample t -test).
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
basal value thereafter. Glucose production as calculated using the compartmental models was similar to that reported in Fig. 7. The largest difference was observed with the standard Steele’s model, which yielded a higher glucose production on the average, but the maximum difference was less that 4 mmol min 1 kg1.
4. Discussion In our experiments, the circulatory model predicted the actual glucose infusion rate with an average error of about 10%. The rising phase of the glucose infusion was precisely matched, while the falling phase was underestimated. The biggest error was observed in this period, but on average the error was about 15%, and in 90% of all the calculated rates of appearance it did not exceed 30%. According to the distribution of the error on the rate of appearance (Fig. 6, Table 2), the performance of the circulatory model was the best. However, none of the models matched the glucose infusion rate accurately, and the evaluation of the relative performance is thus somewhat subjective. Indeed, the accuracy of the two-compartment model and Steele’s model with the larger volume was only slightly inferior to that of the circulatory model. As concerns the circulatory model, the error on the glucose rate of appearance may originate from two major sources. The first is the use of mean instead of individual parameters. Non-steady-state calculations are made using a fixed parameter set for all rats, while it is expected that the individual model parameters are somewhat different from rat to rat. This may introduce an error in the individual rates of appearance. We could not evaluate the magnitude of this error because experimental limitations prevented us from assessing the basal parameters in each rat used for the non-steady-state experiments. However, it is unlikely that this error is also significant for the mean rate of appearance, as the interindividual parameter differences are random and tend to cancel out in the mean. The second error concerns limited information on the system. The circulatory model is a physical representation of the system, but
277
assumes a fixed cardiac output and a timeinvariant distribution of the transit times p(t) of the periphery because the little is known about the changes that the cardiac output and p (t ) undergo in this and other experimental conditions. During glucose challenges or hyperinsulinemia in conscious human subjects, changes in cardiac output do not exceed 20% [21 /23]. Changes from basal in p(t) have been determined during hyperinsulinemic euglycemic glucose clamps, directly in dog [12] and indirectly in man [8], and in both cases are not substantial. In hyperinsulinemic hyperglycemic glucose clamps in anesthetized rats, we have also found small changes in p (t ) (/20% difference in mean transit time, results not shown). These findings justify the assumption of this study of constant cardiac output and time-invariant p (t). Additional experimental and theoretical work are, however, needed to resolve this issue, and to clarify whether more precise assumptions result in a better model performance, or if factors other than the model are the cause of the observed discrepancies, as pointed out previously (e.g. [5]). Unfortunately, the determination of the time course of cardiac output and the changes of p (t ) is a difficult problem to solve, and this limits the possibilities of model improvement. In regard to cardiac output, we have assumed the most appropriate value for anesthetized rats [18]. In principle, cardiac output cannot be ignored, because it does have a role in glucose kinetics and thus in the calculation of glucose fluxes. However, the calculations are not very sensitive to the cardiac output value. Only with rapid changes of fluxes, concentrations and specific activity, which are not observed in the present experiments, is it expected that the calculated glucose fluxes are influenced by the cardiac output value. It is in fact known that in steady-state, or with negligible changes in specific activity, cardiac output is not relevant for calculating glucose turnover. The compartmental models considered here performed differently. The two-compartment model performance was only slightly inferior to that of the circulatory model. The standard Steele equation markedly distorted the shape of the glucose infusion rate, but the poor performance
278
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
parameters of Tables 1 and 3 to experimental conditions in which the cardiac output value is different from that reported here. In conclusion, we have presented a circulatory model for the calculation of the rate of appearance in non-steady-state conditions and we have assessed its accuracy experimentally. The circulatory model gives a physiological representation for glucose kinetics, and its accuracy is to some extent superior to that of the more standard Steele and two-compartment models.
of this model could be corrected by using a larger glucose volume. This stratagem has theoretical support when the specific activity changes slowly [24]. However, the use of larger glucose volume does not ensure sufficient accuracy in general. For instance, the calculation of glucose production during a standard glucose clamp is not accurate with Steele’s model, regardless of the glucose volume used [24]. Therefore, in this experimental condition compartmental models are sufficiently accurate, but, as demonstrated for the Steele’s model, it is not necessarily true that the error due to the incorrect representation of the physiological system is small in all situations. The circulatory model presented here is suitable for application to general non-steady-state conditions, as is the Steele’s one-compartment model and the two-compartment model. To use the circulatory model, the cardiac output value and the parameters of p (t) are required. It is clearly important that the parameters used are appropriate for the specific experimental condition. If possible, the model parameters should be estimated in each individual experiment, using for instance the tracer equilibration period that normally precedes the actual non-steady-state phase. In this work, however, a single set of parameters (Table 1) was used for all rats because individual estimates could not be obtained. These values are appropriate for anesthetized rats. As it is useful to provide analogous values for humans, we report in Table 3 a set of parameters derived from a study in lean and obese subjects [13]. We also present in Appendix B some results useful to extend the
Acknowledgements This study was supported in part by a grant from the Italian National Research Council.
Appendix A In this appendix, the differential equations of the model are derived. Additional details on the equations can be found in [12], which presents a similar model using the same notation. Basic differential equation The key for the transformation of the impulse responses (Eqs. (1), (3) and (4)) into differential equations is the representation of the convolution by means of a differential equation. If the impulse response of the heart /lungs block were the single-
Table 3 Model parameters in the basal state for humans (from [13]) CO (ml min1 kg 1)
E
Heart /lungs b (min
66
0.0223
15.0
1
)
Periphery v (min 5.16
1
)
g (min 1)
a1 (min1)
a2 (min 1)
a3 (min1)
w1
10.0
3.05
0.382
0.0388
0.430
w2 0.505 1
CO, cardiac output; E , fractional extraction. For the remaining symbols see Eqs. (1) /(4). The total volume was 223 ml kg . The parameter v was calculated from cardiac output and the volume of the heart /lungs block (Vlung /17 ml kg1) using Eq. (2). The cardiac output value reported in the table is corrected for the distribution of glucose between plasma and red cells (see Section 2). In [13], cardiac output (flow of blood) was assumed to be 78 ml min 1 kg1, and the correction factor 0.84 (CO /0.84/78/66 ml min 1 kg1). In [13], cardiac output was expressed in ml min 1 m 2 (cardiac index), and is converted here into ml min 1 kg 1 assuming a body surface area to body weight ratio of 1.7/70.
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
exponential function k elt , the relationship between the input (glucose concentration in the right atrium, Gra(t)) and the output (arterial glucose concentration, Ga(t)) would be the convolution Ga (t)kelt Gra (t)
x˙ 3 (t)a1 x3 (t)w1 (1E)Ga (t) x˙ 4 (t)a2 x4 (t)w2 (1E)Ga (t) x˙ 5 (t)a3 x5 (t)(1w1 w2 )(1E)Ga (t)
(A1)
This convolution can be represented by means of the differential equation x(t)lx(t)G ˙ ra (t) Ga (t)kx(t)
(A2)
where x (t) is an additional variable (state variable), and the dot represents the time derivative. Eq. (A2) is a state space representation of a linear system, in a canonical form. Its state variables do not represent physical quantities. The initial condition for Eq. (A2) (and the successive equations) depends on whether endogenous glucose, oral glucose or the tracer is considered. For oral glucose and the tracer, which are not present in the system before time 0, x (0) /0. For endogenous glucose, x (0) is the steady-state value, i.e. the solution of the algebraic equations obtained from Eq. (A2) by letting x(t) ˙ 0:/
x˙ 6 (t)gx6 (t)a1 x3 (t)a2 x4 (t)a3 x5 (t) Gmv (t)gx6 (t)
Because the two-exponential impulse response of the heart /lungs block (Eq. (1)) is represented as a convolution of two single-exponential functions, the relationship between Gra(t ) and Ga(t) can be expressed as a cascade of two differential equations like Eq. (A2), i.e.
(A4)
where Gmv(t) is mixed-venous glucose concentration (output of the periphery block), and x3 /x6 are state variables. In Eq. (A4), the three exponential terms in the square brackets in Eq. (4) are represented by a sum of three differential equations, while the last two equations express the convolution with the exponential function geg . t
Whole-body differential equations In the feedback arrangement of Fig. 1, the input of the heart /lungs block is the output of the periphery block [Gmv(t)] plus a term representing the appearance of glucose, which is Ra(t)/CO. Therefore, the whole set of differential equations describing the circulatory model is: x˙ 1 (t)vx1 (t)gx6 (t)
Heart /lung differential equations
Ra (t) CO
x˙ 2 (t)bx2 (t)vx1 (t) x˙ 3 (t)a1 x3 (t)w1 (1E)bx2 (t) x˙ 4 (t)a2 x4 (t)w2 (1E)bx2 (t) x˙ 5 (t)a3 x5 (t)(1w2 w2 )(1E)bx2 (t)
x˙ 1 vx1 (t)Gra (t)
x˙ 6 (t)gx6 (t)a1 x3 (t)a2 x4 (t)a3 x5 (t)
x˙ 2 bx2 (t)vx1 (t)
Ga (t) bx2 (t)
Ga (t)bx2 (t)
279
(A5)
(A3)
where x1 and x2 are state variables. Periphery differential equations Similarly, the impulse response of the periphery block (Eqs. (3) and (4)) is represented by the four differential equations:
Numerical solution Eq. (A5) is a standard representation in state space form of a linear system of differential equations. The solution of the system yields the glucose concentration (tracer, exogenous or endogenous) that corresponds to a given rate of
280
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
appearance. Eq. (A5) is valid both in steady and non-steady-state. In steady-state E is constant, and Eq. (A5) represents a linear time-invariant system. In non-steady-state E is time-varying, and the system is time-varying. However, as E(t) is assumed to be piecewise constant, the system is in reality piecewise time-invariant, which facilitates the calculation of the solution. To solve the differential equations numerically, Eq. (A5) was expressed in matrix form as: x(t)Ax(t)BR ˙ a (t) Ga (t)Cx(t)
(A6)
where A, B and C are the standard system matrices, x (t) is the state vector in which the state variables xk (t) are stacked, and the initial value x (0) is as discussed for Eq. (A2). Because E (t) is assumed to be piecewise constant, the matrix A is also piecewise constant. Eq. (A6) was integrated by updating the state vector step by step over the time grid resulting from the union of the data time points and the time grid for E (t) and Ra(t). In this way, in each step A, B and C are constant matrices, and x (tk1) was updated from x (tk) after calculation of the matrix exponential eA(tk1tk ) (MATLAB function expm), using the standard equations for state update with piecewise constant inputs [25]. Except for Ra(t ) and E (t), all the parameters of Eq. (A5) are reported in Table 1 (Table 3 for humans). /
Appendix B This appendix briefly addresses the problem of adapting the parameters of Tables 1 and 3, which are appropriate for anesthetized rats or lean and obese humans, to experimental conditions in which the cardiac output value is different from that reported here. The glucose distribution volume in the periphery block is the product of cardiac output, 1-fractional extraction and the mean transit time of the periphery [8]. The glucose distribution volume is likely to be substantially independent from the cardiac output value. If the circulatory model is used in conditions in which cardiac output is different (e.g. conscious vs. anesthetized rats), it is necessary to change the periphery mean transit time in order that the glucose volume remains constant. This is achieved by adjusting the exponents ak of the periphery impulse response as follows. The glucose volume is calculated from the periphery impulse response (Eq. (4)) as the sum of volume terms Vk , one for each exponential (the volume contribution due to g is negligible): Vk CO(1E)
(B1)
To maintain the total glucose volume constant, it is logical to increase in each volume term the exponents ak in proportion to the cardiac output value. Thus, in an experimental condition in which cardiac output is CO?, the new exponents, ak? ?, will be:
Specific activity derivative
ak? The numerical solution of Eq. (A5) provides a smooth time course of Ga(t) (glucose and tracer concentration) and also of its derivative, G˙ a (t): G˙ a (t) is in fact proportional to x˙ 2 (t) (Eq. (A5), last equation), and is thus a linear combination of x1(t) and x2(t) (Eq. (A5), second equation), which are smooth. From the glucose and tracer smooth Ga(t) and G˙ a (t); a smooth time course of the specific activity derivative can be obtained. The smoothed variables were used with the compartmental models (Eqs. (5) and (6)).
wk ak
CO? ak CO
(B2)
where CO and ak are the values reported in Table 1 or Table 3. Eq. (B2) is an approximate correction, which assumes that cardiac output does not affect substantially the distribution of glucose. If there are reasons to think that cardiac output affects glucose distribution, the model parameters should be determined experimentally. The cardiac output problem is not apparent using compartmental models. However, this does
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
not imply that the compartmental model parameters are unaffected by cardiac output. It is in fact expected that in conditions in which cardiac output is different, also the compartmental model parameters (e.g. k12 and k21) differ [8].
References [1] R. Steele, Influences of glucose loading and of injected insulin on hepatic glucose output, Annals of the New York Academy of Sciences 82 (1959) 420 /430. [2] J. Radziuk, K.H. Norwich, M. Vranic, Experimental validation of measurements of glucose turnover in nonsteady-state, American Journal of Physiology 234 (1978) E84 /E93. [3] D.T. Finegood, R.N. Bergman, M. Vranic, Estimation of endogenous glucose production during hyperinsulinemiceuglycemic glucose clamps. Comparison of unlabeled and labeled exogenous glucose infusates, Diabetes 36 (1987) 914 /924. [4] J. Proietto, F. Rohner-Jeanrenaud, E. Ionescu, J. Terrettaz, J.F. Sauter, B. Jeanrenaud, Non-steady-state measurement of glucose turnover in rats by using a onecompartment model, American Journal of Physiology 252 (1987) E77 /E84. [5] A. Mari, J. Wahren, R.A. DeFronzo, E. Ferrannini, Glucose absorption and production following oral glucose: comparison of compartmental and arterio-venous difference methods, Metabolism 43 (1994) 1419 /1425. [6] K. Zierler, A critique of compartmental analysis, Annual Review of Biophysics and Bioengineering 10 (1981) 531 / 562. [7] K. Zierler, Whole body glucose metabolism, American Journal of Physiology 276 (1999) E409 /E426. [8] A. Mari, Circulatory models of intact-body kinetics and their relationship with compartmental and noncompartmental analysis, Journal of Theoretical Biology 160 (1993) 509 /531. [9] A. Mari, Calculation of organ and whole-body uptake and production with the impulse response approach, Journal of Theoretical Biology 174 (1995) 341 /353. [10] P. Meier, K.L. Zierler, On the theory of the indicatordilution method for measurement of blood flow and volume, Journal of Applied Physiology 6 (1954) 731 /744. [11] K.L. Zierler, Theory of the use of arteriovenous concentration differences for measuring metabolism in steady and non-steady-states, Journal of Clinical Investigation 40 (1961) 2111 /2125.
281
[12] O.P. McGuinness, A. Mari, Assessment of insulin action on glucose uptake and production during a euglycemichyperinsulinemic clamp in dog: a new kinetic analysis, Metabolism 46 (1997) 1116 /1127. [13] A. Natali, A. Gastaldelli, S. Camastra, A.M. Sironi, E. Toschi, A. Masoni, E. Ferrannini, A. Mari, Dose-response characteristics of insulin action on glucose metabolism: a non-steady-state approach, American Journal of Physiology 278 (2000) E794 /E801. [14] C.J. Nolan, J. Proietto, The feto-placental glucose steal phenomenon is a major cause of maternal metabolic adaptation during late pregnancy in the rat, Diabetologia 37 (1994) 976 /984. [15] A. Thorburn, J. Muir, J. Proietto, Carbohydrate fermentation decreases hepatic glucose output in healthy subjects, Metabolism 42 (1993) 780 /785. [16] A. Mari, Determination of the single-pass impulse response of the body tissues with circulatory models, IEEE Transactions on Biomedical Engineering 42 (1995) 304 / 312. [17] N.A. Lassen, W. Perl, Tracer Kinetic Methods in Medical Physiology, Raven Press, New York, 1979. [18] M.W. Crawford, J. Lerman, V. Saldivia, F.J. Carmichael, Hemodynamic and organ blood flow responses to halothane and sevoflurane anesthesia during spontaneous ventilation, Anesthesia and Analgesia 75 (1992) 1000 / 1006. [19] A. Mari, Estimation of the rate of appearance in the nonsteady-state with a two-compartment model, American Journal of Physiology 263 (1992) E400 /E415. [20] C. Cobelli, G. Toffolo, E. Ferrannini, A model of glucose kinetics and their control by insulin, compartmental and noncompartmental approaches, Mathematical Biosciences 72 (1984) 291 /315. [21] S. Jern, Effects of acute carbohydrate administration on central and peripheral hemodynamic responses to mental stress, Hypertension 18 (1991) 790 /797. [22] A.D. Baron, Hemodynamic actions of insulin, American Journal of Physiology 267 (1994) E187 /E202. [23] K. Masuo, H. Mikami, T. Ogihara, M.L. Tuck, Mechanisms mediating postprandial blood pressure reduction in young and elderly subjects, American Journal of Hypertension 9 (1996) 536 /544. [24] C. Cobelli, A. Mari, E. Ferrannini, Non-steady-state: error analysis of Steele’s model and developments for glucose kinetics, American Journal of Physiology 252 (1987) E679 /E689. [25] C.T. Chen, Linear system theory and design, third ed. (Chapter 4), Oxford University Press, New York, 1999.
A circulatory model for calculating non-steady-state glucose fluxes. Validation and comparison with compartmental models Andrea Mari a,*, L. Stojanovska b, J. Proietto c, A.W. Thorburn c a
Institute of Systems Science and Biomedical Engineering, National Research Council, LADSEB-CNR, Corso Stati Uniti 4, 35127 Padova, Italy b Faculty of Engineering and Science, Victoria University of Technology, PO Box 14428, Melbourne City 8001, Australia c Department of Medicine, Royal Melbourne Hospital, Parkville, Vic. 3040, Australia Accepted 5 July 2002
Abstract This study presents a circulatory model of glucose kinetics for application to non-steady-state conditions, examines its ability to predict glucose appearance rates from a simulated oral glucose load, and compares its performance with compartmental models. A glucose tracer bolus was injected intravenously in rats to determine parameters of the circulatory and two-compartment models. A simulated oral glucose tolerance test was performed in another group of rats by infusing intravenously labeled glucose at variable rates. A primed continuous intravenous infusion of a second tracer was given to determine glucose clearance. The circulatory model gave the best estimate of glucose appearance, closely followed by the two-compartment model and a modified Steele one-compartment model with a larger total glucose volume. The standard one-compartment model provided the worst estimate. The average relative errors on the rate of glucose appearance were: circulatory, 10%; two-compartment, 13%; modified one-compartment, 11%; standard one-compartment, 16%. Recovery of the infused glucose dose was 939/2, 949/2, 929/2 and 859/2%, respectively. These results show that the circulatory model is an appropriate model for assessing glucose turnover during an oral glucose load. # 2002 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Glucose kinetics; Tracer method; Oral glucose test; Mathematical models
1. Introduction The study of glucose metabolism often requires the calculation of glucose fluxes in non-steadystate conditions, for instance during an oral
* Corresponding author. Tel.: /39-049-829-5753; fax: /39049-829-5763 E-mail address: [email protected] (A. Mari).
glucose load or a glucose clamp. This calculation is based on tracer methodology, and requires a mathematical model. Unless the specific activity of the tracer is kept constant, which is often difficult or impossible, the accuracy of the calculated glucose fluxes depends on the model used. The most commonly used models are compartmental models with one or two compartments [1,2]. These models have proved useful in many
0169-2607/02/$ - see front matter # 2002 Elsevier Science Ireland Ltd. All rights reserved. PII: S 0 1 6 9 - 2 6 0 7 ( 0 2 ) 0 0 0 9 7 - 4
270
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
situations. However, discrepancies between the glucose fluxes calculated with these models and with independent methods have been found. Steele’s one-compartment model [1] has been shown to underestimate glucose production during a standard euglycemic hyperinsulinemic glucose clamp [3]. In experiments which simulate an oral load with an exogenous glucose infusion, both Steele’s model and the two-compartment model by Radziuk and colleagues [2] predicted the rate of appearance with some distortion [2,4]. During an oral glucose load, the rate of appearance calculated using the two-compartment model and the arteriovenous difference method were also somewhat different [5]. The problems encountered with the compartmental models may not originate entirely from the inadequacy of the models (cf. [5]). Nevertheless, compartmental models do not represent the physiological system as it is [6,7], and for this reason may introduce errors [8]. To avoid this drawback, circulatory models have been developed [8,9], which more appropriately represent the physiological system and are based on the solid principles of the theory developed for organ kinetics by Zierler [10,11]. Circulatory models have been used in various artificial non-steady-state conditions (e.g. [12,13]), but a model for general use has not been proposed and tested. In this work, we present a circulatory model for the calculation of nonsteady-state glucose fluxes in the general case. We evaluate the model performance in experiments in anesthetized rats that simulate an oral glucose load through an exogenous glucose infusion, by comparing the model-calculated rate of appearance with the known glucose infusion. We also compare the circulatory model performance with that of the more traditional compartmental approaches.
2. Methods 2.1. Experimental procedures 2.1.1. Basal glucose kinetics To determine the parameters for basal glucose kinetics, six adult male Sprague /Dawley rats were anesthetized with pentobarbitone sodium (Nem-
butal, Boehringer Ingelheim, NSW, Australia) administered by intraperitoneal injection (60 mg kg1 body weight). A catheter was placed in the right jugular vein (for tracer infusion) and left carotid artery (for blood sampling) [14]. The rats were given a single injection of [6-3H]glucose (10 mCi) and blood samples were taken at 0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 25, 30, 40, 50 and 60 min for analysis of plasma tracer. 2.1.2. Simulated oral load Another group of five male Sprague /Dawley rats were anesthetized and implanted with catheters as for the basal experiments. At time 0 min, a primed (185 mmol kg1) continuous infusion (3.8 mmol min 1 kg1) of [6,6-D2]glucose (Tracer Technologies, Somerville, MA, USA) in 0.9% NaCl was commenced. Three blood samples were taken at 60, 65 and 70 min for the estimation of the tracer-to-tracee ratio of [6,6-D2]glucose and plasma glucose. At 70 min, an infusion of 15 mCi ml1 [6-3H]glucose in 25% glucose commenced in an algorithm designed to mimic the appearance of gut-derived glucose following the administration of an oral glucose load. From time 70 to 90 min the glucose infusion rate was increased every minute to a maximum infusion rate of /48 mmol min1. This infusion rate was maintained for 15 min (time 90/105 min) after which the infusion rate was decreased every minute until it was stopped at 125 min. Blood samples (300 ml) were taken at 5 min intervals between 70 and 125 min for measurement of plasma glucose, the tracer-to-tracee ratio of [6,6-D2]glucose, and [6-3H]glucose specific activity. At the end of the study, timed collections of the [6,6-D2]glucose infusate were taken for accurate measurement of the constant [6,6-D2]glucose infusion rate. Timed collections of the variable glucose infusate were also taken to determine accurate time courses of the cold glucose delivery rate in each experiment. The [6-3H]glucose specific activity of the variable glucose infusate was also determined. 2.2. Analytical methods Plasma and infusate glucose levels, and plasma [6-3H]glucose tracer radioactivity and [6,6-D2]glu-
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
cose enrichment were measured as described previously [14,15]. [6,6-D2]glucose enrichment was analyzed using a gas chromatograph mass spectrometer (Shimadzu model GCSM-QP2000, Shimadzu Corporation, Kyoto, Japan) on ion monitoring mode to determine the relative intensity of the 98 and 100 molecular weight fragments i.e. (peak area at 100)/(peak area at 100/peak area at 98). A standard curve was run in parallel to convert relative intensity into relative tracer abundance (A ). The tracer-to-tracee ratio was then calculated as A /(1/A ). The study protocol was approved by the Royal Melbourne Hospital’s Animal Ethics Committee. 2.3. Circulatory model 2.3.1. Model structure Glucose kinetics are described using the circulatory model of Fig. 1, the mathematical theory of which has been presented previously [8,9,16]. In the model, the heart chambers and the lungs are represented in the heart /lungs block (upper block), while the remaining tissues are lumped into the periphery block (lower block). Blood flow for both blocks is cardiac output. In the experimental configuration adopted here, glucose is injected at the venous side, while blood is sampled at the arterial side (as indicated in Fig. 1). The two tissue blocks shown in Fig. 1 can be regarded as single inlet /single outlet organs, coupled in a feedback arrangement. The basis of the mathematical description of a block (as for an organ) is its impulse response [17]. The impulse response is the tracer efflux (concentration times blood flow) observed at the organ outlet after a bolus injection of a unit tracer dose at the organ inlet (this definition assumes tracer does not
recirculate). Cardiac output and the impulse responses of the blocks fully determine the circulatory model, and are specified in this section. Cardiac output was expressed as flow of blood per kilogram of body weight. Cardiac output was fixed to the value of 236 ml min 1 kg1, which is appropriate for anesthetized rats [18]. Since glucose concentration was measured in plasma and not in blood, this cardiac output value was corrected using the ratio between glucose concentration in whole blood and in plasma (0.53 in our rats) to ensure that the product of concentration and flow gives the actual glucose mass flux. A constant cardiac output value during the test was assumed (see Section 4). The impulse response of the heart /lungs block [rlung(t)] was modeled as a two-exponential function, which is the minimal representation ensuring a time course of the response in agreement with the experimental data. The impulse response starts from 0, rapidly increases to a peak value, and returns to zero as a single-exponential function. To simplify the notation, it is convenient to represent rlung(t ) as the convolution of two single-exponential functions: rlung (t)bebt vevt
(1)
where the symbol / denotes the convolution operator. As the coefficients of the single-exponential functions equal their exponents, the integral from 0 to of each exponential is 1. This also implies that the integral from 0 to of rlung(t) is 1, i.e. the glucose fractional extraction in the heart / lungs system is 0 [12]. The exponent that characterizes the rising phase of rlung(t) was fixed (b / 15 min 1, cf. [12]). The exponent of the decaying phase (v , min1) was calculated using the property that the mean transit time of rlung(t) is the ratio of volume to blood flow [17]. From this property, the following equation was obtained (cf. [12]): v
Fig. 1. The circulatory model.
271
bCO bVlung CO
(2)
where cardiac output (CO) and the heart /lungs glucose volume (Vlung) were assumed known
272
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
(CO /0.53 /236 ml min 1 kg1, see above; Vlung /17 ml kg1, see [12]). The impulse response of the periphery block [rper(t)] in the basal steady-state period can be expressed as [17]: rper (t) (1E)p(t)
(3)
where E (dimensionless) is the glucose fractional extraction and p (t) is the glucose transit time density function. The transit time density function p (t ) was modeled as a four-exponential function, starting from zero, rapidly increasing to a peak value, and returning to 0 as a three-exponential function. Similarly to rlung(t), we have represented p(t) as a convolution of a single-exponential function with a three-exponential function: p(t) gegt [w1 a1 ea1 t w2 a2 ea2 t (1w1 w2 )a3 ea3 t ]
(4)
where g , a1, a2, a3 (min 1) and w1, w2 (dimensionless) are parameters. In Eq. (4), the integral from 0 to of the threeexponential function in square brackets is w1/ w2/(1/w1/w2), i.e. 1. This ensures that the integral from 0 to of p (t) is also 1, as required for a transit time density function. The parameters w1 and w2 determine the contribution in the total response of the three exponential terms. As b in rlung(t), g determines the initial rising phase of p (t ). The parameter g was fixed (g/10 min 1, [12]), while a1, a2, a3, w1, w2 and E were estimated from the basal kinetic experiments. Eq. (3) can be extended to the non-steady-state [9]. In non-steady-state, the glucose fractional extraction E varies with time under the action of insulin, and the transit time density function p(t) may also change. However, in the present model we have assumed that the transit time density function is not affected by insulin. This is supported by the finding that the main effect of insulin is exerted on the fractional extraction of the periphery, while the transit time density function does not change substantially [8,12]. Thus, Eq. (3) was assumed to be valid also in the non-steady-
state period, with the same p (t) obtained from the basal period, and E (t) time-varying. To simulate the model of Fig. 1, the multiexponential equations for the impulse responses of the heart /lungs and periphery blocks are first represented as ordinary differential equations. These differential equations are then combined according to the feedback arrangement of Fig. 1, and then solved numerically using MATLAB (see Appendix A for details). 2.3.2. Analysis of basal steady-state experiments For the present circulatory model global identifiability of the model parameters is ensured in steady-state [16]. Mean parameters of the periphery transit time density function p (t) and the basal glucose fractional extraction E (six parameters) were obtained by least-squares fit of the mean [3-3H]glucose concentration in the basal kinetic experiments. The estimated parameters were subsequently used in the analysis of the non-steadystate data (a single parameter set for all nonsteady-state experiments). Parameter estimation was performed using a Levenberg /Marquardt least-squares algorithm (MATLAB function leastsq). 2.3.3. Analysis of non-steady-state experiments In non-steady-state, as p(t) was assumed to be time-invariant, the only time-varying parameter to be determined was E (t). E (t ) was assumed constant in the basal period preceding the start of glucose infusion, and represented as a piecewise constant function of time over 2 min time intervals thereafter (28 elements, E1 /E28, were necessary to cover the 55 min non-steady-state period). The Ek values were estimated from [6,6-2H2]glucose concentration using the model and a method traditionally used for deconvolution, as previously described [12]. In brief, the Ek values were estimated by least-squares fit of [6,6-2H2]glucose concentration, including in the least-squares a penalty term dependent on the second derivative of E (t ), calculated by forward differences of the Ek values. This additional term is necessary to obtain a smooth E (t), which would otherwise exhibit spurious oscillations.
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
273
2.3.4. Exogenous glucose appearance By dividing the [3-3H]glucose concentration by the specific activity of the glucose infusate, the glucose concentration component due to the exogenous glucose infusion was calculated (‘exogenous’ glucose concentration). From the exogenous glucose concentration and the model (using the tracer-determined E (t)), the rate of appearance of exogenous glucose was calculated using a deconvolution method [12]. For this purpose, we have approximated the exogenous glucose appearance rate as a piecewise constant function of time over 2 min time intervals, as E (t). 2.3.5. Glucose production The glucose concentration component due to glucose production (‘endogenous’ glucose concentration) was calculated by subtracting the exogenous glucose concentration from the total (measured) glucose concentration. From the endogenous glucose concentration and the model, glucose production was calculated using a deconvolution method, as for exogenous glucose. Glucose production was approximated as a constant value before t/70 min, and a piecewise constant function of time over 2 min time intervals thereafter. The calculated values were corrected for the non-negligible mass flux of the stable-isotope tracer (/3.8 mmol min 1 kg1), as the endogenous glucose concentration component includes the contribution of the tracer.
Fig. 2. Compartmental models. Top: single-compartment model; bottom: two-compartment model. See Eqs. (5) and (6) for the symbols.
We have used two values for the glucose distribution volume VS. The first value represents a standard volume choice (pool fraction of 0.65 with a total glucose volume of 250 ml kg1, i.e. VS /0.65 /250 ml kg1, [19]). The second value is the total glucose volume as calculated with the circulatory model in the basal period (VS /191 ml kg1). The time course of tracer and glucose concentration (exogenous or endogenous), and the derivative of specific activity were obtained from the data fit of the circulatory model (see Appendix A). This provides the necessary data smoothing and ensures that the differences between Steele’s model and the circulatory model are not due to differences in data smoothing and interpolation.
2.4. Compartmental models 2.4.1. One-compartment model The glucose rate of appearance [Ra(t )] was calculated every 2 min (time instants t0 /0, t1 / 2,. . . min) using Steele’s single-compartment model (Fig. 2 top, [1]), according to the equation (see [19]): Ra (tk )
R(t a k) a(tk )
VS
G(tk )a(t ˙ k) a(tk )
(5)
where Ra (tk ) is the tracer infusion rate, G (tk ) is arterial glucose concentration, a (tk ) is the glucose specific activity (the dot represents the time derivative), and VS is the compartment volume.
2.4.2. Two-compartment model The rate of glucose appearance was calculated using the two-compartment model of Fig. 2, according to the equations developed in [19]: Ra (tk )
Ra(tk )
G(tk )a(t ˙ k) k k V1 V1 12 21 a(tx ) a(tk ) k12 g(tk ) g(tk ) a(tk )
g(tk1 )b1 g(tk )b2 G(tk )b3 G(tk1 ); g(t0 )G(t0 )
274
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
g(tk1 ) b1 g(tk )b2 G(tk )b3 G(tk1 ); g(t0 )G(t0 )
(6)
where the star denotes the tracer, V1 is the volume of the first compartment, k21 and k12 are the intercompartmental rate constants (cf. Fig. 2), g(t) and g *(t ) are auxiliary variables (delayed glucose concentration), and the bk ’s (cf. [19]) are constants calculated from k12 and the calculation interval tk1/tk , which was 2 min as for Steele’s and the circulatory model. The parameters of the two-compartment model were obtained from the basal steady-state experiments, as for the circulatory model. The twocompartment model did not fit the initial part of the tracer disappearance curves accurately, as expected [20]. We have thus discarded the first 2 min in the analysis of the basal tracer curves. As for Steele’s model, concentrations of tracer and glucose, and the derivative of specific activity were obtained from the data fit of the circulatory model.
2.5. Evaluation of the model accuracy The accuracy of the models was evaluated by comparing the model-calculated rate of appearance with the actual glucose infusion rate, the error being the difference of the two. Besides the direct comparison of the time course of the mean rate of appearance, we have calculated the following indices of model performance: (1) percent recovery of the infused glucose dose, i.e. the ratio between the integral of the model-calculated rate of appearance and the total quantity of infused glucose, expressed in percent; (2) integral absolute error, i.e. the integral of the absolute value of the error on the rate of appearance; (3) cumulative distribution of the percent error (absolute value of the error divided by the actual glucose infusion rate). The cumulative distribution function was calculated by pooling the error values in all rats and at all time points. The frequency of occurrence was divided into ten percentiles.
3. Results 3.1. Basal glucose kinetics 3.1.1. Circulatory model The parameters of the circulatory model estimated from the tracer bolus injection in the basal state are reported in Table 1. The coefficients of variation of the parameters provided by the leastsquares algorithm were below 23%. The circulatory model predicted the [3-3H]glucose concentration curve accurately. 3.1.2. Two-compartment model The two-compartment model parameters were: V1 /119 ml kg 1; k01 /0.033; k21 /0.073; k12 / 0.11 min1; clearance /3.9 ml min 1 kg1; total volume /197 ml kg1. The coefficients of variation were below 20%. In the first 2 min, which were excluded from the fit, the two-compartment model underestimated [3-3H]glucose concentration, while in the remaining period the fit was accurate. 3.2. Non-steady-state glucose kinetics 3.2.1. Glucose and tracer concentrations Fig. 3 shows the mean [6,6-2H2]glucose concentration, total glucose concentration, and its exogenous and endogenous components. The solid lines represent the circulatory model fit. 3.2.2. Rate of appearance of exogenous glucose The mean rate of appearance of exogenous glucose calculated with the circulatory and the two-compartment model are shown in Fig. 4. The results for Steele’s model, obtained using the two glucose distribution volumes, are shown in Fig. 5. The standard Steele calculation (Fig. 5 top) was the least accurate, while the performance of the other models was not markedly different on average. All the models underestimated the downslope of the glucose rate of infusion. The twocompartment model most accurately predicted the downslope, while the circulatory model most accurately predicted the upslope. The performance of the modified Steele’s model with the volume set to the total glucose volume was similar to that of the circulatory model.
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
275
Table 1 Model parameters in the basal state for the anesthetized rats CO ml min 1 kg 1
125
E
0.0293
Heart /lungs
Periphery
b (min1)
v (min 1)
g (min1)
a1 (min1)
a2 (min 1)
a3 (min 1)
w1
w2
15.0
14.4
10.0
3.69
0.534
0.0788
0.684
0.265
CO, cardiac output; E , fractional extraction. For the remaining symbols see Eqs. (1) /(4). In the basal state, the clearance was 3.7 ml min 1 kg1 and the total volume 191 ml kg1. The parameter v was calculated from cardiac output and the volume of the heart / lungs block (Vlung /17 ml kg 1) using Eq. (2). The cardiac output value reported in the table is corrected for the distribution of glucose between plasma and red cells (CO /0.53 /236 /12 ml min 1 kg1, see Section 2).
Fig. 3. Mean (9/S.E.M.) tracer concentration, total glucose concentration, exogenous glucose concentration and endogenous glucose concentration. The solid line represents the circulatory model fit.
Table 2 reports the recovery of the infused glucose and the average absolute integral error for all models. Depending on the model, the correlation coefficient between the integral of the calculated rate of appearance and the total quantity of infused glucose was /0.85 /0.88, with a borderline significance (P :/0.05 /0.07). These integral accuracy indexes basically confirmed the results of Figs. 4 and 5.
Fig. 4. Mean (9/S.E.M.) exogenous glucose rate of appearance calculated with the circulatory (top) and the two-compartment model (bottom). The solid lines represent the true mean glucose infusion rate (standard errors are omitted for clarity).
The distribution of the pooled relative error is shown in Fig. 6. The curves show the percentage of the rate of appearance values that are affected by an error which does not exceed a specified value. For instance, for the circulatory model 60% of the calculated rates of appearance (read off the ordinate) are affected by an error which is not greater than 12% (read off the abscissa). In this representation, the curves of the more accurate models lie above those of the less accurate models. According to this analysis, which shows the error
276
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
Fig. 6. Cumulative distribution of the pooled relative error. Results for all models are reported, as indicated in the figure. The ordinate shows the percentage of the calculated rates of appearance that do not have an error which exceeds the level shown on the abscissa. For clarity, only 90% of the calculated rates are included. The error value on the abscissa which corresponds to 100% of the calculated rates (maximal error) is in fact much higher: circulatory model, 145%; two-compartment model, 121%; Steele’s model, 107% (V/162.5 ml kg 1) and 134% (V /191 ml kg 1). Fig. 5. Mean (9/S.E.M.) exogenous glucose rate of appearance calculated with Steele’s model and the two glucose volume values indicated in the figure. The solid lines represent the true mean glucose infusion rate (standard errors are omitted for clarity).
in 90% of the calculated points for clarity, the overall performance of the circulatory model is the best and that of the standard Steele’s model the worst. However, with the exception of the standard Steele’s model, the differences are small. 3.2.3. Glucose production Fig. 7 shows mean glucose production calculated using the circulatory model. Glucose production decreased to 0 at 25 min after the start of glucose infusion, and slowly returned towards the
Fig. 7. Mean (9/S.E.M.) glucose production calculated with the circulatory model.
Table 2 Recovery and integral absolute error
Recovery (%) Integral absolute error (mmol kg 1) a
Circulatory model
Two-compartment model
Steele’s model V/162.5 ml kg 1
Steele’s model V/191 ml kg 1
939/2 7099/89
949/2 7179/72
859/2 8629/106a
929/2 7719/93
P B/0.05 or less vs. circulatory model (paired t -test). The recovery was significantly less than 100% for all models (P B/0.05 or less, one-sample t -test).
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
basal value thereafter. Glucose production as calculated using the compartmental models was similar to that reported in Fig. 7. The largest difference was observed with the standard Steele’s model, which yielded a higher glucose production on the average, but the maximum difference was less that 4 mmol min 1 kg1.
4. Discussion In our experiments, the circulatory model predicted the actual glucose infusion rate with an average error of about 10%. The rising phase of the glucose infusion was precisely matched, while the falling phase was underestimated. The biggest error was observed in this period, but on average the error was about 15%, and in 90% of all the calculated rates of appearance it did not exceed 30%. According to the distribution of the error on the rate of appearance (Fig. 6, Table 2), the performance of the circulatory model was the best. However, none of the models matched the glucose infusion rate accurately, and the evaluation of the relative performance is thus somewhat subjective. Indeed, the accuracy of the two-compartment model and Steele’s model with the larger volume was only slightly inferior to that of the circulatory model. As concerns the circulatory model, the error on the glucose rate of appearance may originate from two major sources. The first is the use of mean instead of individual parameters. Non-steady-state calculations are made using a fixed parameter set for all rats, while it is expected that the individual model parameters are somewhat different from rat to rat. This may introduce an error in the individual rates of appearance. We could not evaluate the magnitude of this error because experimental limitations prevented us from assessing the basal parameters in each rat used for the non-steady-state experiments. However, it is unlikely that this error is also significant for the mean rate of appearance, as the interindividual parameter differences are random and tend to cancel out in the mean. The second error concerns limited information on the system. The circulatory model is a physical representation of the system, but
277
assumes a fixed cardiac output and a timeinvariant distribution of the transit times p(t) of the periphery because the little is known about the changes that the cardiac output and p (t ) undergo in this and other experimental conditions. During glucose challenges or hyperinsulinemia in conscious human subjects, changes in cardiac output do not exceed 20% [21 /23]. Changes from basal in p(t) have been determined during hyperinsulinemic euglycemic glucose clamps, directly in dog [12] and indirectly in man [8], and in both cases are not substantial. In hyperinsulinemic hyperglycemic glucose clamps in anesthetized rats, we have also found small changes in p (t ) (/20% difference in mean transit time, results not shown). These findings justify the assumption of this study of constant cardiac output and time-invariant p (t). Additional experimental and theoretical work are, however, needed to resolve this issue, and to clarify whether more precise assumptions result in a better model performance, or if factors other than the model are the cause of the observed discrepancies, as pointed out previously (e.g. [5]). Unfortunately, the determination of the time course of cardiac output and the changes of p (t ) is a difficult problem to solve, and this limits the possibilities of model improvement. In regard to cardiac output, we have assumed the most appropriate value for anesthetized rats [18]. In principle, cardiac output cannot be ignored, because it does have a role in glucose kinetics and thus in the calculation of glucose fluxes. However, the calculations are not very sensitive to the cardiac output value. Only with rapid changes of fluxes, concentrations and specific activity, which are not observed in the present experiments, is it expected that the calculated glucose fluxes are influenced by the cardiac output value. It is in fact known that in steady-state, or with negligible changes in specific activity, cardiac output is not relevant for calculating glucose turnover. The compartmental models considered here performed differently. The two-compartment model performance was only slightly inferior to that of the circulatory model. The standard Steele equation markedly distorted the shape of the glucose infusion rate, but the poor performance
278
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
parameters of Tables 1 and 3 to experimental conditions in which the cardiac output value is different from that reported here. In conclusion, we have presented a circulatory model for the calculation of the rate of appearance in non-steady-state conditions and we have assessed its accuracy experimentally. The circulatory model gives a physiological representation for glucose kinetics, and its accuracy is to some extent superior to that of the more standard Steele and two-compartment models.
of this model could be corrected by using a larger glucose volume. This stratagem has theoretical support when the specific activity changes slowly [24]. However, the use of larger glucose volume does not ensure sufficient accuracy in general. For instance, the calculation of glucose production during a standard glucose clamp is not accurate with Steele’s model, regardless of the glucose volume used [24]. Therefore, in this experimental condition compartmental models are sufficiently accurate, but, as demonstrated for the Steele’s model, it is not necessarily true that the error due to the incorrect representation of the physiological system is small in all situations. The circulatory model presented here is suitable for application to general non-steady-state conditions, as is the Steele’s one-compartment model and the two-compartment model. To use the circulatory model, the cardiac output value and the parameters of p (t) are required. It is clearly important that the parameters used are appropriate for the specific experimental condition. If possible, the model parameters should be estimated in each individual experiment, using for instance the tracer equilibration period that normally precedes the actual non-steady-state phase. In this work, however, a single set of parameters (Table 1) was used for all rats because individual estimates could not be obtained. These values are appropriate for anesthetized rats. As it is useful to provide analogous values for humans, we report in Table 3 a set of parameters derived from a study in lean and obese subjects [13]. We also present in Appendix B some results useful to extend the
Acknowledgements This study was supported in part by a grant from the Italian National Research Council.
Appendix A In this appendix, the differential equations of the model are derived. Additional details on the equations can be found in [12], which presents a similar model using the same notation. Basic differential equation The key for the transformation of the impulse responses (Eqs. (1), (3) and (4)) into differential equations is the representation of the convolution by means of a differential equation. If the impulse response of the heart /lungs block were the single-
Table 3 Model parameters in the basal state for humans (from [13]) CO (ml min1 kg 1)
E
Heart /lungs b (min
66
0.0223
15.0
1
)
Periphery v (min 5.16
1
)
g (min 1)
a1 (min1)
a2 (min 1)
a3 (min1)
w1
10.0
3.05
0.382
0.0388
0.430
w2 0.505 1
CO, cardiac output; E , fractional extraction. For the remaining symbols see Eqs. (1) /(4). The total volume was 223 ml kg . The parameter v was calculated from cardiac output and the volume of the heart /lungs block (Vlung /17 ml kg1) using Eq. (2). The cardiac output value reported in the table is corrected for the distribution of glucose between plasma and red cells (see Section 2). In [13], cardiac output (flow of blood) was assumed to be 78 ml min 1 kg1, and the correction factor 0.84 (CO /0.84/78/66 ml min 1 kg1). In [13], cardiac output was expressed in ml min 1 m 2 (cardiac index), and is converted here into ml min 1 kg 1 assuming a body surface area to body weight ratio of 1.7/70.
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
exponential function k elt , the relationship between the input (glucose concentration in the right atrium, Gra(t)) and the output (arterial glucose concentration, Ga(t)) would be the convolution Ga (t)kelt Gra (t)
x˙ 3 (t)a1 x3 (t)w1 (1E)Ga (t) x˙ 4 (t)a2 x4 (t)w2 (1E)Ga (t) x˙ 5 (t)a3 x5 (t)(1w1 w2 )(1E)Ga (t)
(A1)
This convolution can be represented by means of the differential equation x(t)lx(t)G ˙ ra (t) Ga (t)kx(t)
(A2)
where x (t) is an additional variable (state variable), and the dot represents the time derivative. Eq. (A2) is a state space representation of a linear system, in a canonical form. Its state variables do not represent physical quantities. The initial condition for Eq. (A2) (and the successive equations) depends on whether endogenous glucose, oral glucose or the tracer is considered. For oral glucose and the tracer, which are not present in the system before time 0, x (0) /0. For endogenous glucose, x (0) is the steady-state value, i.e. the solution of the algebraic equations obtained from Eq. (A2) by letting x(t) ˙ 0:/
x˙ 6 (t)gx6 (t)a1 x3 (t)a2 x4 (t)a3 x5 (t) Gmv (t)gx6 (t)
Because the two-exponential impulse response of the heart /lungs block (Eq. (1)) is represented as a convolution of two single-exponential functions, the relationship between Gra(t ) and Ga(t) can be expressed as a cascade of two differential equations like Eq. (A2), i.e.
(A4)
where Gmv(t) is mixed-venous glucose concentration (output of the periphery block), and x3 /x6 are state variables. In Eq. (A4), the three exponential terms in the square brackets in Eq. (4) are represented by a sum of three differential equations, while the last two equations express the convolution with the exponential function geg . t
Whole-body differential equations In the feedback arrangement of Fig. 1, the input of the heart /lungs block is the output of the periphery block [Gmv(t)] plus a term representing the appearance of glucose, which is Ra(t)/CO. Therefore, the whole set of differential equations describing the circulatory model is: x˙ 1 (t)vx1 (t)gx6 (t)
Heart /lung differential equations
Ra (t) CO
x˙ 2 (t)bx2 (t)vx1 (t) x˙ 3 (t)a1 x3 (t)w1 (1E)bx2 (t) x˙ 4 (t)a2 x4 (t)w2 (1E)bx2 (t) x˙ 5 (t)a3 x5 (t)(1w2 w2 )(1E)bx2 (t)
x˙ 1 vx1 (t)Gra (t)
x˙ 6 (t)gx6 (t)a1 x3 (t)a2 x4 (t)a3 x5 (t)
x˙ 2 bx2 (t)vx1 (t)
Ga (t) bx2 (t)
Ga (t)bx2 (t)
279
(A5)
(A3)
where x1 and x2 are state variables. Periphery differential equations Similarly, the impulse response of the periphery block (Eqs. (3) and (4)) is represented by the four differential equations:
Numerical solution Eq. (A5) is a standard representation in state space form of a linear system of differential equations. The solution of the system yields the glucose concentration (tracer, exogenous or endogenous) that corresponds to a given rate of
280
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
appearance. Eq. (A5) is valid both in steady and non-steady-state. In steady-state E is constant, and Eq. (A5) represents a linear time-invariant system. In non-steady-state E is time-varying, and the system is time-varying. However, as E(t) is assumed to be piecewise constant, the system is in reality piecewise time-invariant, which facilitates the calculation of the solution. To solve the differential equations numerically, Eq. (A5) was expressed in matrix form as: x(t)Ax(t)BR ˙ a (t) Ga (t)Cx(t)
(A6)
where A, B and C are the standard system matrices, x (t) is the state vector in which the state variables xk (t) are stacked, and the initial value x (0) is as discussed for Eq. (A2). Because E (t) is assumed to be piecewise constant, the matrix A is also piecewise constant. Eq. (A6) was integrated by updating the state vector step by step over the time grid resulting from the union of the data time points and the time grid for E (t) and Ra(t). In this way, in each step A, B and C are constant matrices, and x (tk1) was updated from x (tk) after calculation of the matrix exponential eA(tk1tk ) (MATLAB function expm), using the standard equations for state update with piecewise constant inputs [25]. Except for Ra(t ) and E (t), all the parameters of Eq. (A5) are reported in Table 1 (Table 3 for humans). /
Appendix B This appendix briefly addresses the problem of adapting the parameters of Tables 1 and 3, which are appropriate for anesthetized rats or lean and obese humans, to experimental conditions in which the cardiac output value is different from that reported here. The glucose distribution volume in the periphery block is the product of cardiac output, 1-fractional extraction and the mean transit time of the periphery [8]. The glucose distribution volume is likely to be substantially independent from the cardiac output value. If the circulatory model is used in conditions in which cardiac output is different (e.g. conscious vs. anesthetized rats), it is necessary to change the periphery mean transit time in order that the glucose volume remains constant. This is achieved by adjusting the exponents ak of the periphery impulse response as follows. The glucose volume is calculated from the periphery impulse response (Eq. (4)) as the sum of volume terms Vk , one for each exponential (the volume contribution due to g is negligible): Vk CO(1E)
(B1)
To maintain the total glucose volume constant, it is logical to increase in each volume term the exponents ak in proportion to the cardiac output value. Thus, in an experimental condition in which cardiac output is CO?, the new exponents, ak? ?, will be:
Specific activity derivative
ak? The numerical solution of Eq. (A5) provides a smooth time course of Ga(t) (glucose and tracer concentration) and also of its derivative, G˙ a (t): G˙ a (t) is in fact proportional to x˙ 2 (t) (Eq. (A5), last equation), and is thus a linear combination of x1(t) and x2(t) (Eq. (A5), second equation), which are smooth. From the glucose and tracer smooth Ga(t) and G˙ a (t); a smooth time course of the specific activity derivative can be obtained. The smoothed variables were used with the compartmental models (Eqs. (5) and (6)).
wk ak
CO? ak CO
(B2)
where CO and ak are the values reported in Table 1 or Table 3. Eq. (B2) is an approximate correction, which assumes that cardiac output does not affect substantially the distribution of glucose. If there are reasons to think that cardiac output affects glucose distribution, the model parameters should be determined experimentally. The cardiac output problem is not apparent using compartmental models. However, this does
A. Mari et al. / Computer Methods and Programs in Biomedicine 71 (2003) 269 /281
not imply that the compartmental model parameters are unaffected by cardiac output. It is in fact expected that in conditions in which cardiac output is different, also the compartmental model parameters (e.g. k12 and k21) differ [8].
References [1] R. Steele, Influences of glucose loading and of injected insulin on hepatic glucose output, Annals of the New York Academy of Sciences 82 (1959) 420 /430. [2] J. Radziuk, K.H. Norwich, M. Vranic, Experimental validation of measurements of glucose turnover in nonsteady-state, American Journal of Physiology 234 (1978) E84 /E93. [3] D.T. Finegood, R.N. Bergman, M. Vranic, Estimation of endogenous glucose production during hyperinsulinemiceuglycemic glucose clamps. Comparison of unlabeled and labeled exogenous glucose infusates, Diabetes 36 (1987) 914 /924. [4] J. Proietto, F. Rohner-Jeanrenaud, E. Ionescu, J. Terrettaz, J.F. Sauter, B. Jeanrenaud, Non-steady-state measurement of glucose turnover in rats by using a onecompartment model, American Journal of Physiology 252 (1987) E77 /E84. [5] A. Mari, J. Wahren, R.A. DeFronzo, E. Ferrannini, Glucose absorption and production following oral glucose: comparison of compartmental and arterio-venous difference methods, Metabolism 43 (1994) 1419 /1425. [6] K. Zierler, A critique of compartmental analysis, Annual Review of Biophysics and Bioengineering 10 (1981) 531 / 562. [7] K. Zierler, Whole body glucose metabolism, American Journal of Physiology 276 (1999) E409 /E426. [8] A. Mari, Circulatory models of intact-body kinetics and their relationship with compartmental and noncompartmental analysis, Journal of Theoretical Biology 160 (1993) 509 /531. [9] A. Mari, Calculation of organ and whole-body uptake and production with the impulse response approach, Journal of Theoretical Biology 174 (1995) 341 /353. [10] P. Meier, K.L. Zierler, On the theory of the indicatordilution method for measurement of blood flow and volume, Journal of Applied Physiology 6 (1954) 731 /744. [11] K.L. Zierler, Theory of the use of arteriovenous concentration differences for measuring metabolism in steady and non-steady-states, Journal of Clinical Investigation 40 (1961) 2111 /2125.
281
[12] O.P. McGuinness, A. Mari, Assessment of insulin action on glucose uptake and production during a euglycemichyperinsulinemic clamp in dog: a new kinetic analysis, Metabolism 46 (1997) 1116 /1127. [13] A. Natali, A. Gastaldelli, S. Camastra, A.M. Sironi, E. Toschi, A. Masoni, E. Ferrannini, A. Mari, Dose-response characteristics of insulin action on glucose metabolism: a non-steady-state approach, American Journal of Physiology 278 (2000) E794 /E801. [14] C.J. Nolan, J. Proietto, The feto-placental glucose steal phenomenon is a major cause of maternal metabolic adaptation during late pregnancy in the rat, Diabetologia 37 (1994) 976 /984. [15] A. Thorburn, J. Muir, J. Proietto, Carbohydrate fermentation decreases hepatic glucose output in healthy subjects, Metabolism 42 (1993) 780 /785. [16] A. Mari, Determination of the single-pass impulse response of the body tissues with circulatory models, IEEE Transactions on Biomedical Engineering 42 (1995) 304 / 312. [17] N.A. Lassen, W. Perl, Tracer Kinetic Methods in Medical Physiology, Raven Press, New York, 1979. [18] M.W. Crawford, J. Lerman, V. Saldivia, F.J. Carmichael, Hemodynamic and organ blood flow responses to halothane and sevoflurane anesthesia during spontaneous ventilation, Anesthesia and Analgesia 75 (1992) 1000 / 1006. [19] A. Mari, Estimation of the rate of appearance in the nonsteady-state with a two-compartment model, American Journal of Physiology 263 (1992) E400 /E415. [20] C. Cobelli, G. Toffolo, E. Ferrannini, A model of glucose kinetics and their control by insulin, compartmental and noncompartmental approaches, Mathematical Biosciences 72 (1984) 291 /315. [21] S. Jern, Effects of acute carbohydrate administration on central and peripheral hemodynamic responses to mental stress, Hypertension 18 (1991) 790 /797. [22] A.D. Baron, Hemodynamic actions of insulin, American Journal of Physiology 267 (1994) E187 /E202. [23] K. Masuo, H. Mikami, T. Ogihara, M.L. Tuck, Mechanisms mediating postprandial blood pressure reduction in young and elderly subjects, American Journal of Hypertension 9 (1996) 536 /544. [24] C. Cobelli, A. Mari, E. Ferrannini, Non-steady-state: error analysis of Steele’s model and developments for glucose kinetics, American Journal of Physiology 252 (1987) E679 /E689. [25] C.T. Chen, Linear system theory and design, third ed. (Chapter 4), Oxford University Press, New York, 1999.