A Physiologically Based Approach for the Estimation of Recirculatory Parameters

Indicator dilution methods have been very useful in circulatory studies of physiological systems because they can be used to provide estimates of blood flow and volumes of the intra- and extravascular space (Chinard et al., 1962; Goresky, 1962; Zierler, 1962, 1965, 2000; Maseri et al., 1970; Lund-Johansen, 1990; van der Hoeven and Olsman, 2000). The techniques used rely on an intravenous bolus administration of a tracer at a specific site of the circulatory system, followed by frequent blood sampling at a distal site. In reality, the concentration-time (C,t) profile of the first pass of the tracer through the circulation is required for the valid estimation of the physiological circulatory parameters. However, this profile is not available in practice due to recirculation that distorts the downslope of the observed curve; this leads to an overestimation of the area under the concentration-time curve, which is the key parameter in this field of studies. To overcome this drawback, various methods have been proposed. The Stewart-Hamilton extrapolation is a common technique assuming a monoexponential decay of the concentration after a time point (Millard, 1997), whereas other procedures suggest the fitting of the declining part of the curve with the gamma, lognormal, and local density random walk probability distributions (Wise, 1966). There have also been proposed a number of easy-to-apply methods and deconvolution approaches (Stephenson, 1948; Bradley, 1969; Millard, 1997). To avoid the second passage of the indicator from the sampling point, various techniques using nonrecirculatory “tracers” have been developed, namely, inert gases, like hydrogen and freon-22, as well as Li⁺ ions (Wittenberg et al., 1971; Jones et al., 1991; Band et al., 1997). In the same vein, thermodilution constitutes another alternative where a certain amount of iced solution is injected instead of a tracer (Sibbald et al., 1983; Mihm et al., 1987; Conway and Lund-Johansen, 1990; Jansen, 1994; Picker et al., 2001). Additionally, various noninvasive techniques have been developed for the measurement of cardiac output (Ludman et al., 1993; Kasravi et al., 1998).

In this study, we present a novel approach for the estimation of physiological recirculatory model parameters. Our proposed methodology uses a dispersion-convection partial differential equation (PDE), which was derived from the fractal geometry of the vascular tree and the hydrodynamics of the blood flow (Dokoumetzidis and Macheras, 2003). The PDE was found to adequately describe the initial spatiotemporal profile of tracers in the circulatory system (Dokoumetzidis and Macheras, 2003). In the present study, the PDE is used to obtain virtual C,t data for an “ideal, nonrecirculatory tracer” following a two-step procedure. First, the PDE is fitted to the actual experimental C,t data of the tracer to derive estimates for the PDE model parameters. Subsequently, the PDE is applied again, using the parameter values derived previously, but now the model is properly altered to generate virtual C,t data corresponding to those that would have been obtained if the tracer was nonrecirculatory, i.e., if the experimental conditions were ideal. The appropriate alterations of the PDE are implemented by assigning the injection and sampling points in accord with the ideal “experimental” conditions; besides, the recirculatory effect is deliberately suppressed imposing full elimination of the tracer after its passage from the sampling point. Then, the so-derived single-pass C,t curve is used for the estimation of physiological circulatory parameters. To check the validity of the proposed approach, the two-step procedure was applied to eight C,t data sets obtained from literature (Krejcie et al., 1996, 1997; Niemann et al., 2000).

Materials and Methods

Cardiac Output, Partial Blood Volumes, and Mean Transit Times

Indicator dilution studies can be used to calculate a variety of recirculatory physiological parameters such as cardiac output (CO), partial blood volume (V, “needle to needle”), tissue volume, mean transit time (MTT), and liver function.

A common method for the estimation of CO is based on eq. 1 (Zierler, 1962, 2000; Lund-Johansen, 1990) where Dose is the amount of the injected tracer (e.g., indocyanine green) and AUC is the area under the concentration-time curve, which ideally corresponds to the first passage of the tracer from the sampling point. However, several assumptions are required for the valid use of eq. 1 (Goresky, 1962; Zierler, 1962). It is necessary that flow and volume are constant during the period of measurement. Besides, eq. 1 presupposes mass conservation, i.e., the entire amount of the injected indicator should pass from the sampling point sooner or later.

Also, MTT of a tracer when passing from injection to the sampling point are calculated using eq. 2 where AUMC and AUC refer to area under the first moment concentration-time curve and area under the curve, respectively.

The needle to needle volume (V), i.e., the volume between the injection and the sampling point is calculated with the application of the classic relationship (Zierler, 1962, 2000; Godje et al., 1998, 2000): where F is the blood flow. Also, concomitant administration of two or more tracers can lead to the determination of the disposition for each one of them based on the corresponding MTT values (Godje et al., 1998, 2000; Xu et al., 1990).

A valid application of the eqs. 1, 2, 3 presupposes that the concentration-time profile of the tracer comprises only its first passage through the sampling point. However, this profile is not available in practice due to recirculation that distorts the downslope of the observed curve and leads to an overestimation of the AUC. Various methods have been published regarding this drawback, as mentioned in the Introduction.

Mathematical Model for the Initial Mixing of Tracers

Recently, a model for the mammalian circulatory system, built in entirely physiological terms was introduced (Dokoumetzidis and Macheras, 2003). The model consists of a ring shaped, one-dimensional tube that corresponds to the arterial and venular trees successively. It appropriately describes the heterogeneous concentration-time profile of the initial mixing phase of a tracer in the vascular tree. The model is based on the scaling laws of the fractal structure of the vascular tree, which is reduced to an equivalent one-dimensional tube. The two ends of the tube are actually considered to be connected and correspond to the left part of the heart (Fig. 1). The spatiotemporal profile of the intravascular tracer C(x,t) is described by a dispersion-convection PDE with piecewise constant coefficients: where x is the dimensionless spatial coordinate that takes values from 0 to 1, which correspond to the various regions of the circulatory system. D(x) is the normalized dispersion coefficient, with units per second, which is a piecewise constant, taking the values x_c is the location of the capillaries and x_p is the location of the right part of the heart (Fig. 1). K represents the uptake rate constant; Θ (x) = 1 only along the region where uptake is considered and Θ (x) = 0 otherwise. U is the normalized blood velocity with units per second.

It is worth mentioning that the spatial coordinate, x, does not correspond to actual length because it is not expressed in linear scale but comes from a volume preserving coordinate transformation. The latter progressively stretches the lengths at the lower levels of the tree (toward capillaries), allowing the velocity, as well as the dispersion coefficient in the new coordinate system, to be constant throughout the entire tree, instead of being spatially dependent as they are in the ordinary coordinate system before the transformation. In a way, it is like using a logarithmic scale for the spatial coordinate. The entire derivation is presented in Dokoumetzidis and Macheras (2003).

The boundary condition that implements the recirculation of the blood connecting the two ends of the tube and forming a ring shaped tube is The initial condition describing the intravenous bolus injection of the tracer is a “thin” Gaussian function, namely, where Q(0) is the tracer dose, V is the blood volume, x_inj is the location of the injection site, and b is a shape parameter with a large value (b = 10⁵).

Equation 4 can be used to simulate C(x,t) of the concentration of a tracer when a dose Q(0) is injected at the site x_inj. However, a more useful output is the C(x_sam,t), which is the time profile of the concentration at the sampling site, x_sam. This last function can be used to fit the model to experimental data to estimate its parameters. Alternatively, it may be used to generate C,t data for a given set of parameters, considering a variety of virtual sampling and injection sites, implementing various protocol designs.

A Two-Step Procedure for the Estimation of Recirculatory Parameters

The estimation of CO, MTT, and partial blood volumes is based on eqs. 1, 2, 3 in conjunction with C,t data of a tracer after its injection to a specific point and sampling at a distal point in the circulatory system. However, the recirculation of tracer affects the downslope of the C-t curve and leads to erroneous estimates for the recirculatory parameters. In the present study, this drawback is abolished with the application of the proposed two-step procedure using eq. 4, which allows the generation of virtual C,t data without the presence of a second passage of the tracer through the sampling point. Initially, eq. 4 is fitted to the experimental C,t data of the tracer. This step allows the estimation of the necessary parameters of eq. 4, which are held fixed in the subsequent step. Then, eq. 4 is used again to generate virtual C,t data, which would correspond to a nonrecirculatory tracer. For this purpose, the parameters of eq. 4 are altered appropriately to suppress the second passage of the tracer through the sampling point, namely, 1) the location of the injection and sampling sites are altered in accord with the ideal “experimental” conditions, and 2) the recirculatory effect is suppressed imposing total elimination of the tracer after its passage from the sampling point. The so-derived single-pass C,t curve is then used for the estimation of the circulatory parameters.

It should be noted that eq. 4 possesses favorable characteristics for the estimation of recirculatory parameters (CO, MTT, partial blood volumes). One of the basic elements of the PDE model is that flow, F, is considered to be constant and equal to cardiac output. For this reason, when we refer to flow we mean the cardiac output and we will use CO throughout. Also, eq. 4 offers the possibilities that the injection and sampling points can be set in accord with an ideal implementation of a virtual experiment. For example, specific values are assigned to the injection, x_inj and sampling, x_sam points to mimic injection to right or left atrium and sampling to the ascending aorta, respectively, when the PDE is used to generate data for CO measurement. In parallel, to avoid the recirculation effect we set a high value to the uptake rate constant, as in the case of an indicator that is completely removed after the first passage from the sampling point. The generated concentration-time curve corresponds to the case of an indicator that was injected/detected in ideal points and in the same time does not recirculate. Thus, the required conditions for the estimation of cardiac output are simulated relying only on the PDE model and one easy-to-apply experimental process. Multiplying the reciprocal of the first pass area (AUC) of the simulated C,t curve with Dose, an estimate for CO is obtained (eq. 1).

Fittings of the PDE to Tracer Data. To perform the tasks described previously, the methodology presented was applied to C,t data obtained from literature characterizing the initial mixing of a variety of tracers injected to humans and dogs (Krejcie et al., 1996, 1997; Niemann et al., 2000). Before the application of the two-step procedure, a preliminary investigation was performed by fitting eq. 4 to all available data to identify parameters that can be held fixed in the model. This exercise was applied to reduce the number of model parameters that need to be estimated in the fittings. The parameters of the model are the site of the right heart, x_p; the site of the peripheral capillaries, x_c; the total blood volume, V; the velocity of the blood inside the tube, U; the dispersion coefficients, D₁ and D₂, of the arterial and pulmonary trees, respectively; and the uptake rate constant, K. Parameters that were assigned arbitrary values are the length of the uptake site, 0.02, and the shape parameter of the initial condition b = 10⁵. The PDE model was fitted to all available data and estimates for all the seven parameters were derived. This investigation revealed that the parameters x_c and x_p can depend only on the physiology of the species and may be kept constant for each species. Thus, two models were determined one for humans and one for dogs, each one left with five parameters to be estimated. The only difference of the two models is the values of the parameters x_c and x_p.

Eight data sets taken from literature were used to apply the methodology. The data were obtained by digitization from the published figures. The data sets are divided into three groups of different subjects enumerated according to their source: I and II; III, IV, and V; and VI, VII, and VIII. They include four different tracers: indocyanine green (ICG) in humans (data sets I and II) and dogs (data sets III and VI), lidocaine (data set V), antipyrine (data sets IV and VII) and inulin (data set VIII) in dogs (Krejcie et al., 1996, 1997; Niemann et al., 2000). Each tracer has different properties regarding the extent of distribution and the sites of uptake. For extravascular tracers, the term uptake does not necessarily mean elimination from the body, it only means the blood concentration diminution of the tracer. For fitting purposes, the model was adapted appropriately for each one of the tracers examined. The variants of the model used for each specific tracer are presented in Table 1. Each model variant may be considered either for humans or dogs by changing the aforementioned x_c and x_p values accordingly.

ICG is distributed almost exclusively within the intravascular space and exhibits liver uptake. According to the model used for the description of ICG kinetics (Table 1, model F1), the injection site is located at a major vein (Fig. 1, site γ) and the sampling site is in a major artery (Fig. 1, site α₂). The hepatic uptake of ICG takes place at the capillaries (Fig. 1, site β). Because of its solely intravascular distribution, ICG was used to estimate CO, partial blood volumes, and MTT.

Inulin is a polysaccharide that is distributed to extravascular fluids by free water diffusion (Krejcie et al., 1996). However, before its recirculation, inulin is confined to the intravascular space. The variant of the model used is identical to the one used for ICG (Table 1, model F1). Only MTT is estimated from inulin data.

Lidocaine is a basic amine used as a tissue water tracer, which undergoes extensive pulmonary uptake (Krejcie et al., 1997). The model variant used for lidocaine (Table 1, model F2) considers injection at point γ, sampling at point α₂ and uptake from region ϵ (lungs). Lidocaine data were used to estimate MTT only.

Antipyrine like other lipid-soluble drugs is a total body water tracer and can easily penetrate cell membranes (Krejcie et al., 1996, 1997). The model variant used for antipyrine (Table 1, model F3) considers injection at γ, sampling at α₂, and uptake from sites β and ϵ (Fig. 1). Like with lidocaine, only MTT values were estimated from antipyrine data.

Virtual Experiments of Ideal Tracers. Based on the fittings results using the real data of the above-mentioned tracers, it is possible to use other variants of the model to simulate virtual experiments of ideal nonrecirculatory tracers that are injected and sampled at ideal sites for the specific flow, volume, and MTT that need to be calculated. All the rest model parameters are kept to the values estimated by the fitting. The suppression of recirculation property is achieved by making the entire quantity of the virtually injected tracer disappear from the system right after the first pass by attenuating the uptake rate constant as well as by extending the elimination region which is situated either at site β or ϵ (Fig. 1). More specifically, five model variants (Table 1, S1–S5) were used to conduct virtual experiments, depending on the quantities desired to be calculated. Thus, a C,t curve is simulated, which represents the first passage of the tracer. Multiplication of the reciprocal of the area under this curve with the value of the administered dose leads to the estimation of CO (eq. 1). The mean transit time can be calculated using eq. 2, whereas the respective needle to needle volume is calculated from eq. 3. However, the uncertainty regarding the exact needle positions (x_inj, x_sam) sets limitations when the role of these positions is crucial, e.g., models S2 and S3 (Fig. 1). Accordingly, the models S2 and S3 are mostly presented to exhibit the flexibility of the approach.

Model S1 is used for the calculation of the total cardiac output, CO_S1; the pulmonary volume of distribution, including the heart, V_S1, of the tracer; and the pulmonary mean transit time, including the heart, MTT_S1. It includes injection a little before right atrium (at γ of Fig. 1), sampling a little after left ventricle (at α₁ of Fig. 1) and attenuated elimination along the extended region β_ext (Fig. 1). This model was applied to ICG data for the estimation of CO_S1, V_S1, and MTT_S1 and to all other tracer data for the estimation of MTT_S1.

Model S2 is used to calculate the cardiac output of the right heart, CO_S2, the volume, V_S2 of the right heart, and the respective MTT_S2 for the tracer. It includes injection just before right atrium (at γ of Fig. 1), sampling just after right ventricle (at δ of Fig. 1), and attenuated elimination along the extended region β_ext (Fig. 1). This model was applied to ICG data for the estimation of CO_S2, V_S2, and MTT_S2 and to all other tracer data for the estimation of MTT_S2.

Model S3 is the equivalent of S2 for the left heart. It includes injection a little before left atrium (at ζ of Fig. 1), sampling a little after right ventricle (at α₁ of Fig. 1), and attenuated elimination along the extended region β_ext (Fig. 1). This model was applied to ICG data for the estimation of the cardiac output of the left heart CO_S3, the volume of the left heart, V_S3, and the respective MTT_S3, and to all other tracer data for the estimation of MTT_S3.

Model S4 is used solely for the calculation of peripheral MTT_S4 of the tracers. It includes injection a little after left ventricle (at α₁ of Fig. 1), sampling a little before right atrium (at γ of Fig. 1), and attenuated elimination along the extended region β_ext (Fig. 1). This model was applied to all tracer data.

Model S5 is used solely for the calculation of pulmonary mean transit time of the tracer excluding the heart, MTT_S5. It includes injection a little after right ventricle (at δ of Fig. 1), sampling a little before left atrium (at ζ of Fig. 1), and attenuated elimination at region ϵ. This model was applied to all tracer data.

The algorithm used for the fittings is a classical least-squares optimizer (LMDIF subroutine of the MINPACK Fortran optimization package; http://www.netlib.org/minpack/). The PDE is solved numerically using a Crank-Nicolson implicit finite-difference scheme (Burden and Faires, 1993), which reduces the solution of a PDE to the solution of a system of linear algebraic equations.

Results

Initially, eq. 4 was fitted to all eight data sets and estimates for the seven parameters were derived. The fittings for human and dog data resulted in similar values for the parameters x_c and x_p for each species. The values are x_c = 0.40 and x_p = 0.90 for the human, and x_c = 0.15 and x_p = 0.65 for the dog data. These sets of values were fixed for the later use of the model, thus eliminating two parameters and leaving the model with five remaining parameters for estimation.

The first step in our analysis was the fitting of the appropriate model F1, F2, or F3 of Table 1 to the eight sets of experimental data. A representative sample of fittings is shown in Fig. 2A for the data set I with the model F1. Visual inspection of Fig. 2A reveals that the model describes the experimental data nicely. The values of the five parameters along with the relevant standard errors derived from all fittings are quoted in Table 2. The parameter values of Table 2 that were derived from the optimization were used subsequently to generate C,t data without the effect of recirculation using the models S1 to S5 listed in Table 1. A typical plot for data set I where recirculation was deliberately suppressed is shown in Fig. 2B.

Simulations with models S1 to S5 using the model parameters that correspond to data set I and various injection and sampling scenarios are shown in Fig. 3. The exercise was applied for all data sets and MTT estimates were calculated using eq. 2. The values for MTT, derived for data sets I to VIII, are listed in Table 3. Also, in the case of ICG administration (data sets I, II, III, and VI), apart from the MTT, the cardiac outputs CO_S1, CO_S2, and CO_S3, as well as the respective volumes (V_S1, V_S2, and V_S3) of ICG distribution were estimated, using models S1, S2, and S3, respectively. The results of these estimations are listed in Table 4.

Discussion

The values of the estimated model parameters (Table 2) derived from the fitting of eq. 2 to the data sets I to VIII are considered physiologically reasonable. Although the velocity (U) values are normalized, they can be transformed to classical units (centimeters per second) by multiplying them with a length factor. This factor is estimated from the volume V and assuming that the cross-sectional area of the aorta is 4.5 cm² (Labarbera, 1990). By doing so, the normalized velocity value U = 0.0298 s^-1 is transformed to 0.30 cm/s, which is a reasonable finding. Other remarks can be also pointed out for the velocity values. The velocity of data set II is smaller than the one of data set I, a fact that can be attributed to the presence of propranolol in the experiment of data set II (Niemann et al., 2000). The comparison of the values of the uptake rate constant, K, cannot be made easily because this quantity has different physical meaning for each tracer and is also dependent on the length of the uptake region assigned for each tracer. Also, difficult to assess are the absolute values of the dispersion coefficients D₁ and D₂; however, comparative remarks are possible, e.g., D₂, which corresponds to the single pulmonary tree, is always smaller than D₁, which corresponds to the multiple peripheral trees. The partial volumes, V, correspond to the volumes of distribution of the tracers during the initial mixing phase. However, the absolute values of the volumes of data sets VI, VII, and VIII, have no physical meaning, because these data sets were published in normalized units and not in concentration units (Krejcie et al., 1996).

The C,t plot presented in Fig. 2A indicate that when recirculation is present, apart from the presence of the second peak of the plot, the lower part of the declining phase of the first peak is attenuated. This is not the case when recirculation is not present. An indicative C-t curve from data set I, but without the recirculation is shown in Fig. 2B exhibiting that recirculation has been suppressed completely.

Figure 3 shows the C-t curves derived after simulations with models S1 to S5 in case of certain injection and sampling scenarios. Visual inspection of these plots reveals that the results are reasonable. For example, the shape of the curves that correspond to longer blood route, are wider and less sharp (Fig. 3D, model S4) than the curves that correspond to an injection/sampling scenario of shorter distance (Fig. 3C, model S3).

The values of the various MTT estimates, derived from the application of the models S1 to S5 to data sets I to VIII, are quoted in Table 3. It should be noted that in the case of data sets VI, VII, and VIII, the concentration values reported in the relevant article are provided in a normalized scale (Krejcie et al., 1996). This is reflected in the estimates of CO and V but not in the estimates of MTT because the latter are not affected by the normalized concentration values. To transform the normalized values of data set VI to concentration units, we used a scale factor derived from data set III because data sets III and VI have identical experimental design, whereas data set III is available both in classical and normalized units (Krejcie et al., 1997). Data sets VII and VIII are only used for the MTT estimation, so transformation was not necessary.

By inspecting the values of MTT in Table 3, interesting observations can be made. The peripheral MTT_S4 is always greater than the pulmonary MTT_S5, a physiologically sound finding that is more attenuated in humans than in dogs. This difference seems to be dependent on the size of the species and is probably due to the higher inhomogeneity of the tracer distribution in the peripheral vessel trees compared with pulmonary trees. Due to the presence of propranolol, the MTT values of data set II are larger than the respective values of data set I. In principle, the values of MTT depend on whether the tracer is intravascular or extravascular. The more a tracer is distributed extravascularly, the higher the MTT value is. Consequently, comparing data sets of the same subjects (data sets III versus V and VI versus VIII), we can observe that the MTT values for data set VI (antipyrine) are higher than the values of data set III (ICG), and the values of data set V (lidocaine) are even higher. Accordingly, the MTT values for data set VI (ICG) are equivalent to data set VIII values (inulin), but lower than those of data set VII (antipyrine).

Reasonable findings were also obtained for the cardiac output and the partial volumes (Table 4). Data set II exhibits a significantly lower value for CO than data set I, as it is expected due to the administration of propranolol. Besides, the calculated absolute values of CO, using the proposed methodology, were similar to those derived from the thermodilution technique (Krejcie et al., 1996, 1997; Niemann et al., 2000). This can be shown concisely if one compares the estimates 8.4, 4.1, 2.52, and 2.51 l/min for CO_S1 listed in Table 4 with the reported values 10.6 ± 2.5, 4.1 ± 0.7, 3.22 ± 0.81, and 2.49 ± 0.93 l/min, respectively (Krejcie et al., 1996, 1997; Niemann et al., 2000). Rational estimates were also obtained for the partial blood volumes. The proportion of blood volume being in the pulmonary circulation (heart plus pulmonary tree) was found to be 20 and 23% for data set I and II, respectively, in contrast to the approximate value of 16% reported in literature (Guyton and Hall, 2000). It should be noted that the estimates for MTT and partial blood volumes derived from the analysis of the model variants S2 and S3 are heavily dependent on the specific values used for x_inj and x_sam because of the short distance between the injection and sampling sites. However, this limitation does not affect the AUC estimates and the corresponding CO values.

Overall, reasonable estimates for the physiological variables can be obtained using the present method, when appropriate indicator dilution data are available. As expected, the most appropriate data used in our analysis are the ICG data that correspond to strictly intravascular distribution. Thus, the strength and the flexibility of the approach was demonstrated when valid estimates for the mean transit time, the cardiac output and the respective volume were derived from the analysis of ICG data. When tracers with extravascular distribution were used, the estimations were limited to the mean transit time.

Conclusions

In this study, we presented a novel approach to analyze indicator dilution data to estimate various recirculation parameters. The procedure consists of two steps. In a first step, the mathematical model is fitted to the C,t data, to estimate the various parameters of the model. Then, in a second step, a simulation is performed using the model with the set of parameters estimated from the first step, but suppressing the recirculation, i.e., conducting a virtual indicator dilution experiment. In this way, the derived area under the curve gives the opportunity to estimate various recirculation parameters. Moreover, the model allows the freedom to set the administration and sampling sites of the model to appropriate positions depending on desired recirculatory parameters.

The procedure was applied successfully to indicator dilution data of various tracers from humans and dogs. The intravascular tracer indocyanine green gave valid estimates for CO, MTT, and partial volumes; however, the meaningful estimates acquired using other tracers that distribute outside the vascular bed, are limited to the MTT. The results presented show that the method developed works well as an analysis tool for indicator dilution data. Using this method the recirculation effect of a tracer can be suppressed, which is the main task of all the alternative previously published techniques. However, here this is done not in an arbitrary manner but in physiological basis. Moreover, entire virtual experiments may be implemented that include different administration and sampling sites, that suit the needs of the desired measurement, offering in this way additional flexibility. The method, although it includes a complicated mathematical formulation, it can be regarded as relatively easy to apply, because the model can be tuned to be applied as a black box without detailed knowledge of the mathematics or the algorithms used for the solution of the PDE.