Unmixing dynamic PET images with variable specific binding kinetics

Unmixing dynamic PET images with variable specific binding kinetics

Yanna Cruz Cavalcanti, , Thomas Oberlin, ,
Nicolas Dobigeon, , Simon Stute, Maria Ribeiro, Clovis Tauber, 
Y. C. Cavalcanti, T. Oberlin and N. Dobigeon are with University of Toulouse, IRIT/INP-ENSEEIHT, CNRS, 2 rue Charles Camichel, BP 7122, 31071 Toulouse Cedex 7, France (e-mail: {Yanna.Cavalcanti, Thomas.Oberlin, Nicolas.Dobigeon}@enseeiht.fr).S. Stute is with UMRS Inserm U1023 IMIV-CEA SHFJ, 91400 Orsay, France (e-mail: simon.stute@cea.fr).M. Ribeiro and C. Tauber are with UMRS Inserm U930 - Université de Tours, 37032 Tours, France (e-mail: {maria.ribeiro, clovis.tauber}@univ-tours.fr).Part of this work has been supported by CAPES.

Objective: To analyze dynamic positron emission tomography (PET) images, various generic multivariate data analysis techniques have been considered in the literature, such as principal component analysis (PCA), independent component analysis (ICA), factor analysis and nonnegative matrix factorization (NMF). Nevertheless, these conventional approaches neglect any possible nonlinear variations in the time activity curves describing the kinetic behavior of tissues with specific binding, which limits their ability to recover a reliable, understandable and interpretable description of the data. This paper proposes an alternative analysis paradigm that accounts for spatial fluctuations in the exchange rate of the tracer between a free compartment and a specifically bound ligand compartment. Methods: The method relies on the concept of linear unmixing, usually applied on the hyperspectral domain, which combines NMF with a sum-to-one constraint that ensures an exhaustive description of the mixtures. The spatial variability of the signature corresponding to the specific binding tissue is explicitly modeled through a perturbed component. Results: The performance of the method is assessed on both synthetic and real data and is shown to compete favorably when compared to other conventional analysis methods. Conclusion: The proposed method improved both factor estimation and proportions extraction for specific binding. Significance: Modeling the variability of the specific binding factor has a strong potential impact for dynamic PET image analysis.

Dynamic PET image, unmixing, brain imaging, factor analysis, matrix factorization, NMF.


I Introduction

Dynamic positron emission tomography (PET) is a non-invasive nuclear imaging technique that allows biological processes to be quantified and organ metabolic functions to be evaluated through the three-dimensional measure of the radiotracer concentration over time.

The analysis of dynamic PET images, in particular the quantification of the kinetic properties of the tracer, requires the extraction of tissue time-activity-curves (TACs) in order to estimate the parameters from compartmental modeling [1]. Nevertheless, PET images are corrupted by a prominent statistical noise and elementary TACs are mixed up due to the partial volume effect. Therefore, inferring rigorous and reliable information from these images still remains a challenging issue.

Several generic methods have been applied to estimate elementary TACs and their corresponding proportions from dynamic PET images. These techniques have different denominations depending on the application context, but all aim at tackling blind source separation (BSS) problems. For instance, Barber [2] and Cavailloles et al. [3] proposed matrix factorization-based PET analysis techniques, referred to as factor analysis of dynamic structures (FADS). An improvement of FADS taking into account the nonnegative physiological characteristic of PET images was proposed subsequently by Wu et al. [4] and Sitek et. al. [5]. Nonnegative matrix factorization (NMF) techniques pursue the same objective, under nonnegativity constraints on the latent factors to be recovered, and have been intensively applied in dynamic PET studies [6, 7, 8]. The works of Sitek et al. [9] and El Fakhri et al. [10] improved nonnegative FADS with a penalization that promoted non-overlapping regions in each voxel.

The approach proposed in this paper follows the same line as NMF or nonnegative FADS. It aims at decomposing each PET voxel TAC into a weighted combination of pure physiological factors, representing the elementary TAC associated with the different tissues present within the voxel. This factor modeling is enriched with a sum-to-one constraint to the factor proportions, so that they can be interpreted as tissue percentages within each voxel. In particular, this additional constraint explicitly solves the scaling ambiguity inherent to any NMF models, which has proven to increase robustness as well as interpretability. This BSS technique, referred to as unmixing or spectral mixture analysis, originates from the geoscience and remote sensing literature [11] and has proven its interest in other applicative contexts, such as microscopy [12] and genetics [13].

The linearity assumption for voxel decomposition underlined by the above-mentioned PET analysis methods can be envisioned in the light of compartment modeling, a tool widely employed to describe the kinetic behavior occurring within the voxel. The selection of a given model should be based on the radiotracer under study [14, 1] but most of the models consider that the measured signal in a given voxel is the sum of the comprising compartments.

However, factor TACs to be recovered cannot always be assumed to have constant kinetic patterns, as implicitly considered in conventional methods. Figure 1 depicts an example of a 2-tissue compartmental model, where the radioligand is assumed to move between three compartments: represents the radioligand concentration in arterial plasma, represents the free plus non-specific compartment, and represents the specifically bound compartment. The exchange between compartments are subject to rate constants () [15]. Considering both the 2-tissue and reference compartment models, the assumption of constant kinetic patterns seems appropriate for the blood compartment as well as non-specific binding tissues, since they present some homogeneity besides some perfusion difference (e.g. white matter versus gray matter). Therefore their contribution to the voxel TAC should be fairly proportional to the fraction of this type of tissue in the voxel. However, things get different regarding the specific binding class, as the TAC associated with this tissue is nonlinearly dependent on both the perfusion and the concentration of the radiotracer target. The spatial variation on target concentration is in part governed by differences in the and kinetic parameters, which nonlinearly modify the shape of the TAC characterizing this particular class. Muzi et al. [16] discussed the accuracy of parameter estimates for tumor regions and underlined high errors for the parameters related to specific binding, namely 26% for and 49% for . These results were further confirmed by Schiepers et al. [17]. More specifically, they studied the kinetics of lesioned regions that were tumor and treatment change predominant, showing that variations on and may allow for differentiation. Bai et al. [18] further discussed nonuniformity in intratumoral uptake and its impact on predicting treatment response and tumor aggressiveness. Nonetheless, this fluctuation phenomenon has not been taken into account by the decomposition models from the literature.


Fig. 1: Configuration of the classic three-compartment kinetic model used in many imaging studies.

The main motivation of this paper is to propose a more accurate description of the tissues and kinetics composing the voxels in dynamic PET images, in particular for those affected by specific binding. To this end, this work proposes to explicitly model the nonlinear variability inherent to the TAC corresponding to specific binding, by allowing the corresponding factor to vary spatially. This variation is approximated by a linear expansion over the atoms of a dictionary, which have been learned beforehand by conducting a principal component analysis on a learning dataset. Note that part of this work has been previously presented at the European Signal Processing Conference (EUSIPCO) [19].

The sequel of this paper is organized as follows. The proposed mixing-based analysis model is described in Section II. Section III presents the corresponding unmixing algorithm able to recover the factors, their corresponding proportions in each voxel and the variability maps. Simulation results obtained with synthetic data are reported in Section IV. Experimental results on real data are provided in Section V. Section VI concludes the paper.

Ii Problem Statement

Ii-a Specific binding linear mixing model (SLMM)

Consider voxels of a 3D dynamic PET image acquired at successive time-frames. First, we omit the spatial blurring induced by the point spread function (PSF) of the instrument and any measurement noise. The TAC in the voxel () over the time-frames is denoted . Akin to various BSS techniques and following the linear mixing model (LMM) for instance advocated in the PET literature by Barber et al. [2], each TAC is assumed to be a linear combination of elementary factors


where denotes the pure TAC of the tissue type and is the factor proportion of the tissue in the voxel. The factors () correspond to the kinetics of the radiotracer in a particular type of tissue in which they are supposed spatially homogeneous. For instance, the experiments conducted in this work and described in Sections IV and V consider types of tissues that fall into this category: the blood, the non-specific gray matter and the white matter.

Additional constraints regarding these sets of variables are assumed. First, since the elementary TACs are expected to be nonnegative, the factors are constrained as


Moreover, nonnegativity and sum-to-one constraints are assumed for all the factor proportions ()


For a given voxel indexed by , this sum-to-one constraint (3) enforces the mixing coefficients () to be interpreted as concentrations [20]. Therefore, all factor proportion vectors () lie inside the unit -simplex. Similarly, the TAC belongs to the convex set whose vertices are the columns of , represented by a -simplex in [11]. Fig. 2 shows an example of a simplex defined by elementary factors.

Fig. 2: Illustration of the simplex for a mixing matrix of 3 factors. The blue circles represent the vertices of the simplex, corresponding to the factors and red circles are the TACs.

More importantly, when factors are affected by possibly nonlinear and spatially varying fluctuations within the image, the conventional NMF-like linear mixing model (1) no longer provides a sufficient description of data. Therefore, over recent years, factor variability has received increased interest in the hyperspectral imagery literature as it allows changes on lightening and environment to be taken into account [21, 22]. Recently, Thouvenin et al. have proposed a perturbed LMM (PLMM) to further address this problem [23]. In the dynamic PET image framework, factor variability is expected to mainly affect the TAC associated with specific binding, denoted , while the possible variabilities in the TACs () related to tissues devoid of a specifically bound compartment are supposed weaker and neglected in this study. Since this so-called specific binding factor (SBF) is assumed to vary spatially, it will be spatially indexed. Thus, adapting the PLMM approach to our problem, the SBF in a given voxel will be modeled as a spatially-variant additive perturbation affecting a nominal and common SBF


where the additive term describes its spatial variability over the image. However, recovering the spatial fluctuation in each image voxel is a high-dimensional problem. To reduce this dimensionality, the variations will be assumed to lie inside a subspace of small dimension . As a consequence, similarly to the strategy followed in [24], the additive terms () are supposed to be approximated by the linear expansion


where the variability basis elements can be chosen beforehand, e.g., by conducting a PCA on a learning set composed of simulated or measured SBFs. The PCA aims at extracting the main variability patterns, while allowing for dimension reduction. Thus, the set of coefficients quantify the amount of variability in the voxel. The nominal SBF is also roughly estimated from the lowest values of this dataset and further fixed so that all the variability coefficients are nonnegative, in order to reduce correlation between the variability elements and the other tissue factors.

Combining the linear mixing model (1), the perturbation model (4) and its linear expansion (5), the voxel TACs are described according to the following so-called specific binding linear mixing model (SLMM)


To be fully comprehensive and motivated by the findings in [25], this work also proposes to explicitly model the PET scan point spread function (PSF), combining a deconvolution step jointly with parameter estimation. We will denote by the linear operator that computes the 3D convolution by some known and spatially invariant PSF, which leads to


where is a matrix containing the factor TACs, is a matrix composed of the factor proportion vectors, “” is the Hadamard point-wise product, is the matrix , is the matrix containing the basis elements used to expand the spatial variability of the SBF, is the matrix containing the intrinsic proportions, and is the matrix accounting for noise and mismodeling. Note that if and , the model in (7) reduces to the conventional linear mixing model generally assumed by factor model techniques like NMF and ICA.

While noises associated with count rates are traditionally modeled by a Poisson distribution [26], postprocessing corrections and filtering operated by modern PET systems significantly alter the nature of the noise corrupting the final reconstructed images. Modeling the noise on this final data is a highly challenging task [27]. However, as demonstrated in [28], pre-corrected PET data can be sufficiently approximated by a Gaussian distribution. As a consequence, in this work, the noise vectors () are assumed to be normally distributed. Moreover, without loss of generality, all vector components ( and ) will be assumed to be independent and identically distributed. This assumption seems to evade any spatial and temporal correlations that may characterize the noise generally affecting the reconstructed PET images [29]. However, the proposed model can be easily generalized to handle colored noise by weighting the model discrepancy measure according to the noise covariance matrix, as in [28]. Alternatively, after diagonalizing the noise covariance matrix, the PET image to be analyzed can undergo a conventional whitening pre-processing step [30, 31, 32].

In addition to the nonnegativity constraints applied to the elementary factors (2) and factor proportions (3), the intrinsic variability proportion matrix is also assumed to be nonnegative, mainly to avoid spurious ambiguity, i.e.,


where denotes the -matrix made of ’s and stands for a component-wise inequality. We accordingly fix the nominal SBF with a robust estimation of the TAC chosen as a lower bounding signature of a set of previously generated or measured SBF TACs. This means that a negative bias on the SBF is artificially introduced to model the spatially-varying SBF TACs (). This is alternatively compensated by a variability that is distorted by the same quantity but positively. This constraint is chosen to avoid a high correlation between the other factor TACs and when is allowed to be negative. Capitalizing on this model, the unmixing-based analysis of dynamic PET images is formulated in the next paragraph.

Ii-B Problem formulation

The SLMM (7) and constraints (2), (3) and (8) can be combined to formulate a constrained optimization problem. In order to estimate the matrices , , , a proper cost function is defined. The data-fitting term is defined as the Frobenius norm of the difference between the dynamic PET image and the proposed data modeling . This corresponds to the negative log-likelihood under the assumption of Gaussian noise. Since the problem is ill-posed and non-convex, additional regularizers become essential. In this paper, we propose to define penalization functions , and to reflect the available a priori knowledge on , and , respectively. The optimization problem is then defined as




where the parameters , and control the trade-off between the data fitting term and the penalties , and , described hereafter.

Ii-B1 Factor proportion penalization

The factor proportions representing the amount of different tissues are assumed to be spatially smooth, since neighboring voxels may contain the same tissues. We thus penalize the energy of the spatial gradient


where is the operator computing the first-order spatial finite differences. More details are reported in [33].

Ii-B2 Factor penalization

The chosen factor penalization benefits from the availability of rough factor TACs estimates . Thus, we propose to enforce similarity (in term of mutual Euclidean distances) between these primary estimates and the factor TACs to be recovered


Ii-B3 Variability penalization

The SBF variability is expected to affect only a small number of voxels, those belonging to the region containing the SBF. As a consequence, we propose to enforce sparsity via the use of the -norm, also known as the LASSO regularizer [34]


where is the norm. This penalty forces to be outside the high-uptake region, thus reducing overfitting.

Iii A PALM-based algorithm

Given the nature of the optimization problem (9), which is genuinely nonconvex and nonsmooth, the adopted minimization strategy relies on the proximal alternating linearized minimization (PALM) scheme [35]. PALM is an iterative, gradient-based algorithm which generalizes the Gauss-Seidel method. It consists in iterative proximal gradient steps with respect to , and and ensures convergence to a local critical point , and . The principle of PALM is briefly recalled in [33]. It is specifically instantiated for the unmixing-based kinetic component analysis considered in this paper. The resulting SLMM unmixing algorithm, whose main steps are described in the following paragraphs, is summarized in Algo. 1.

Input: , ,
1 while stopping criterion not satisfied do
3      proxλLBk∥.∥1(P+(Bk-γLBk∇BJ(Mk+1,Ak+1,Bk)))
Result: , ,
Algorithm 1 SLMM unmixing: main algorithm

Iii-a Optimization with respect to

A direct application of [35] under the constraints defined by (2) leads to the following updating rule


where is the projector onto the nonnegative set and the required gradient writes


Iii-B Optimization with respect to

Similarly to paragraph III-A, the factor proportion update is defined as the following


where is the projection on the set defined by the factor proportion constraints (3), which can be computed with efficient algorithms, see, e.g., [36]. The gradient can be computed as

with .

Iii-C Optimization with respect to

Finally, the updating rule for the variability coefficients can be written as

Iv Evaluation on Synthetic Data

Iv-a Synthetic data generation

To illustrate the accuracy of our algorithm, experiments are first conducted on synthetic data for which the ground truth of the main parameters of interest (i.e., factor TACs and factor proportion maps) is known. In the clinical PET framework, ground truth concerning the tracer kinetics and uptake is never completely known. Meanwhile, simulations benefit from an entire knowledge of the patient properties and kinetics, and their degree of complexity and details can be selected according to the purpose of the study. Furthermore, several simulations can be performed in a reasonable time.

Thus, experimentations are conducted on set of synthetic images of size . In these images, each voxel is constructed as a combination of pure TACs representative of the brain, which is the organ of interest in the present work: specific gray matter, pure blood or veins, pure white matter and non-specific gray matter. First, a high resolution dynamic PET numerical phantom with labeled regions of interest (ROIs) [37], has been used to create the ground truth for factor and proportions. In this phantom, all the distributions of the tracer per region have been extracted from a real dynamic PET acquisition and are observed in times of acquisition ranging from to minutes in a total period of minutes. Using a phantom extracted from real-life acquisitions brings to the synthetic image the complexity of real data. For instance, the pure TACs constitutive of this phantom image are created from the averaging of visually segmented regions and therefore are still mixed up due to partial volume. To simulate realistic variability of the SBF, a set of synthetic TACs generated through a realistic compartment-based model is used. More details can be found in [33].

The overall generation process is as follows:

  • The dynamic PET phantom has been first linearly unmixed using the N-FINDR [38] and SUnSAL [39] algorithms to select the ground-truth non-specific factor TAC and factor proportions , respectively. These factor TACs and corresponding factor proportion maps are depicted in Fig. 4 (left) and Fig. 5 (left), respectively.

  • Following a 2-tissue compartment-based model [40], a large database of SBF TACs has been generated by randomly varying the parameter (representing the specific binding rate of the radiotracer in the tissue). A PCA was conducted on this dataset, and an analysis of the eigenvalues led to the choice of a unique variability basis element (i.e., ), depicted in Fig. 3 (left).

  • The nominal SBF TAC is then chosen as the TAC of minimum area under the curve (AUC) among all the TACs of this database. This TAC is depicted in Fig. 3 (right, red curve).

  • The st row of the factor proportion matrix , namely was designed to locate the region associated with specific binding. Then, the matrix mapping the SBF variability in each voxel has been artificially generated. The high-uptake region was divided into subregions with non-zero coefficients , as shown in Fig. 6 (left), while these coefficients are set to outside the region affected with SBF. In each of these subregions, the non-zero coefficients have been drawn according to Gaussian distributions with a particular mean value and small variances. The spatially-varying SBFs in each region are then generated according to the model in (5) and (4). Some resulting typical SBF TACs are showed in Fig. 3.

Fig. 3: Left: variability basis element identified by PCA. Right: generated SBFs (blue) and the nominal SBF signature (red).

After this primary generation process, a PSF defined as a space-invariant and isotropic Gaussian filter with FWHMmm is applied to the output image. In brain imaging using a clinical PET scanner, this is an acceptable approximation, since the degradation of the scanner resolution mainly affects the borders of the field-of-view [41, 42]. Finally the measurements have been corrupted by a zero-mean white Gaussian noise with a signal-to-noise ratio dB, in agreement with a preliminary study conducted on the realistic replicas of [37], which showed that the SNR ranges from approximately dB on the earlier frames to dB on the latter ones.

Iv-B Compared methods

The results of the proposed algorithm have been compared to those obtained with several classical linear unmixing methods and other BSS techniques. The methods are recalled below with their most relevant implementation details.

NMF (no variability)

The NMF algorithm herein applied is based on multiplicative update rules using the Euclidean distance as cost function [43]. The stopping criterion is set to . To obtain a fair comparison mitigating scale ambiguity inherent to matrix factorization-like problem, results provided by the NMF have been normalized by the maximum value for the abundance, i.e.,


where denotes the th row of the estimated factor proportion matrix .

VCA (no variability)

The factor TACs are first extracted using the vertex component analysis (VCA) which requires pure voxels to be present in the analyzed images [44]. The factor proportions are subsequently estimated by sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL) [39].

LMM (no variability)

To appreciate the interest of explicitly modeling the spatial variability of the SBF, a depreciated version of the proposed SLMM algorithm is considered. More precisely, it uses the LMM (1) without allowing the SBF to be spatially varying. The stopping criterion, defined as , is set to . The values of the regularization parameter are reported in Table I.

SLMM (proposed approach)

As detailed in Section II-A, matrix is constrained to be nonnegative to increase accuracy. Consequently, the nominal SBF TAC is initialized as the TAC with the minimum AUC learned from the generated database to ensure a positive . The regularization parameters have been tuned to the values reported in Table I. As for the other approaches, the stopping criterion is set to .

Since the addressed problem is non-convex, these algorithms require an appropriate initialization. In this work, the factor TACs have been initialized as the outputs of a K-means clustering conducted on the PET image. These K-means TACs estimates are also considered for performance comparison.

0.010 0.010
0.010 0.010
- 0.020
0.001 0.001
TABLE I: factor proportion, factor and variability penalization parameters for LMM and SLMM with SNRdB

The performance of the algorithms has been accessed through the use of a normalized mean square error (NMSE) computed for each variable


where is the estimated variable and the corresponding ground truth. The NMSE has been measured for the following parameters: the factor proportions corresponding to the high-uptake region, the remaining factor proportions , the SBFs affected by the variability , the non-specific factor TACs and finally the variability factor proportion matrix .

Iv-C Hyperparameter influence

Considering the significant number of hyperparameters to be tuned in both LMM and SLMM approaches (i.e., , , ), a full sensitivity analysis is a challenging task, which is further complexified by the non-convex nature of the problem. To alleviate this issue, each parameter has been individually adjusted while the others have been set to zero. Several simulations empirically showed that the result is not very sensitive to the choice of parameters. The parameters have been tuned such that the total percentage of their corresponding term in the overall objective function does not surpass % of the total value of the function. Given the high level of noise corrupting the PET images, the hyperparameter associated with the factor proportions has been set so as to reduce the noise impact while avoiding too much smoothing. The factor TAC penalization hyperparameter results from a trade-off between the quality of the initial factor TAC estimates and the flexibility required by PALM to reach more accurate estimates. Finally the variability penalization has been tuned to achieve a compromise between the risks of capturing noise into the variability term (i.e., overfitting) and of losing information. While there are more automatized ways to choose the hyperparameter values (e.g., using cross-validation, grid search, random search and Bayesian estimation), these hyperparameter choices have seemed to be sufficient to assess the performance of the proposed method. The hyperparameter values used in LMM and SLMM are finally reported in Table I.

Iv-D Results

Fig. 4: Factor proportion maps of the 15th time-frame obtained for SNR=15dB corresponding to the specific gray matter, white matter, gray matter and blood, from top to bottom. The first 3 lines show a transaxial view while the last one shows a sagittal view. All images are in the same scale in
Fig. 5: TACs obtained for dB. For the proposed SLMM algorithm, the represented SBF TAC corresponds to the empirical mean of the estimated spatially varying SBFs .
Fig. 6: Ground-truth (left) and estimated (right) SBF variability.

The algorithms have been applied to different realizations of the noise to get reliable performance measures. Table II presents the NMSE averaged over these realizations for all algorithms and variables of interest while Table III presents their corresponding variances. The factor proportion maps recovered by the compared algorithms are shown in Fig. 4. Each row corresponds to a specific factor: SBF, white matter factor, non-specific gray matter factor, blood factor (from top to bottom, respectively). The six columns contain the factor proportion maps of the ground truth, and those estimated by K-means, NMF, VCA, LMM and the proposed SLMM (from left to right, respectively). A visual comparison suggests that the factor proportion maps obtained with LMM and SLMM are more consistent with the expected localization of each factor in the brain than VCA. Meanwhile, they are less noisy than the maps obtained by NMF. The estimated LMM and SLMM proportions maps are closer to the ground truth than both of them, particularly in the region affected by specific binding, as quantitatively shown in Table II. It can also be observed that the factor proportion maps obtained with the proposed SLMM approach present a higher contrast compared to LMM and other approaches, especially in the high-uptake region.

The maps of SLMM are also sharper compared to LMM. Additionally, it is also possible to see that NMF results for white matter are sharper but also more noisy than both LMM and SLMM approaches. However, for the specific gray matter, both LMM and SLMM approaches show sharper estimated factor proportion maps. Note that the sharpness of the factor proportions is not necessarily a good criterion of comparison. Indeed, factor analysis-based methods do expect to recover smooth maps that take into account the spilling part of partial volume effect, which is not considered within deconvolution. The aim of unmixing is not hard-clustering or classification.

The corresponding estimated factor TACs are shown in Fig. 5 where, for comparison purposes, the SBF depicted for SLMM is the empirical average over the whole set of spatially varying SBFs, as it is also the case for the SBF ground truth TACs. The best estimate of the SBF TAC seems to be obtained by the proposed SLMM approach, for which the TAC has been precisely recovered, as opposed to K-means, VCA and NMF. K-means provide the best estimate of the white matter TAC, closely followed by SLMM while NMF highly overestimates it. The best estimate of the non specific gray matter TAC has been obtained by VCA and NMF, even though it is slightly overestimated. It can be observed that SLMM and LMM have underestimated this factor TAC, which has been compensated with higher values in the corresponding factor proportion map. This underestimation results from the poor initialization of these algorithms by the K-means outputs. A more powerful initialization for the gray matter may provide better results. The factor TAC associated with blood is correctly estimated by SLMM, LMM, VCA and NMF.

The quantitative results in Table II confirm the preliminary findings drawn from the visual inspection of Fig. 4 and 5. The proposed method outperforms all the others for the estimation of , and . In particular, SLMM provides a very precise estimation of the mean SBF TAC with an NMSE of . In Fig. 5, the mean of the estimated SBF TACs is very close to the ground truth for LMM and SLMM but the the individual errors computed for each voxel demonstrate better performance obtained by SLMM. It also shows better results than K-means and NMF for , even though it is less effective but still competitive when compared to LMM and VCA.

K-means 0.567 0.669 0.120 0.442 -
VCA 0.547 0.481 0.517 0.248 -
NMF 0.512 0.558 0.517 0.133 -
LMM 0.437 0.473 0.349 0.148 -
SLMM 0.359 0.495 0.009 0.128 0.259
TABLE II: Mean of Normalized Mean Square Errors of estimated variables for K-means, VCA, NMF, LMM and SLMM with SNR=15dB
K-means () 0.043 0.225 0.015 6.136 -
VCA () 0.667 1.817 0.093 1.301 -
NMF () 0.010 0.381 0.448 1.486 -
LMM () 3.791 0.043 0.602 1.469 -
SLMM () 1.270 0.312 0.003 0.098 2.261
TABLE III: Variance of Normalized Mean Square Errors of estimated variables for K-means, VCA, NMF, LMM and SLMM with SNR=15dB

Taking into account the SBF variability allows the estimation of to be improved up to . Fig. 6 compares the actual variability factor proportions and those estimated by the proposed SLMM. This figure shows that the estimated non-zeros coefficients are correctly localized in the subregions characterized by some SBF variability. These non-zero values seem to be affected by some estimation inaccuracies, mainly due to the deconvolution. However, the estimation error still stays close to .

To summarize, in dynamic PET imaging, the TACs associated with the different types of tissues are always highly correlated. Consequently, the conventional LMM approach converges to poor local optima and is only guided by the factor similarity penalization defined by 12. On the contrary, SLMM shows better results since it explicitly takes the spatial variability of the SBF into account. VCA associated with SUnSAL shows acceptable results, in particular given the high level of noise affecting the PET images. Finally, as expected, SLMM presents better results for the voxels containing specific gray matter tissue.

V Evaluation on Real Data

V-a PET data acquisition

Fig. 7: factor proportion maps of the real PET image with []DPA-714 of a subject with a stroke. The first 3 lines show a transaxial view while the last one shows a sagittal view. From top to bottom: the specific gray matter, white matter, non-specific gray matter and blood.

To assess the behavior of the proposed approach when analyzing real dynamic PET images, the different methods have been applied to a dynamic PET image with [18F]DPA-714 of a stroke subject. Cerebral stroke is a severe and frequently occurring condition. While different mechanisms are involved in the stroke pathogenesis, there is an increasing evidence that inflammation, mainly involving the microglial and the immune system cells, account for its pathogenic progression. The [18F]DPA-714 is a ligand of the 18-kDa translocator protein (TSPO) for in vivo imaging, which is a biomarker of neuroinflammation. The subject was examined using an Ingenuity TOF Camera from Philips Medical Systems, seven days after the stroke.

The PET acquisition was reconstructed into a -voxels dynamic PET image with time-frames. The PET scan image registration time ranged from seconds to minutes over a minutes period. The voxel size was of mm. As for the experiments conducted on simulated data, the voxel TACs have been assumed to be mixtures of types of elementary TAC: specific binding associated with inflammation, blood, the non-specific gray and white matters. The K-means method was applied to the images to mask the cerebrospinal fluid and to initialize NMF, LMM and SLMM algorithms. A ground truth of the high-uptake tissue was manually labeled by an expert based on a magnetic resonance imaging (MRI) acquisition. The stroke region was segmented on this registered MRI image to define a set of voxels used to learn the variability descriptors by PCA. The nominal SBF has been fixed as the empirical average of the corresponding TACs with AUC comprised between the 5th and 10th percentile. The choice to use the average of a percentile instead of the minimum AUC TAC is motivated by the fact that, in this case, the learning set is corrupted by noise and partial volume effects.

Fig. 8: TACs obtained by estimation from the real image.

V-B Results

Figure 7 depicts the factor proportion maps estimated by the compared methods. The corresponding estimated factor TACs are shown in Fig. 8. The LMM and SLMM algorithms estimate four distinct TACs associated with different tissues, as expected. In the first row of Fig. 7, corresponding to SBF, all methods seem to correctly recover the main localization of the stroke area. However, the proposed SLMM approach identifies a significantly larger area. This result seems to be in better agreement with the stroke area identified in the MRI acquisition of the same patient (see st column of Fig. 9). Moreover, the specific gray matter factor proportion maps estimated by SLMM and K-means shows high values in the thalamus, which is a region known to present specific binding of [18F]DPA-714. Another remarkable result is the factor proportion maps for the blood. The sagittal view represented in the last row is in the exact center of the brain. Both NMF and SLMM recovers factor proportion maps that are in very good agreement with the superior sagittal sinus vein that passes on the higher part of the brain. On the contrary, VCA estimates two factors which seem to be mixtures of the vein TACs and other region TACs.

Fig. 9 depicts three different views of the AVC area identified by the expert on MRI acquisition (st column), the estimated specific gray matter factor proportions (nd-th columns) and the estimated corresponding variability (th column). This figure shows that SLMM provides sharper and more accurate maps, characterized by a larger area with high uptake, as in the MRI ground-truth. This figure also compares the map recovered by K-means, which is used as the initialization of the proposed SLMM, with the final SLMM estimate. It is possible to note an interesting improvement of the final estimate when compared to its initialization. This demonstrates that the method converges to an estimation of the specifically bound gray matter that is more accurate with the proposed model.

Fig. 9: From left to right: MRI ground-truth of the AVC area, SBF coefficient maps estimated by K-means, NMF, VCA, LMM, SLMM and SBF variability estimated by SLMM.

Vi Conclusion and Future Works

This paper introduced a new model to conduct factor analysis of dynamic PET images. It relied on the unmixing concept accounting for specific binding TACs variation. The method was based on the hypothesis that the variations within the SBF can be described by a small number of basis elements and their corresponding proportions per voxel. The resulting optimization problem is extremely non-convex with highly correlated factors and variability basis elements, which leads to a high number of spurious local optima for the cost function. However, the experiments conducted on synthetic data showed that the proposed method succeeded in estimating this variability, which improved the estimation of the specific binding factor and the corresponding proportions. For the other quantities of interest, the proposed approach compared favorably with state-of-the-art unmixing techniques. The proposed approach has many potential applications in dynamic PET imaging. It could be used for the segmentation of a region-of-interest, classification of the voxels, creation of subject-specific kinetic reference regions or even simultaneous filtering and partial volume correction. Besides exploring such applications of the method, future works should focus on the introduction of a Poisson-fitting measure of divergence used in the cost function, e.g. Kullback-Leibler, to better model noise frequently encountered in low rate PET data.


  • [1] R. B. Innis et al., “Consensus nomenclature for in vivo imaging of reversibly binding radioligands,” J. Cereb. Blood Flow Metab, vol. 27, no. 9, pp. 1533–1539, May 2007.
  • [2] D. C. Barber, “The use of principal components in the quantitative analysis of gamma camera dynamic studies,” Physics Med. Biol., vol. 25, no. 2, pp. 283–292, 1980.
  • [3] F. Cavailloles, J. P. Bazin, and R. Di Paola, “Factor analysis in gated cardiac studies,” J. Nuclear Med., vol. 25, p. 1067–1079, 1984.
  • [4] H.-M. Wu et al., “Factor analysis for extraction of blood time-activity curves in dynamic FDG-PET studies,” J. Nuclear Med., vol. 36, no. 9, pp. 1714–1722, September 1995.
  • [5] A. Sitek, E. V. R. Di Bella, and G. T. Gullberg, “Factor analysis with a priori knowledge—application in dynamic cardiac spect,” Physics Med. Biol., vol. 45, p. 2619–2638., 2000.
  • [6] J. S. Lee et al., “Non-negative matrix factorization of dynamic images in nuclear medicine,” in Proc. IEEE Nuclear Sci. Symp (NSS), 2001.
  • [7] M. Padilla, P.and López et al., “Nmf-svm based cad tool applied to functional brain images for the diagnosis of alzheimer’s disease,” in IEEE Trans. Med. Imag., vol. 31, no. 2, Feb. 2012, pp. 207–216.
  • [8] D. Schulz, “Non-negative matrix factorization based input function extraction for mouse imaging in small animal PET - comparison with arterial blood sampling and factor analysis,” J. Molecular Imag. Dynamics, vol. 02, 2013.
  • [9] A. Sitek, G. T. Gullberg, and R. H. Huesman, “Correction for ambiguous solutions in factor analysis using a penalized least squares objective,” IEEE Trans. Med. Imag., vol. 21, no. 3, p. 216–225, 2002.
  • [10] G. El Fakhri et al., “Quantitative dynamic cardiac pet using generalized factor and compartment analyses,” J. Nuclear Med., vol. 46, no. 8, pp. 1264–1271, August 2005.
  • [11] J. M. Bioucas-Dias et al., “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” IEEE J. Sel. Topics Appl. Earth Observations Remote Sens., vol. 5, no. 2, pp. 354–379, Apr 2012.
  • [12] Y. Huang et al., “Temporal dynamics of host molecular responses differentiate symptomatic and asymptomatic influenza A infection,” PLoS Genetics, vol. 8, no. 7, Aug. 2011.
  • [13] N. Dobigeon and N. Brun, “Spectral mixture analysis of EELS spectrum-images,” Ultramicroscopy, vol. 120, pp. 25–34, Sept. 2012.
  • [14] R. N. Gunn et al., “Parametric imaging of ligand-receptor binding in PET using a simplified reference region model,” NeuroImage, vol. 6, pp. 279–287, 1997.
  • [15] I. Häggström, B. J. Beattie, and C. R. Schmidtlein, “Dynamic PET simulator via tomographic emission projection for kinetic modeling and parametric image studies,” Medical Physics, vol. 43, no. 6, Part 1, pp. 3104–3116, 2016.
  • [16] M. Muzi et al., “Kinetic modeling of 3′-deoxy-3′-fluorothymidine in somatic tumors: Mathematical studies,” J. Nuclear Med., vol. 46, no. 2, pp. 371–380, 2005.
  • [17] C. Schiepers et al., “18F-fluorothymidine kinetics of malignant brain tumors,” Eur. J. Nuclear Med. Molecular Imag., vol. 34, no. 7, pp. 1003–1011, Feb 2007.
  • [18] B. Bai, J. Bading, and P. S. Conti, “Tumor quantification in clinical positron emission tomography,” Theranostics, vol. 3, no. 10, pp. 787–801, 2013.
  • [19] Y. C. Cavalcanti et al., “Unmixing dynamic PET images with a PALM algorithm,” in Proc. European Signal Process. Conf. (EUSIPCO), Aug. 2017, pp. 425–429.
  • [20] N. Keshava, “A survey of spectral unmixing algorithms,” Lincoln Laboratory Journal, vol. 14, pp. 56–78, 2003.
  • [21] A. Zare and K. Ho, “Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing,” IEEE Signal Process. Mag., vol. 31, no. 1, pp. 95–104, Jan 2014.
  • [22] A. Halimi, N. Dobigeon, and J.-Y. Tourneret, “Unsupervised unmixing of hyperspectral images accounting for endmember variability,” IEEE Trans. Image Process., vol. 24, no. 12, pp. 4904–4917, Dec 2015.
  • [23] P.-A. Thouvenin, N. Dobigeon, and J.-Y. Tourneret, “Hyperspectral unmixing with spectral variability using a perturbed linear mixing model,” IEEE Trans. Signal Process., vol. 64, no. 2, pp. 525–538, 2016.
  • [24] S.-U. Park, N. Dobigeon, and A. O. Hero, “Variational semi-blind sparse deconvolution with orthogonal kernel bases and its application to MRFM,” Signal Process., vol. 94, pp. 386–400, Jan. 2014.
  • [25] S. Henrot et al., “Does deblurring improve geometrical hyperspectral unmixing?” IEEE Trans. Image Process., vol. 23, no. 3, pp. 1169–1180, 2014.
  • [26] L. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imag., pp. 113–122, May 1982.
  • [27] D. Wilson, B. Tsui, and H. Barrett, “Noise properties of the em algorithm. ii. monte carlo simulations,” Physics Med. Biol., vol. 39, pp. 847–871, 1994.
  • [28] J. A. Fessler, “Penalized weighted least-squares image reconstruction for positron emission tomography,” IEEE Trans. Med. Imag., vol. 13, no. 2, pp. 290–300, June 1994.
  • [29] O. Tichy and V. Smidl, “Bayesian blind separation and deconvolution of dynamic image sequences using sparsity priors,” IEEE Trans. Med. Imag., 2015.
  • [30] T. Thireou et al., “Blind source separation for the computational analysis of dynamic oncological PET studies,” Oncol Rep, vol. 15, p. 1007–1012, 2006.
  • [31] E. Bullmore et al., “Colored noise and computational inference in neurophysiological (fMRI) time series analysis: Resampling methods in time and wavelet domains,” Human Brain Mapping, vol. 12, pp. 61–78, 2001.
  • [32] F. Turkheimer et al., “A linear wavelet filter for parametric imaging with dynamic PET,” IEEE Trans. Med. Imag., vol. 22, no. 3, pp. 289–301, Mar 2003.
  • [33] Y. C. Cavalcanti et al., “Unmixing dynamic PET images with variable specific binding kinetics – Supporting materials,” University of Toulouse, IRIT/INP-ENSEEIHT, France, Tech. Rep., July 2017. [Online]. Available: http://dobigeon.perso.enseeiht.fr/papers/Cavalcanti_TechReport_2017.pdf
  • [34] R. Tibshrani, “Regression shrinkage and selection via the lasso,” J. Roy. Stat. Soc., vol. 58, no. 1, pp. 267–288, 1996.
  • [35] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized minimization for nonconvex and nonsmooth problems,” Math. Programming, vol. 146, no. 1-2, pp. 459–494, Jul 2013.
  • [36] L. Condat, “Fast projection onto the simplex and the -ball,” Math. Programming, vol. 158, no. 1-2, pp. 575–585, Sep 2015.
  • [37] S. Stute et al., “Analytical simulations of dynamic PET scans with realistic count rates properties,” in Proc. IEEE Nuclear Sci. Symp. Med. Imag. Conf. (NSS-MIC), Nov. 2015.
  • [38] M. E. Winter, “N-findr: an algorithm for fast autonomous spectral end-member determination in hyperspectral data,” in Proc. SPIE Imaging Spectrometry V, M. R. Descour and S. S. Shen, Eds., vol. 3753, no. 1.   SPIE, 1999, pp. 266–275.
  • [39] J. M. Bioucas-Dias and M. A. T. Figueiredo, “Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing,” in Proc. IEEE GRSS Workshop Hyperspectral Image SIgnal Process.: Evolution in Remote Sens. (WHISPERS), 2010.
  • [40] M. E. Phelps, J. C. Mazziotta, and H. R. Schelbert, Positron Emission Tomography and Autoradiography (Principles and Applications for the Brain and Heart).   Raven, New York, 1986.
  • [41] A. Rahmim, J. Qi, and V. Sossi, “Resolution modeling in pet imaging: theory, practice, benefits, and pitfalls,” Med. Phys., vol. 40, no. 6, Jun 2013.
  • [42] A. Mehranian et al., “Pet image reconstruction using multi-parametric anato-functional priors,” Physics Med. Biol., vol. 62, no. 15, p. 5975, 2017.
  • [43] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in Proc. Neural Info. Process. Syst. (NIPS), 2000.
  • [44] J. Nascimento and J. Dias, “Vertex component analysis: a fast algorithm to unmix hyperspectral data,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 4, pp. 898–910, Apr 2005.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description