Journal of Healthcare Engineering

Volume 2018, Article ID 5942873, 14 pages

https://doi.org/10.1155/2018/5942873

## Direct Parametric Maps Estimation from Dynamic PET Data: An Iterated Conditional Modes Approach

^{1}Dipartimento di Ingegneria dell’Informazione, University of Pisa, Pisa, Italy^{2}Fondazione Toscana G. Monasterio, Via G. Moruzzi 1, 56124 Pisa, Italy^{3}CNR Institute of Clinical Physiology, Via G. Moruzzi 1, 56124 Pisa, Italy

Correspondence should be addressed to Maria Filomena Santarelli; ti.rnc.cfi@leratnas

Received 20 December 2017; Revised 27 March 2018; Accepted 8 May 2018; Published 8 July 2018

Academic Editor: Emiliano Schena

Copyright © 2018 Michele Scipioni et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We propose and test a novel approach for direct parametric image reconstruction of dynamic PET data. We present a theoretical description of the problem of PET direct parametric maps estimation as an inference problem, from a probabilistic point of view, and we derive a simple iterative algorithm, based on the Iterated Conditional Mode (ICM) framework, which exploits the simplicity of a two-step optimization and the efficiency of an analytic method for estimating kinetic parameters from a nonlinear compartmental model. The resulting method is general enough to be flexible to an arbitrary choice of the kinetic model, and unlike many other solutions, it is capable to deal with nonlinear compartmental models without the need for linearization. We tested its performance on a two-tissue compartment model, including an analytical solution to the kinetic parameters evaluation, based on an auxiliary parameter set, with the aim of reducing computation errors and approximations. The new method is tested on simulated and clinical data. Simulation analysis led to the conclusion that the proposed algorithm gives a good estimation of the kinetic parameters in any noise condition. Furthermore, the application of the proposed method to clinical data gave promising results for further studies.

#### 1. Introduction

Dynamic positron emission tomography (PET) consists of the acquisition of a sequence of 3D images over time to follow the uptake and washout of an injected radiotracer in the imaged object. Hence, dynamic PET provides additional temporal information about tracer kinetics compared to static PET. The pharmacokinetic analysis allows, then, to estimate biologically relevant kinetic parameters from the measured concentration in tissue over time, that is, the tissue’s time activity curve (TAC). These parameters can provide a greater insight into the diagnosis of several diseases, helping clinicians in differentiating among different kind of lesions, and can assist during the follow-up of treatment response in a more specific way than the simple evaluation of the standard uptake value at a single time point [1].

There are two main approaches to characterize tracer kinetics from dynamic PET data: region-of-interest (ROI) kinetic modeling and parametric imaging [2, 3]. The ROI-based approach fits a kinetic model to the average TAC of a selected ROI. On the contrary, parametric imaging aims at estimating kinetic parameters for every voxel, thus providing a representation of their spatial distribution as parametric maps.

Parametric images have been found useful to enhance characterization of the regional heterogeneity. However, this technique is more computationally demanding, and more sensitive to noise in dynamic PET data than the ROI-based kinetic modeling [4]. The main difficulties related to estimating these maps are that a single voxel’s TAC is noisier than the average TAC used for ROI-based modeling and that the voxel-wise estimate needs to be performed in a sequential way, without the possibility to exploit information about spatial proximity as regularization constraint during the fitting of complex kinetic models. This can lead to different optimal parameter sets even for voxels belonging to the same tissue and supposedly sharing the same kinetic behavior; this produces potentially noisy and biased parametric maps. There is some research about the effects of spatial regularization on kinetic parameter estimation [5], but this approach is not commonly used because it requires huge parallelization power to process all the voxel in the 3D volume at the same time and because it could easily lead to biased maps, if the spatial smoothing is too strong.

The standard approach to process dynamic PET scans (4D PET) starts from the independent reconstructions of series PET images, acquired in a certain number of consecutive time frames. This can potentially lead to low count statistics, according to the chosen time scheme, which determines the length of each time frame [6]. Due to the ill-conditioning of the tomographic problem, iterative reconstruction algorithms, such as the maximum likelihood expectation maximization (MLEM) method, are a common choice, even if the filtered backprojection (FBP) method is sometimes preferred. The pharmacokinetic modeling is typically applied later, and thus, it is sensitive to noise in the previously reconstructed PET images. An efficient and unbiased estimate of the kinetic parameters of the chosen model requires properly taking into account the noise distribution of the reconstructed activity images [7] or a parametric images denoising by means of a spatial regularization [5]. However, for the frequently used iterative image reconstruction algorithms, this noise distribution is rarely known precisely and very hard to model with a known mathematical form because it is related to the number of iterations employed and other features of the reconstruction algorithm [8–10].

Postprocessing techniques such as filtering or ROI-based analysis (that, e.g., clusters voxel with similar kinetics) exist. However, all postreconstruction modeling approaches are fundamentally limited as they operate only indirectly on the reconstructed images so that the noise statistics cannot be accurately modeled. Inaccurate modeling of noise statistics may result in unreliable estimation of voxel-wise kinetic parameters, particularly if there are high levels of noise in the acquired data.

Given these considerations, the growing interest in fully 4D image reconstruction techniques can be easily understood. 4D reconstruction methods try to address the problems of noise characterization and limited counts in dynamic emission tomography by incorporating a theoretical model of the temporal behavior of the radiotracer directly into the image reconstruction algorithm [7, 11, 12]. This approach allows for adding physiologically meaningful constraints to the reconstruction process itself, and it provides reliable estimations of voxel-wise kinetic parameters directly from raw data by exploiting the same mathematical models normally used on a postreconstruction basis [2, 3].

It is well known that PET projection data can be described as independent Poisson realizations. Directly estimating kinetic parameters from sinogram data facilitates accurate modeling of noise distribution which follows Poisson statistics and gives an accurate compensation of the noise propagation, from projection measurement to kinetic fitting. This peculiar property has been shown [13–15] to provide a better bias-variance characteristic with respect to the conventional indirect approaches, both using linear and nonlinear kinetic models.

In this work, we propose a theoretical description of the problem of PET direct parametric maps estimation as an inference problem, from a probabilistic point of view. This assumption allows us to derive a simple two-step iterative optimization algorithm, based on the Iterated Conditional Mode (ICM) framework [16]. The resulting method is general enough to be flexible to an arbitrary choice of the kinetic model, and unlike many other solutions [17–19], it is natively capable to deal with nonlinear compartmental models without the need for linearization. For the testing purpose, we chose to couple the proposed optimization method with a two-tissue compartment model [2, 3], including an analytical solution to the kinetic parameters evaluation problem, based on an auxiliary parameter set, with the aim of reducing computational errors and approximations. The new method is tested on simulated and clinical data.

#### 2. Materials and Methods

##### 2.1. PET Data Model

Let us approximate the radiotracer activity within the region-of-interest of the subject’s body using a set of point sources , displaced on a regular voxel grid. Since emission events in the same voxel are not time correlated, their rate can be described as a Poisson distribution, with the expected value . The geometry of the acquisition system and attenuation determines the probability of a photon emitted by voxel being detected by line of response (LOR) . From the sum property of Poisson distribution, counts recorded in are, again, Poisson distributed, with the expected value , and measurements in each detector bin are independent, conditionally to activity [20, 21]. Keeping in mind that in dynamic PET imaging both activity image and sinogram counts are functions of time, we can express the probability to observe given as

The aim of direct parametric PET reconstruction is to generate parametric images , with number of model parameters, directly from the measured raw dynamic data and to use them to guide the image reconstruction process [7, 12, 13, 15, 22, 23]. If we define a vector , which contains the parameters related to voxel , and we identify with a generic kinetic model, which provides a theoretical representation of the time activity curve (TAC) for voxel over time, the relationship between the model and the reconstructed voxel TAC can be expressed as follows:

Irrespectively of the chosen kinetic model , (2) encodes the assumption that voxel intensities can be modeled as noisy realizations of a hidden dynamic process, and each voxel’s TAC can be parameterized using a kinetic model, defined on a set of parameters. Given (2), we decided to model the uncertainty using a probability distribution: for each time point , the corresponding value of has a Gaussian distribution with a mean equal to value of the model curve and standard deviation . Thus, we have

Given (1) and (2), we can explicit the connection between measured sinogram counts and model parameters through the following reinterpretation of the standard affine model of PET data:where is the estimate of random and scattered events occurring in the same LOR *i*, during time frame *m*. This allows us to directly express the Poisson likelihood in (1) as a , function of model parameters , instead of activity values *x*.

##### 2.2. ICM-Based 4D Reconstruction

Instead of directly deriving an optimization algorithm to maximize the likelihood , as done in other works about direct PET maps estimation [7, 12, 13, 15, 22, 23], we made use of the Iterated Conditional Modes (ICM) framework, introduced by Besag [16]. ICM consists of splitting the optimization problem by maximizing, in turn, each conditional probability, given the provisional estimates of the other variable: this procedure defines a single cycle of an iterative algorithm to estimate all the variables.

Considering the Poisson likelihood in (1) and the Gaussian likelihood in (3), we alternatively estimate (i) the parameters of the kinetic model with the highest probability given the activity and (ii) the activity with highest probability given the parameters and the PET projection data.(i)Given the provisional estimate of activity , we can adopt a maximum likelihood approach to maximize the logarithm of the likelihood function, which has arisen from the assumption of Gaussian noise distribution in (3). It is easy to see how maximizing the Gaussian log-likelihood is equivalent to minimizing the *sum of squares error function* between model and voxel’s TAC:which can be done using any conventional nonlinear least-squares (NLS) method (e.g., the Levenberg–Marquardt algorithm [23]).(ii)Given the updated parametric maps (i.e., parameter vector for each voxel of the volume of interest), we now want to maximize of (1) incorporating the result from the previous step. Inspired by the approach proposed by Wang and Qi [24, 25] (in their optimization transfer expectation maximization (OT-EM) algorithm for PET direct reconstruction) for transferring optimization of the Poisson likelihood from the image to parameters domain, we chose to maximize this conditional *pdf* as follows:where is the mean value of activity in voxel at time frame , given the updated estimate of model parameters performed in the previous step.

For the rest of this work, we will denote this concurrent optimization algorithm as Iterated Conditional Modes based Expectation Maximization (ICM-EM) reconstruction.

##### 2.3. Kinetic Modeling with a Two-Tissue Compartmental Model

The proposed inference framework was described in the previous section as independent of the choice of the kinetic model. However, we now need a proper expression for the theoretical model in (2) and (5), able to describe the behavior of the tracer in tissue. While linear models are often the preferred choice because of their computational efficiency and robustness, nonlinear ones could provide a greater insight into the biochemical properties of the various tissues [2, 3]. Most of these nonlinear models are based on compartments to describe the behavior of the tracer: each compartment identifies either a distinct physical location or a different chemical state of the tracer. The unknown parameters of the model are constant transfer rates, which describe the flow of the tracer among different compartments.

Following the compartmental model theory [3], the total tracer concentration in a tissue region can be modeled as follows:where is the fractional volume of blood in tissue, is the radioactive decay constant of the chosen tracer (e.g., for [18F]FDG, ), is the measured tracer concentration in plasma, is the concentration in the whole blood, and represents the expected impulse response function (IRF) in the compartment. The relationship between compartment concentrations is usually described by a set of ordinary differential equations (ODEs):where , and are the model constant parameters matrices, and denotes the system’s input function. In the common case of a two-tissue compartment model, we havewhere with and we distinguish between the concentration of the free and bound compartments of a two-tissue model. The system of ODE in (8) can be solved analytically in case of an impulsive input , and the impulse response function (IRF) iswhere . As shown by Gunn et al. [3], this solution can be further simplified introducing a set of auxiliary parameters :so that we can express the tissue IRF as a sum of two exponential functions:

##### 2.4. Modeling of the Input Function

Once we have the tissue IRF, from system theory we know that the output of our dynamic system with respect to a generic input can be computed by convolving (12) with the system input function, . Combining (7) with (12), we obtain

This means that tracer kinetic modeling with PET requires to measure the tracer time activity curves in both plasma and tissue to estimate physiological parameters.

In the present work, we decided to adopt a theoretical representation also for the input function, using a combination of exponential terms, that for [18F]FDG tracer is given by the following equation [26]:

##### 2.5. Analytical Tissue Compartment Modeling

The convolution operation in (13) is usually performed numerically, for each voxel in our volume, and this operation can be computationally and time expensive, given the size of the 4D volumes we are dealing with. It is interesting to note how this operation can be computed more efficiently when discretization is avoided [27].

Using Feng’s input function model (14), with an amplitude vector , and inverse time constant vector , we can derive an analytical solution to the convolution operation, so to avoid entirely the numerical convolution and temporal sampling of the input function. Expanding the convolution operator in the output function (13), substituting the input function model (14), and performing a series of algebraic manipulations, the output function can be expressed analytically as follows:where .

Adding to equation (15), the correction factors for radioactive decay and blood fraction in tissue, as in (7), we obtain the final version of our two-tissue analytic compartment model:

##### 2.6. Algorithm Implementation

The ICM-EM algorithm here proposed was implemented in Matlab (The MathWorks, Inc., Natick, Massachusetts, United States).

For the first step of the algorithm (5), we used the *lsqcurvefit* function, provided by Matlab with the Optimization Toolbox, for nonlinear least-squares fitting with a Levenberg–Marquardt algorithm [23]: we chose to constrain the search space for the kinetic constants between 0 and 10, and for the parameter, between 0 and 1. Default stopping conditions were used.

We, then, used the *Image Reconstruction Toolbox* (IRT), provided by Fessler [28], to implement the second step (6) of updating the activity estimate given the new estimate of parametric maps and sinogram data. We exploited this tool for system matrix generation, based on a square pixel basis and strip-integral detector model, and we applied the required modification to the IRT’s MLEM algorithm in order to correctly implement (6) for 2D single-slice dynamic PET data. Before reconstructing, images were initialized by setting all voxels inside the field of view (FOV) to 1.0 and all parameters to 0.1.

##### 2.7. Simulation

A numerical PET phantom was created starting from the Zubal brain phantom [29] image of 111 × 111 points, which includes white matter and gray matter. The field of view was 70 cm. We simulated an ideal measurement system with ideal detectors with equal sensitivity, efficiency, and detection cross section.

Dynamic data were generated by applying on each pixel the TAC curve that follows the kinetic parameter behavior according to the two-compartment model. The kinetic parameter values used for the simulation of both gray matter and white matter are shown in Table 1. The blood input function was modeled according to (14), with the following parameter values: , , and , and , , and .