Optimized dynamic framing for PET-based myocardial blood flow estimation

We specify the input function in terms of a gamma variate function, which is capable of accurately reproducing TACs in the cavity of the left ventricle (LV) following a bolus injection. For convenience in estimating the input function model parameters, we describe the gamma variate Γ(t) in terms of parameters (C_max, t_max) corresponding to the peak's amplitude and time, α for shape and t₀ for time delay (Madsen 1992).

$$\begin{equation} {\rm \Gamma }( {\rm t} ) = {\rm C}_{{\rm max}} {\rm *}\left( {\frac{{t - t_0 }}{{t_{{\rm max}} - t_0 }}} \right)^\alpha {\rm exp}\left( {\alpha {\rm *}\left( {1 - \frac{{t - t_0 }}{{t_{{\rm max}} - t_0 }}} \right)} \right)\;{\rm for}\;t > t_0 \end{equation} \tag{ 1 }$$

and where Γ(t) = 0 for t ⩽ t₀. Whereas Γ(t) models in part the arterial activity concentration time course following a perfect bolus injection, we are interested in describing the arterial activity following ⁸²Rb administered using an ⁸²Sr/⁸²Rb generator infusion system (Bracco Diagnostics Inc., Princeton, NJ), which typically delivers rubidium chloride over 5 to 25 s. To account for this, we convolve Γ(t) with a rectangular infusion pulse with unit area and duration t_Inf. We also add a rising asymptotic term, which represents a combination of equilibrium activity in the blood and spillover from the myocardial uptake due to partial volume and limited spatial resolution (Golish et al 2001), in order to obtain input function C_a(t) representing the time course of radioactivity concentration in arterial blood

$$\begin{equation} C_a ( t ) = \left[ {{\rm rect}\left( {\frac{t}{{t_{{\rm Inf}} }}} \right) \otimes {\rm \Gamma }( {\rm t} )} \right] + C_0 \left[ {1 - {\rm exp}\left( { - \frac{{( {t - t_0 } )}}{\beta }} \right)} \right]\end{equation} \tag{ 2 }$$

where rect(t/t_Inf) equals 1/t_Inf on the interval t ∈ [0, t_inf] and 0 otherwise, C₀ is the steady-state value of the measured blood activity, and C_a(t) = 0 for t ⩽ t₀. In clinical practice, we estimate C_a(t) from a TAC for a volume drawn in the LV cavity. The model-predicted PET measurement for arterial radioactivity concentration in the jth frame M_{Ca, j} is determined as the value of C_a(t) averaged over the duration of frame j, which spans the time interval from t_j to t_j+1.

$$\begin{equation} M_{Ca,j} = \frac{1}{{t_{j + 1} - t_j }}\mathop \int_{t_j }^{t_{j + 1} } C_a ( t )\,{\rm d}t\end{equation} \tag{ 3 }$$

To achieve clinically realistic input functions, we determined the values of parameters C_max, t_max, C₀, α, β, t₀ and t_inf by fitting (3) to the data obtained from five human, resting ⁸²Rb scans. In these scans, PET data were acquired in list-mode on a Gemini TF-16 PET/CT (Philips Healthcare, Cleveland, OH, USA) (Surti et al 2007). The acquisitions began concurrently with infusion of 307 ± 60 MBq of ⁸²Rb in 6 ± 2 mL of saline. At the peak count rate of each scan, we estimate that no more than 30% of the events were lost to system dead time and that with dead time correction the bias in the values was less than 5% of the true value.

Starting at the time at which activity arrived in the field of view, as judged by the singles rate exceeding 5% of background, PET data were binned into a dynamic time series with 20 frames each of 1 s duration (abbreviated 20  ×  1 s), followed by 8  ×  5 s, 17  ×  15 s, and 2  ×  30 s. The frames were individually reconstructed using a blob-based, list-mode OSEM algorithm (Popescu et al 2004) with correction for attenuation, scatter, randoms, dead time and radioactive decay. TACs for the LV cavity were generated from the data by taking the mean activity concentration over a 0.22 mL volume of interest drawn at the base of the heart in a short-axis view. To fit the TACs using (3), first values for C_max, t_max and C₀ were calculated from a fit to the peak and tail of the measured data, and then α, β, t₀ and t_inf were iteratively estimated using lsqcurvefit in MATLAB (Optimization Toolbox 6.2.1; The Mathworks, Inc., Natick, MA, USA). The means of estimated parameter values were used to construct the input function shape for subsequent simulation of tissue TACs; the variety of LV cavity TAC shapes among clinical patients was therefore not simulated.

To create simulated data for MBF estimation, we used the one-tissue-compartment physiologic model described by Hutchins et al (1990), who used ammonia as the radiotracer, and later applied to ⁸²Rb by Coxson et al (1995). With noisy data, this model is a good choice as it facilitates a more reliable fit than multiple-compartment models (Coxson et al 1997). The model equation representing the activity concentration contributing to the myocardial PET signal is

$$\begin{equation} C_T ( t ) = f_v C_a ( t ) + ( {1 - f_v } )\mathop \int _0^t {\rm e}^{\{ { - k_2 ( {t - t'} )} \}} K_1 C_a ( {t'} )\,{\rm d}t'\end{equation} \tag{ 4 }$$

where C_a(t) is the activity concentration in the blood and the convolution term represents the accumulated radioactivity concentration in the myocardium. The measured tissue activity was modelled as a function of parameters for uptake K₁, washout k₂, and a third parameter f_v that accounts for both the vascular fraction and spillover from the LV cavity to the myocardium tissue region. The model-predicted activity in the jth contiguous image frame is calculated as the activity concentration time-averaged over the frame interval,

$$\begin{equation} M_j = \frac{1}{{t_{j + 1} - t_j }}\mathop \int _{t_j }^{t_{j + 1} } C_T ( t ){\rm d}t.\end{equation} \tag{ 5 }$$

The integrals in (3)–(5) were numerically evaluated using CVODES (Serban and Hindmarsh 2005) as implemented in COMKAT (Muzic and Cornelius 2001, Fang et al 2010). Noise Σ_j was then added to the noise-free, model-predicted tissue concentration M_j to generate a simulated data value D_j. The noise in each frame was calculated by sampling from the standard normal distribution and scaling to achieve the specified standard deviation (SD; σ).

$$\begin{equation} D_j = M_j + {\rm \Sigma }_j ;\;{\rm \Sigma }_j \sim N\big( {0,\sigma _j^2 } \big). \end{equation} \tag{ 6 }$$

In our noise model, variance $\sigma _j^2$ was proportional to the number of counts contributing to that frame, which is calculated from the frame duration and a decay term (t_H = 1.273 min) as the PET data were decay-corrected during image reconstruction.

$$\begin{equation} \sigma _j^2 = \xi ^2 *M_j *( {t_{j + 1} - t_j } )^{ - 1} *2^{(t_j /t_H )} . \end{equation} \tag{ 7 }$$

The value of proportionality constant ξ was determined such that the average SD matches a value specified (e.g. 10% of M_j) for each set of simulated TAC.

To determine what noise level appropriately models that of clinical imaging, we scanned a torso phantom filled with 7.2 kBq mL⁻¹ of ¹⁸F containing a cardiac insert filled with 100 kBq mL⁻¹ (Data Spectrum, Hillsborough, NC, USA). Statistical replicate images, each with approximately the same number of counts as one frame of a human study, were created by splitting list-mode data into 1 ms intervals and dealing successive intervals into 20 data files to create 20 independent datasets with, on average, the desired number of counts. Each replicate data file was reconstructed, a volume of interest was drawn in the myocardium insert, and the SD among the replicate images of the volume's mean was calculated. Simulations were then focused on these clinically relevant noise levels.

Values for the three parameters, K₁, k₂, and f_v, of the physiologic model were estimated by minimizing a weighted least-squares (WLS) objective function:

$$\begin{equation} O_{{\rm WLS}} ( {K_1 ,k_2 ,f_v } ) = \frac{1}{2}\mathop \sum \limits_{j = 1}^{nf} w_j ( {M_j ( {K_1 ,k_2 ,f_v } ) - D_j } )^2 \end{equation} \tag{ 8 }$$

where M_j and D_j are the values of the model output and measured (or simulated) data at frame j of the number nf of total frames and w_j are weights, specified as reciprocals of the variance (Muzic and Christian 2006) and estimated from the number of counts in the frame, i.e. ξ²/σ² in (7). For experiments in which parameter estimation is performed infrequently and the quality of the estimate can be verified by careful inspection, the WLS objective function can be minimized using nonlinear optimization. A more automated, efficient and robust method is desired for routine clinical use. Accordingly, we implemented a means to minimize the WLS objective function that takes advantage of certain properties of the model formulation. In particular, in the PET model equation, i.e. combining (4) and (5),

$$\begin{equation} M_j = \frac{1}{{t_{j + 1} - t_j }}\mathop \int _{t_j }^{t_{j + 1} } f_v C_a ( t)\,{\rm d}t + \frac{1}{{t_{j + 1} - t_j }}\mathop \int _{t_j }^{t_{j + 1} } \left[ {( {1 - f_v } )\mathop \int_0^t {\rm e}^{\{ { - k_2 ( {t - t'} )} \}} K_1 C_a ( {t'} )\,{\rm d}t'\,{\rm d}t} \right]\end{equation} \tag{ 9 }$$

K₁ can be factored outside the integral and it becomes evident that the equation is nonlinear in only the k₂ parameter. In considering the values of functions time-averaged over each frame, we can write this as a vector-matrix operation,

$$\begin{equation} M = \left[ {\begin{array}{@{}c@{}} {M_1 } \\ {M_2 } \\ \vdots \\ {M_{nf} } \\ \end{array}} \right] = X\beta \end{equation} \tag{ 10 }$$

with

$$\begin{eqnarray} X = \left[ {\begin{array}{@{}cc@{}} M_{Ca,1} & \frac{1}{{t_2 - t_1 }}\int_{t_1 }^{t_2 } {\left( {\int_{t' = t_0 }^t {C_a (t')\,{\rm e}^{\{ - k_2 (t - t')\}}\, {\rm d}t'} } \right)\,{\rm d}t} \\ M_{Ca,2} & \frac{1}{{t_3 - t_2 }}\int_{t_2 }^{t_3 } {\left( {\int_{t' = t_0 }^t {C_a (t')\,{\rm e}^{\{ - k_2 (t - t')\}}\, {\rm d}t'} } \right)\,{\rm d}t} \\ \vdots &\vdots\\ M_{Ca,nf} & \frac{1}{{t_{nf + 1} - t_{nf} }}\int_{t_{nf} }^{t_{nf + 1} } {\left( {\int_{t' = t_0 }^t {C_a (t')\,{\rm e}^{\{ - k_2 (t - t')\}} \,{\rm d}t'} } \right)\,{\rm d}t} \\ \end{array}} \right],\quad \beta = \left[ {\begin{array}{@{}c@{}} f_v \\ ( {1 - f_v } )K_1 \end{array}} \right].\nonumber\\ \end{eqnarray} \tag{ 11 }$$

Given a specified arterial input function and value of k₂, β was determined from the measured data by solving, in a least-squares sense,

$$\begin{equation} WM = WX\beta \end{equation} \tag{ 12 }$$

where W is a diagonal matrix holding the weights w_j for each frame. Indeed, the value of β that minimizes O_WLS is available in closed form as

$$\begin{equation} \beta = ( {WX} )^{ - 1} WM.\end{equation} \tag{ 13 }$$

Therefore, for a given value of k₂, f_v = β₁ and K₁ = β₂/(1 − β₁) could be determined. This calculation was repeated for a list, i.e. a one-dimensional grid, of k₂ values for each of which the value of O_WLS was calculated. The value of k₂ that gave the lowest value of the WLS objective function, as well as its corresponding K₁ and f_v, were taken as the optimal estimate of the parameter vector. This method of estimating the parameter values is similar to the linearized least squares approach (Koeppe et al 1985), which was extended to include a spillover term in studying blood flow in tumours (Lodge et al 2000). To simulate the use of an image-derived input function, we modelled the PET measurement of the input function for each framing schedule using (3), and then applied linear interpolation of the frame values as the 'measured' input function for parameter estimation.

The goal of optimal experiment design is to maximize the sensitivity of the data to the parameters of interest in order to maximize the precision of the parameter estimates. Various metrics have been applied to quantify this sensitivity. We used a common one called the D-optimal criterion which entails maximizing the determinant of the Hessian (H), the matrix of second derivatives of the objective function with respect to the parameter values (Muzic et al 1996). The Hessian provides a measure of the steepness with which the objective function increases as one moves away from its minimum value with greater steepness yielding better precision of the parameter estimates. Maximizing the determinant of the Hessian corresponds to minimizing the volume of the confidence region. As the Hessian depends on the values of the parameters, to evaluate it we assumed representative values for the physiologic parameters, i.e. [K₁, k₂ and f_v] = [1.0 mL g⁻¹ min⁻¹, 0.3 min⁻¹, 0.3 (unitless)], and adjusted the values of the experiment design parameters, such as the infusion duration or framing schedule, as to maximize the determinant of the Hessian of the WLS objective function. We then tested such designs using Monte Carlo simulations wherein K₁, k₂ and f_v values are estimated for a large set of simulated TACs.

The Hessian of the WLS objective function can be approximated in terms of the model sensitivity functions and weighting matrix W (Bard 1974).

$$\begin{equation} J_{j,i} = \frac{{\partial M_j }}{{\partial \theta _i }}{\rm }\end{equation} \tag{ 14 }$$

$$\begin{equation} {\rm det}( H ) = {\rm det}( {J^T WJ} )\end{equation} \tag{ 15 }$$

where θ = [K₁, k₂, f_v]' is the parameter vector. We maximized det(H) using two approaches. In the first, the schedule was parameterized as a list of independent frame durations with no gaps or overlaps between frames, and the optimal schedule was determined by minimizing $1/\det ( {\rm H} )$ using an active-set algorithm available in MATLAB (fmincon, Optimization Toolbox 6.2.1; The Mathworks, Inc., Natick, MA, USA). We observed that when the sensitivity functions change only slightly between frames, the optimization methods had difficulty. As the results were sensitive to initial guesses and converged to local minima, we tried a second approach, an exhaustive search in which the frame durations were parameterized in clusters and also constrained to integer numbers of seconds. Three different frame durations were allowed in a framing design we denote as (n × X) + (m × Y) + (p × Z) + (1 × R), where n, m and p are the numbers of frames of respective durations X, Y, and Z seconds. A final frame of the remaining duration R was included to achieve a fixed, 6 min acquisition length. Naturally, n, m and p must be integers, and we restricted them to be values from 0 to 20, based on a cursory analysis showing diminishing returns with more than 60 total frames. While frame durations could be real numbers, in practice we considered integer values from 0 to 30 and also constrained Z ⩾ Y ⩾ X. Given the discrete problem, we calculated det(H) for each possible schedule in the search space. To gauge the potential of our approach for improving current practice, we evaluated published experimental protocols for MBF estimation (Lortie et al 2007, El Fakhri et al 2009, Raylman et al 1993, Sherif et al 2011). In the published schedules, we extended the final frame or removed frames as needed in order to match a fixed 6 min total acquisition duration. We then calculated det(H) for the published sampling schedules and for a schedule with the same number of frames, that had we optimized.

In the previous section, the framing schedule was optimized for a particular input function and assumed values of blood flow and other model parameters. In contrast, a clinical population represents a range of parameter values and an experiment design that works well in patients with normal physiologic function might work poorly in pathologic conditions. We addressed this by expanding the approach described in the previous section by considering four sets of physiologic parameter values denoted by A, B, C and D in table 1. Simulated subjects A and B were based on fits to stress and rest data obtained from healthy volunteers (Lortie et al 2007). Simulated subject C was based on fits to a canine model of ischemia under pharmacologic stress (Herrero et al 1992), and subject D was based on fits to the rest data described in section 2.1. The different spillover fractions likely reflect the variation in heart size, scanner resolution, motion, and image reconstruction methods. A 10% noise SD was used, as calculated from the phantom study described in section 2.2.

Using a minimax optimization, robustness is achieved by determining the schedule where the minimum objective function value is maximized—that is, minimax optimization identifies the schedule with the best worst-case performance, in our case according to the D-optimal criterion (Salinas et al 2007). As in the previous section, we calculated det(H) using an exhaustive search: for each framing schedule in the search space, we computed the determinant of the Hessian for each simulated subject. For each schedule, we identified the minimum value of det(H) using that schedule among the subjects; the schedule with the highest minimum det(H), or equivalently the lowest, maximum value of the objective function, was the most robust.

For ⁸²Rb-PET, a slow infusion has been proposed as an alternative to a bolus-like infusion, in order to minimize both the scanner's dead time (deKemp et al 2008) and the temporal sampling requirements needed to capture the rapid dynamics (Raylman et al 1993). We modelled the duration of the tracer infusion by changing t_Inf in (2), covering a range of durations suggested in the literature, from a bolus (El Fakhri et al 2009) to a 4 min infusion (deKemp et al 2006) while holding the total injected activity constant. The D-optimality criterion was calculated for each duration, and Monte Carlo simulations were performed for selected durations to measure the precision of the parameter estimates. In this way we quantitatively evaluated the anticipated benefit of a rapid infusion. Infusion duration was studied at a fixed, 6 min acquisition.

We also studied the effect of modifying the duration of PET data acquisition in order to select a study duration that achieves high precision without excessively extending the scan duration into a region of diminishing returns. Using a fixed, 10 s infusion, the D-optimality criterion was calculated for data durations ranging from 0.5 to 20 min, using an optimized framing schedule for each duration. Selected durations were examined using Monte Carlo simulation, and the bias and precision of parameter estimates were calculated.

The input function model was applied to TACs from the LV cavity seen in five human scans in order to obtain a typical curve shape. Table S1 (available from stacks.iop.org/PMB/58/5783/mmedia) shows the range of parameters describing the curves derived from the five images; figure S1 (available from stacks.iop.org/PMB/58/5783/mmedia) shows the model fit for the first subject. The normalized root-mean-square of the residual differences between the fits and data for the five subjects ranged from 6%–14%, suggesting that the input function model approximates the data well. The infusion durations estimated from the input function, i.e. 0.11–0.38 min, fell in the range typically reported by the infusion system. Based on the quality of these fits, we concluded that the mean curve, labelled 'Model' in table S1 (available from stacks.iop.org/PMB/58/5783/mmedia), was a plausible representation of a human measurement, and we used it as the input function for subsequent simulations. For simulating arterial input functions, we fixed the delay parameter at 3 s and the infusion length at zero to initially study a bolus-like infusion.

The model sensitivity functions J_{i, j} and the D-optimal criterion $\det ( H )$ were first calculated for nominal values for the physiological parameter values, [K₁, k₂ and f_v] = [1.0 mL g⁻¹ min⁻¹, 0.3 min⁻¹, 0.3 (unitless)]. As expected, the sensitivity functions had the greatest amplitudes at early times and had distinct shapes, thus suggesting linear independence hence good local identifiability. The sensitivities are shown together and scaled to facilitate visual comparison of their shapes in figure 1(a). The D-optimal criterion was calculated from the sensitivity functions, for sampling schedules with uniform frame durations ranging from 0.5 to 15 s (720 to 24 equally-spaced frames), figure 1(b). As expected, more frequent sampling leads to an increased value of det(H) with very small improvements beyond approximately 180 frames, which corresponds to frame durations of 2 s.

Zoom In Zoom Out Reset image size

Noise-free and noisy TACs were then simulated with uniformly spaced frames of 5 s duration. Noise was added according to (6) with a specified average SD ranging from 1% to 50%. Estimates for the physiologic model parameters were calculated as described in section 2.3 using a uniformly-spaced grid of 200 k₂ values ranging from 0 to 1 min⁻¹. Parameter estimation for 12 noise levels, each with 30,000 TACs, took 14 min on a standard desktop PC (running on four cores at 2.80 GHz). Using this number of TACs and 0.5% resolution in k₂ was deemed to be more than sufficient to estimate the mean and SD of parameter estimates on the basis of cursory studies comparing 10⁶ and progressively fewer TACs at high and low noise levels.

As shown in figure 2, the parameter estimates had low biases, i.e. bars are centred near zero error, and variances, i.e. the widths of the bars, decrease with noise level. The bias of K₁, was less than 1% of the true value with noise levels up to 20% SD. The bias at the highest noise level (50% SD) was (3.4 ± 0.4)% of the true value (99% CI). This indicated that our noisy simulations and parameter estimates behaved as expected and that there was no evidence of programming errors.

Zoom In Zoom Out Reset image size

To determine an appropriate noise level for subsequent simulations, we analysed the cardiac phantom data described in section 2.2. Volumes of interest were drawn in the cardiac insert portion of the phantom, corresponding to single image voxels (0.01 mL), and small, average and large segments (1.7, 3.5 and 5.0 mL) in a standard 17-segment cardiac model (Cerqueira 2002), as well as the entire LV myocardium (64 mL). As expected, image noise decreased smoothly with an increase in the number of coincidences or region size, figure S2 (available from stacks.iop.org/PMB/58/5783/mmedia). For subsequent simulations, we chose a noise level ξ = 0.1 which, when applied following (7), gave simulated TACs with a frame SD in the range of 4% (for a 15 s frame at 120 s) to 15% (for a 1 s early frame or a 15 s frame at 6 min).

In the first optimization method, the schedule was parameterized as a list of 13 contiguous frame durations. The quantity $\left[ {\det ( H )} \right]^{ - 1}$ was used as the objective function and was minimized using MATLAB's fmincon with an initial guess of uniform (18 s) frames (figure 3(a)) and constrained such that the sum of all frames equalled the 6 min study duration. As expected, the iterative optimization preferred numerous short-duration frames early after the arrival of activity and longer duration frames later, when the curve shape was not changing as rapidly (figure 3(b)). Iterative optimization was feasible with up to 13 frames with careful attention given to initial guesses, tolerances, and other optimization algorithm parameters in order to avoid local minima. Next, we exhaustively searched 13-frame designs following the parameterization method described in section 2.4. The optimal 13-frame schedule from the exhaustive search is shown in figure 3(c); its corresponding det(H) 8.71×10¹², slightly lower than the iteratively optimized 13-frame schedule (det(H) = 9.61×10¹²) shown in figure 3(b). The slightly higher value with iterative optimization is due to the iterative method having access to parameter values between those on the exhaustive method's grid. As expected, both optimization approaches yielded higher values for the objective function, and are therefore expected to provide better estimate precision, than the initial condition (figure 3(a)) with uniform frame durations (det(H) = 2.95×10¹²).

Zoom In Zoom Out Reset image size

The experiment design in figure 3(c) is one of many considered in the exhaustive search process. Figure 4 shows the det(H) calculated for the range of the exhaustive search. Among each column of points, which represents a fixed number of frames, the top-most point has the highest det(H) and is the optimal design for that number of frames. We found that for any number of frames, the optimal schedule gave a higher det(H) than any optimal schedule with fewer frames. The more-optimal schedules tend to have several short frames clustered early, when the signal is changing rapidly and near the sensitivity functions' maxima, and longer frames at later times, following the overall pattern seen in figures 3(b) and (c). As examples, the optimal schedules for 5, 10, 20 and 50 total frames are noted in figure 4 by stars. Those schedules are

We then compared det(H) for our optimized schedules to experiment schedules described in literature. The schedules described in literature are described in table 2 and indicated in figure 4 by triangles. In all cases our design optimization was able to determine sampling schedules that improved the previously published schedules.

Zoom In Zoom Out Reset image size

To investigate the correspondence of det(H) and the precision of the flow estimate, MC simulation was performed for the a subset of the most promising schedules. Figure 5(a) shows that the optimal five-frame schedule performed poorly, thus resulting in a wide and non-normal distribution of parameter estimates, and that the optimal 10, 20 and 50 frame schedules performed progressively better, with estimates clustered closer to the true parameter values with a larger number of frames. We expanded the analysis using a larger number of schedules, selected the 44 schedules with the highest det(H) for each total number of frames from 17 to 60, and constructed a scatter plot of the precision of the K₁ estimate versus det(H) (figure 5(b)). This demonstrated that the standard deviation of K₁ estimates decreased with increasing det(H) and thus provided evidence that det(H) is a valid predictor of estimate precision and appropriate for use in experiment design.

Zoom In Zoom Out Reset image size

We calculated det(H) and performed Monte Carlo simulations for a range of schedules, infusion lengths and physiologic parameters, with the goal of determining a single schedule applicable to the patient population. We examined the 30-frame schedules in our exhaustive search parameterization, having observed that much of the available det(H) in the entire grid is reached with 30 frames (figure 4) and that 30-frame images are manageable in research and clinical use. Figure 6 shows schedules, ordered in descending values of det(H) for subject C since this had the poorest det(H) therefore is the controlling subject for the robust design. The first schedule thus is the minimax or optimum for the robust design and consists of 29 1 s frames followed by a single, 331 s frame (29×1 + 1×331). In comparing two schedules, generally one which performed poorer for subject C also performed poorer for the other subjects although there are exceptions to this as indicated by jumps in the det(H) curve such that it is not monotonically decreasing. However, focusing on the first 1000 schedules, figure 6(b), reveals that the rank-ordering of the best-performing schedules shows similarities across subjects. That is, many of the best schedules for subject C were also among the best for A, B, and D.

Zoom In Zoom Out Reset image size

Up to this point we considered an infusion of minimal duration (t_inf = 0) because we expected short infusions to require the most demanding framing schedule and because requiring a different framing schedule for each infusion duration might not be practical for clinical use. We explored the validity of this approach by evaluating the effect of infusion length, determined by calculating det(H) for a range of values of t_Inf in (2). The D-optimal criterion was first calculated using fine, i.e. uniform 1 s, frame sampling to ensure that there was sufficient time resolution in order to capture the dynamics and decouple the problem of optimizing the infusion length from the problem of optimizing the framing schedule. As expected, the optimality was highest for the shortest infusion duration (figure 7), although infusions shorter than 3 s offered only slight additional benefit, i.e. an infusion of 3 s achieves 98% of the maximum det(H) of a 0.25 s infusion. This result suggests that a short, 3 s infusion is nearly as effective as a bolus injection. Next, the det(H) calculation was repeated with 1, 5, 15, 30 and 60 s infusion lengths, with a range of schedules in order to determine a robust schedule for use with different infusion durations. Subject A's physiologic parameters were used. For the 1 s infusion, a schedule with rapid initial sampling was optimal, i.e. (20×1 + 20×2 + 20×10 + 1×100). Schedules with slightly longer initial frames were optimal for slower infusions; for the 60 s infusion, the optimal schedule was (40×3 + 20×10 + 1×40). As with the physiologic parameter results, the rank order of schedules for each infusion duration was similar, and many similar schedules gave similar det(H) values for a given infusion duration. As shown in figure 7, the five schedules that are optimum for each of the five infusion durations are nearly optimal when applied to data simulated with each of the other infusion durations. Applying the optimal schedule for a 60 s duration to all infusion durations, in order to maximize the minimum det(H), still gives at worst 98.5% of the optimal det(H) and an increase in K₁ SD from 3.5% to 3.6% of the its true value. Bias was examined using simulation and was not significantly affected by the infusion duration; this is consistent with our choice of 1 s time sampling for this part of the study.

Zoom In Zoom Out Reset image size

The preceding simulations were performed using a fixed, 6 min acquisition duration. We studied the effect of acquisition duration by simulating acquisition durations ranging from 30 s to 20 min, using 10% noise, uniform 1 s frame durations, and a clinically realistic, 10 s infusion duration,. Acquisitions of 4 and 6 min' duration achieved 95% and 99%, respectively, of the maximum det(H) of the longest (20 min) acquisition. The standard deviation of the estimated K₁ and k₂ values from simulated data followed a similar pattern, with little change between the 6 min and 20 min acquisitions. Bias was detected in K₁ and k₂, although it remained less than 1.5% of the true parameter values.