Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2901873

Formats

Article sections

Authors

Related links

IEEE Nucl Sci Symp Conf Rec (1997). Author manuscript; available in PMC 2010 July 12.

Published in final edited form as:

IEEE Nucl Sci Symp Conf Rec (1997). 2008; 4774103: 3625–3628.

doi: 10.1109/NSSMIC.2008.4774103PMCID: PMC2901873

NIHMSID: NIHMS111590

PET center, Yale University, New Haven, CT 06520 USA

Jianhua Yan: ude.elay@nay auhnaij; Beata Planeta-Wilson: ude.elay@nosliw-atenalp.ataeb; Richard E. Carson: ude.elay@nosrac.e.drahcir

See other articles in PMC that cite the published article.

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 (*V*_{T}).

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 [1] 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 [2] 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×10^{9} 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 [5] and Wang et al [6] 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 [15] from static to 4D. The physical model for the projection measurement *Y _{it}* on line of response (LOR)

$${Y}_{\mathit{it}}\sim \mathit{Poisson}\phantom{\rule{0.2em}{0ex}}\left(\mathrm{\Delta}t\mathrm{\sum}_{j}{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{\lambda}_{\mathit{jt}}\right)$$

(1)

where *c _{ij}* is the probability of an emission from voxel

$${\lambda}_{\mathit{jt}}={K}_{1j}\underset{0}{\overset{t}{\int}}{C}_{p}(\tau ){e}^{-{k}_{2j}(t-\tau )}d\tau ={K}_{1j}\mathrm{\sum}_{\tau}{C}_{p}(\tau ){e}^{-{k}_{2j}(t-\tau )}\mathrm{\Delta}\tau ={K}_{1j}\mathrm{\sum}_{\tau}{P}_{\tau}{e}^{-{k}_{2j}(t-\tau )}$$

(2)

where P_{τ} is the input function and the parameters are *K*_{1} and *k*_{2}. The log likelihood function is:

$$\mathrm{log}\phantom{\rule{0.2em}{0ex}}f=\mathrm{\sum}_{t}\left\{-\mathrm{\Delta}t\mathrm{\sum}_{j}\mathrm{\sum}_{t}{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{\lambda}_{\mathit{jt}}+{Y}_{\mathit{it}}\phantom{\rule{0.2em}{0ex}}\mathrm{log}\left(\mathrm{\Delta}t\mathrm{\sum}_{j}\mathrm{\sum}_{t}{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{\lambda}_{\mathit{jt}}\right)\right\}$$

(3)

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 *X _{ijtτ}* to be the counts collected along LOR

$${X}_{\mathit{ijt\tau}}\sim \mathit{Poisson}\phantom{\rule{0.2em}{0ex}}\left(\mathrm{\Delta}t{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{K}_{1j}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2j}(t-\tau )}\right)$$

(4)

From the complete data, the EM algorithm for the estimation of the parameters can be derived.

E-step, at iteration *n*:

$${N}_{\mathit{ijt\tau}}^{(n)}=E\left\{{X}_{\mathit{ijt\tau}}|{K}_{1j}^{(n)},{k}_{2j}^{(n)},Y\right\}={Y}_{\mathit{it}}\phantom{\rule{0.2em}{0ex}}\frac{\mathrm{\Delta}t{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{K}_{1j}^{(n)}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2j}^{(n)}(t-\tau )}}{\mathrm{\Delta}t\mathrm{\sum}_{j}\mathrm{\sum}_{\tau}{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{K}_{1j}^{(n)}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2j}^{(n)}(t-\tau )}}$$

(5)

M-step: maximization of conditional data likelihood,

$$\mathrm{log}\phantom{\rule{0.2em}{0ex}}g=\mathrm{\sum}_{i}\mathrm{\sum}_{j}\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\left\{-E\left({X}_{\mathit{ijt\tau}}\right)+{N}_{\mathit{ijt\tau}}\phantom{\rule{0.2em}{0ex}}\mathrm{log}\left(E\left({X}_{\mathit{ijt\tau}}\right)\right)\right\}+\mathrm{constant}$$

(6)

Define the sensitivity image,

$${Q}_{j}=\mathrm{\sum}_{i}{c}_{\mathit{ij}}$$

(7)

Taking partial derivatives with respect to *K*_{1j} and *k*_{2j}, and setting to zero, for *K*_{1j} yields:

$$\frac{\partial \mathrm{log}\phantom{\rule{0.2em}{0ex}}g}{\partial {K}_{1j}}=\mathrm{\sum}_{i}\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\left\{-{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2j}(t-\tau )}+\frac{{N}_{\mathit{ijt\tau}}^{(n)}}{{K}_{1j}}\right\}=0$$

(8)

which leads to the *K*_{1j} update equation

$${K}_{1j}^{(n+1)}=\frac{\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\mathrm{\sum}_{i}{N}_{\mathit{ijt\tau}}^{(n)}}{{Q}_{j}\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\mathrm{\Delta}t{P}_{\tau}{e}^{-{k}_{2j}^{(n+1)}(t-\tau )}}$$

(9)

For *k*_{2j},

$$\frac{\partial \mathrm{log}\phantom{\rule{0.2em}{0ex}}g}{\partial {k}_{2j}}=\mathrm{\sum}_{i}\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\left\{{c}_{\mathit{ij}}\phantom{\rule{0.2em}{0ex}}{K}_{1j}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2j}(t-\tau )}(t-\tau )-(t-\tau )\phantom{\rule{0.2em}{0ex}}{N}_{\mathit{ijt\tau}}^{(n)}\right\}=0$$

(10)

The *K*_{1j} update (Eq. 9) is inserted into Eq.10, which leads to

$$\frac{\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\phantom{\rule{0.2em}{0ex}}(t-\tau ){P}_{\tau}{e}^{-{k}_{2j}(t-\tau )}}{\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2j}(t-\tau )}}=\frac{\mathrm{\sum}_{i}\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\phantom{\rule{0.2em}{0ex}}(t-\tau )\phantom{\rule{0.2em}{0ex}}{N}_{\mathit{ijt\tau}}^{(n)}}{\mathrm{\sum}_{i}\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\phantom{\rule{0.2em}{0ex}}{N}_{\mathit{ijt\tau}}^{(n)}}$$

(11)

The left hand side is a function of *k*_{2} 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 *k*_{2}:

$$H({k}_{2})=\frac{\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2}(t-\tau )}(t-\tau )}{\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\phantom{\rule{0.2em}{0ex}}{P}_{\tau}{e}^{-{k}_{2}(t-\tau )}}$$

(12)

and obtain the *k*_{2j} update equation

$${k}_{2j}^{(n+1)}={H}^{-1}\phantom{\rule{0.2em}{0ex}}\left(\frac{\mathrm{\sum}_{i}\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\phantom{\rule{0.2em}{0ex}}(t-\tau )\phantom{\rule{0.2em}{0ex}}{N}_{\mathit{ijt\tau}}^{(n)}}{\mathrm{\sum}_{t}\mathrm{\sum}_{\tau}\mathrm{\sum}_{i}\phantom{\rule{0.2em}{0ex}}{N}_{\mathit{ijt\tau}}^{(n)}}\right)$$

(13)

where *H*(*k*_{2}) can be pre-computed when the input function is known. To assess the proposed algorithm, 2D tests and 4D tests were performed.

To update *k*_{2}, the inverse of *H*(*k*_{2}) is required. *H*(*k*_{2}) is independent of measured data, so it can be computed when a plasma input is given, as is shown in Fig.1A. *H*(*k*_{2}) was computed for a 30-min simulation (Fig. 1B). A fifth order polynomial log function can fit *H*(*k*_{2})’*s* inverse very well (Fig.1 C). The fitting error is shown in Fig.1D; errors of less then 1×10^{-4} in *k*_{2} were found. This is determined by precomputing *H*(*k*_{2}) and then fitting the inverse functions with a polynomial log function.

A: Measured plasma input, B: *H*(*k*_{2}), C: Inverse log function: *k*_{2}~*Log*(*H*(*k*_{2})), D: *k*_{2} fitting error with polynomial (note change of scale)

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 *K*_{1} and *k*_{2} 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×10^{5} 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 (*K _{1}*=0.55,

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 *V*_{T}, 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 *K*_{1} 30 sub-iterations image, there is little difference between the two methods. But for *V*_{T}, 4D is clearly better than frame-based method, especially for the high *V*_{T} value region BG and low *V*_{T} value region WM. However, for the *K _{1}* 60 sub-iterations image, frame-based method is less noisy, but for the

The results of two methods after 30 (left side: A, B, C and D) and 60 sub-iterations (left side: E, F, G and H). A. *K*_{1} image (0~2) for 4D B. *K*_{1} image (0~2) for frame-based method, C. *V*_{T} mean image (0~30) for 4D, D. *V*_{T} mean image (0~30) for frame-based **...**

COV versus Bias curves for 4D (solid line) and frame-based method (circled line). As sub-iterations goes up (30~90), the curves go up

For *K _{1}*, 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.

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*(*k*_{2}) 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.

1. Snyder DL. Parameter estimation for dynamic studies in emission-tomography systems having list-mode data. IEEE Trans Nucl Sci. 1984;31:925–31.

2. Carson RE, Lange K. The EM parametric image reconstruction algorithm. J Am Statist Assoc. 1985;80:20–22.

3. Kamasak ME, Bouman CA, Morris ED, Sauer K. Direct reconstruction of kinetic parameter images from dynamic PET data. IEEE Trans Med Imaging. 2005;24:636–650. [PubMed]

4. Yetik S, Qi J. Direct estimation of kinetic parameters from the sinogram with an unknown blood function. IEEE Biomed Img Symp Macro to Nano. 2006;1:295–298.

5. Tsoumpas C, Turkheimer FE, Thielemans K. Study of direct and indirect parametric estimation methods of linear models in dynamic positron emission tomography. Med Phys. 2008;35:1299–309. [PubMed]

6. Wang G, Fu L, Qi J. Maximum a posteriori reconstruction of the Patlak parametric image from sinograms in dynamic PET. Phys Med Biol. 2008;53:593–604. [PubMed]

7. Matthews J, Bailey DL, Price P, Cunningham VJ. The direct calculation of parametric images from dynamic PET data using maximum-likelihood iterative reconstruction. Phys Med Biol. 1997;42:1155–73. [PubMed]

8. Meikle SR, Matthews JC, Cunningham VJ, Bailey DL, Livieratos L, Jones T, Price P. Parametric image reconstruction using spectral analysis of PET projection data. Phys Med Biol. 1998;43:651–66. [PubMed]

9. Nichols TE, Qi J, Asma E, Leahy RM. Spatiotemporal reconstruction of list-mode PET data. IEEE Trans Med Imaging. 2002;21:396–404. [PubMed]

10. Reader AJ, Sureau FC, Comtat C, Trebossen R, Buvat I. Joint estimation of dynamic PET images and temporal basis functions using fully 4D ML-EM. Phys Med Biol. 2006;51:5455–74. [PubMed]

11. Li Q, Asma E, Ahn S, Leahy RM. A fast fully 4D incremental gradient reconstruction algorithm for list mode PET data. IEEE Trans Med Imaging. 2007;26:58–67. [PubMed]

12. Verhaeghe J, D’Asseler Y, Vandenberghe S, Staelens S, Lemahieu I. An investigation of temporal regularization techniques for dynamic PET reconstructions using temporal splines. Med Phys. 2007;34:1766–78. [PubMed]

13. Reader AJ, Matthews JC, Sureau FC, Comtat C, Trebossen R, Buvat I. Fully 4D image reconstruction by estimation of an input function and spectral coefficients. IEEE Nucl Sci Symp Med Im Conf. 2007;4:3260–7.

14. Verhaeghe J, Ville DVD, Khalidov I, D’Asseler Y, Lemahieu I, Unser M. Dynamic PET reconstruction using wavelet regularization with adapted basis functions. IEEE Trans Med Imaging. 2008;27:943–59. [PubMed]

15. Lange K, Carson RE. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomogr. 1984;8:306–316. [PubMed]

16. Carson RE, Barker WC, Johnson JL. Design of a motion-compensation OSEM list-mode algorithm for resolution-recovery reconstruction for the HRRT. IEEE Nucl Sci Symp Med Im Conf. 2003;5:3281–3285.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |