|Home | About | Journals | Submit | Contact Us | Français|
We present a direct method for producing images of kinetic parameters from list mode PET data. The time-activity curve for each voxel is described by a one-tissue compartment, 2-parameter model. Extending previous EM algorithms, a new spatiotemporal complete data space was introduced to optimize the maximum likelihood function. This leads to a straightforward parametric image update equation with moderate additional computation requirements compared to the conventional algorithm. Qualitative and quantitative evaluations were performed using 2D (x,t) and 4D (x,y,z,t) simulated list mode data for a brain receptor study. Comparisons with the two-step approach (frame-based reconstruction followed by voxel-by-voxel parameter estimation) show that the proposed method can lead to accurate estimation of the parametric image values with reduced variance, especially for the volume of distribution (VT).
Dynamic positron emission tomography (PET) permits quantification of tracer dynamics. The current data processing path for parametric imaging is to reconstruct a time series of images from measured projection data independently and then estimates each voxel’s kinetic parameters from the time-activity curve (TAC), typically with a compartmental model. This frame-based approach requires selection of the duration of each frame, involving a choice between collecting longer scans with good counting statistics but poor temporal resolution, or shorter scans that are noisy but preserve temporal resolution. Optimal use of the dynamic data requires accurate noise estimates for data weighting; this estimation is challenging for nonlinear iterative reconstruction methods because its noise is spatially variant and object dependent.
Direct approaches for creation of parametric images have been in the literature for over 20 years. In 1983, Snyder  developed a list mode expectation-maximization (EM) maximum likelihood (ML) algorithm for estimation of parametric images using inhomogeneous spatial-temporal Poisson processes and a kinetic compartmental model. Carson and Lange  also proposed an EM framework for direct parametric image reconstruction algorithm with a one tissue model. Subsequently, many direct kinetic estimation methods [3-8] for sinogram data have been produced. However, for a high-resolution scanner such as the HRRT with 4.5×109 potential lines of response, list mode data is preferred; this can reduce the data storage requirements while maintaining highest resolution by storing the measured attributes of each event in the list. Other direct methods [9-14] were developed from linear basis functions, whose the temporal model is linear with respect to the parameters. Although such basis functions can represent the TACs, they have no direct physiological meaning and kinetic parameters must still be calculated from the dynamic images, leading again to a two-step process. The application of kinetic compartment models is more biologically based, becoming the primary goal of dynamic PET. Tsoumpas et al  and Wang et al  proposed direct parametric image reconstruction methods, whose compartment model is linear (Patlak plot) model, and is thus limited to irreversible tracers such as FDG. Here, we present a new EM-based direct parametric image reconstruction algorithm for list mode data by combining nonlinear one tissue compartment model (1T) into the physical model.
We extend the notation of Lange and Carson  from static to 4D. The physical model for the projection measurement Yit on line of response (LOR) i at time t is:
where cij is the probability of an emission from voxel j (activity λjt at time t) detected on LOR i, and Δt is the duration of time bin t (equal to the time resolution of the list mode data). For simplicity, attenuation, normalization, randoms and scatter are not included here, although this is a straightforward addition to the algorithm. The 1T model is:
where Pτ is the input function and the parameters are K1 and k2. The log likelihood function is:
The complete data space recognizes that the activity at any time is the sum of activities delivered to the system at earlier times, each with a different residence time. Thus, define the spatiotemporal complete data space Xijtτ to be the counts collected along LOR i in time bin t emitted from voxel j and where the input was delivered at time τ
From the complete data, the EM algorithm for the estimation of the parameters can be derived.
E-step, at iteration n:
M-step: maximization of conditional data likelihood,
Define the sensitivity image,
Taking partial derivatives with respect to K1j and k2j, and setting to zero, for K1j yields:
which leads to the K1j update equation
The left hand side is a function of k2 and the input function (assumed known) and the right hand side is a function of the data. Define the left hand side expression as a function of k2:
and obtain the k2j update equation
where H(k2) can be pre-computed when the input function is known. To assess the proposed algorithm, 2D tests and 4D tests were performed.
To update k2, the inverse of H(k2) is required. H(k2) is independent of measured data, so it can be computed when a plasma input is given, as is shown in Fig.1A. H(k2) was computed for a 30-min simulation (Fig. 1B). A fifth order polynomial log function can fit H(k2)’s inverse very well (Fig.1 C). The fitting error is shown in Fig.1D; errors of less then 1×10-4 in k2 were found. This is determined by precomputing H(k2) and then fitting the inverse functions with a polynomial log function.
Tests were performed for the case of 100 voxels, where gray matter (GM), white matter (WM) and basal ganglia (BG) were simulated. Fig.2 and Table I show the K1 and k2 values of the 1D phantom, which are typical of receptor kinetics. Fifty 30-min list mode data replicates were created by forward projection of Eq. 2 with a measured input function and with sampling, resolution and sensitivity (Q) comparable to that of the HRRT. Poisson noise was added to each replicate, corresponding to 6.3×105 events. Random and scatter events were not included. Values were reconstructed with the new algorithm and with the 2-step method (30×1 min frames). To fairly compare the two methods, for the 2-step method, nonlinear weighted least square method was utilized to estimate parametric values after the first step reconstruction, with weights are based on Noise Equivalent Counts (NEC).
A 10-cm diameter spherical 3D phantom (Fig.3) was created with 3 embedded regions, corresponding to GM (K1=0.55, VT=6), WM(K1=0.15, VT=3) and BG(K1=0.55, VT=12), with the same kinetics as the 2D test. One list mode file was simulated based on the HRRT and the new 4D method was integrated into the MOLAR software . Results were compared to the frame-based 2-step method (6×30 sec, 3×60 sec, 2×120 sec, 4×300 sec) using 30 subsets for both methods. Images were saved after each subset reconstruction update (sub-iteration).
Table I shows the mean percent bias and coefficients of variation (COV) for GM, WM and BG across the 50 realizations for frame based method and the proposed method after 60 iterations. 4D percent bias was comparable to frame-based values, and slightly smaller for VT, which is more the important parameter for brain receptor studies. 4D’s advantage was clearly demonstrated by COV reduction of ~ 30%.
Fig.4 shows 4D and frame-based results after 30 and 60 sub-iterations. Visually, for the K1 30 sub-iterations image, there is little difference between the two methods. But for VT, 4D is clearly better than frame-based method, especially for the high VT value region BG and low VT value region WM. However, for the K1 60 sub-iterations image, frame-based method is less noisy, but for the VT, 4D is still better. After 1 iteration (30 sub-iterations), the 4D method produced parametric images close to true images, which indicates that it may have faster convergence than frame-based method. Mean bias and COV across each region of interest (ROI’s) voxels were calculated from the one realization. Fig.5 shows the plots of bias vs. variability for the three ROIs, plotted at intervals of 10 sub-iterations (30~90). As sub-iteration number increases, the curves (and the noise) go up.
For K1, 4D has higher noise than frame-based method at the same iteration for all the three regions, and 4D is more biased for the GM region, although these biases were small. K1 bias was larger for the frame method in BG and WM. Like the 2D tests, 4D produce more accurate VT images than frame-based. In addition, unlike the frame-based method, 4D’s bias doesn’t change much with the increasing sub-iterations, which further suggests that it may converge faster.
In terms of computational complexity, the 4D method is somewhat more intensive than the frame method, because the computations in forward- and backprojection steps are more intensive. H(k2) and its inverse add an insignificant amount of computation. In the 4D method, handling of the whole list mode file (usually several billion events) is required, which requires more memory or more disk swapping.
In this paper, we introduced a method for the direct reconstruction of kinetic parameters from PET list mode data, which integrated the one tissue compartment model and the PET physical model and facilitated kinetic parameters estimation with the help of new spatiotemporal complete data space. Substantial more evaluation is required including simulations with motion correction, random correction and scatter correction and human data evaluations. Careful comparison of convergence between 4D and frame-based methods will be essential for fair comparisons of resolution and noise. In addition, it is necessary to optimize the software to decrease the computation time. We believe that this 4D EM method is a promising approach for kinetic parameter estimation directly from measured data.
We acknowledge the support of R01NS058360 (NINDS) and Siemens Medical Systems, Zhongdong Sun for programming and the staffs of Yale PET Center for the studies which formed the basis of this work.