Unmixing dynamic PET images with variable specific binding kinetics
Abstract
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 sumtoone 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.
IEEEexample:BSTcontrol
I Introduction
Dynamic positron emission tomography (PET) is a noninvasive nuclear imaging technique that allows biological processes to be quantified and organ metabolic functions to be evaluated through the threedimensional 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 timeactivitycurves (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 factorizationbased 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 nonoverlapping 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 sumtoone 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 abovementioned 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 2tissue compartmental model, where the radioligand is assumed to move between three compartments: represents the radioligand concentration in arterial plasma, represents the free plus nonspecific compartment, and represents the specifically bound compartment. The exchange between compartments are subject to rate constants () [15]. Considering both the 2tissue and reference compartment models, the assumption of constant kinetic patterns seems appropriate for the blood compartment as well as nonspecific 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.
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 mixingbased 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
Iia Specific binding linear mixing model (SLMM)
Consider voxels of a 3D dynamic PET image acquired at successive timeframes. 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 timeframes 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
(1) 
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 nonspecific 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
(2) 
Moreover, nonnegativity and sumtoone constraints are assumed for all the factor proportions ()
(3) 
For a given voxel indexed by , this sumtoone 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.
More importantly, when factors are affected by possibly nonlinear and spatially varying fluctuations within the image, the conventional NMFlike 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 socalled 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 spatiallyvariant additive perturbation affecting a nominal and common SBF
(4) 
where the additive term describes its spatial variability over the image. However, recovering the spatial fluctuation in each image voxel is a highdimensional 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
(5) 
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 socalled specific binding linear mixing model (SLMM)
(6) 
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
(7) 
where is a matrix containing the factor TACs, is a matrix composed of the factor proportion vectors, “” is the Hadamard pointwise 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], precorrected 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 preprocessing 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.,
(8) 
where denotes the matrix made of ’s and stands for a componentwise 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 spatiallyvarying 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 unmixingbased analysis of dynamic PET images is formulated in the next paragraph.
IiB 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 datafitting 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 loglikelihood under the assumption of Gaussian noise. Since the problem is illposed and nonconvex, 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
(9) 
with
(10) 
where the parameters , and control the tradeoff between the data fitting term and the penalties , and , described hereafter.
IiB1 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
(11) 
where is the operator computing the firstorder spatial finite differences. More details are reported in [33].
IiB2 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
(12) 
IiB3 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]
(13) 
where is the norm. This penalty forces to be outside the highuptake region, thus reducing overfitting.
Iii A PALMbased 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, gradientbased algorithm which generalizes the GaussSeidel 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 unmixingbased 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.
Iiia Optimization with respect to
IiiB Optimization with respect to
IiiC Optimization with respect to
Finally, the updating rule for the variability coefficients can be written as
Iv Evaluation on Synthetic Data
Iva 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 nonspecific 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 reallife 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 compartmentbased 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 NFINDR [38] and SUnSAL [39] algorithms to select the groundtruth nonspecific 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 2tissue compartmentbased 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 highuptake region was divided into subregions with nonzero 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 nonzero coefficients have been drawn according to Gaussian distributions with a particular mean value and small variances. The spatiallyvarying 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.
After this primary generation process, a PSF defined as a spaceinvariant 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 fieldofview [41, 42]. Finally the measurements have been corrupted by a zeromean white Gaussian noise with a signaltonoise 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.
IvB 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 factorizationlike problem, results provided by the NMF have been normalized by the maximum value for the abundance, i.e.,
(17) 
where denotes the th row of the estimated factor proportion matrix .
VCA (no variability)
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 IIA, 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 nonconvex, these algorithms require an appropriate initialization. In this work, the factor TACs have been initialized as the outputs of a Kmeans clustering conducted on the PET image. These Kmeans TACs estimates are also considered for performance comparison.
LMM  SLMM  

0.010  0.010  
0.010  0.010  
  0.020  
0.001  0.001 
The performance of the algorithms has been accessed through the use of a normalized mean square error (NMSE) computed for each variable
(18) 
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 highuptake region, the remaining factor proportions , the SBFs affected by the variability , the nonspecific factor TACs and finally the variability factor proportion matrix .
IvC 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 nonconvex 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 tradeoff 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 crossvalidation, 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.
IvD Results
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, nonspecific 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 Kmeans, 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 highuptake 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 analysisbased 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 hardclustering 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 Kmeans, VCA and NMF. Kmeans 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 Kmeans 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 Kmeans and NMF for , even though it is less effective but still competitive when compared to LMM and VCA.
Kmeans  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 
Kmeans ()  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 
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 nonzeros coefficients are correctly localized in the subregions characterized by some SBF variability. These nonzero 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
Va PET data acquisition
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]DPA714 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]DPA714 is a ligand of the 18kDa 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 timeframes. 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 nonspecific gray and white matters. The Kmeans method was applied to the images to mask the cerebrospinal fluid and to initialize NMF, LMM and SLMM algorithms. A ground truth of the highuptake 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.
VB 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 Kmeans shows high values in the thalamus, which is a region known to present specific binding of [18F]DPA714. 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 (ndth 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 groundtruth. This figure also compares the map recovered by Kmeans, 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.
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 nonconvex 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 stateoftheart unmixing techniques. The proposed approach has many potential applications in dynamic PET imaging. It could be used for the segmentation of a regionofinterest, classification of the voxels, creation of subjectspecific 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 Poissonfitting measure of divergence used in the cost function, e.g. KullbackLeibler, to better model noise frequently encountered in low rate PET data.
References
 [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 timeactivity curves in dynamic FDGPET 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., “Nonnegative 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., “Nmfsvm 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, “Nonnegative 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. BioucasDias et al., “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regressionbased 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 spectrumimages,” Ultramicroscopy, vol. 120, pp. 25–34, Sept. 2012.
 [14] R. N. Gunn et al., “Parametric imaging of ligandreceptor 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â²deoxy3â²fluorothymidine in somatic tumors: Mathematical studies,” J. Nuclear Med., vol. 46, no. 2, pp. 371–380, 2005.
 [17] C. Schiepers et al., “18Ffluorothymidine 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 semiblind 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 leastsquares 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/INPENSEEIHT, 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. 12, pp. 459–494, Jul 2013.
 [36] L. Condat, “Fast projection onto the simplex and the ball,” Math. Programming, vol. 158, no. 12, 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. (NSSMIC), Nov. 2015.
 [38] M. E. Winter, “Nfindr: an algorithm for fast autonomous spectral endmember 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. BioucasDias 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 multiparametric anatofunctional priors,” Physics Med. Biol., vol. 62, no. 15, p. 5975, 2017.
 [43] D. D. Lee and H. S. Seung, “Algorithms for nonnegative 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.