Voxelwise Principal Component Analysis of Dynamic [S-Methyl-11C]Methionine PET Data in Glioma Patients

Simple Summary Recent works on dynamic amino acid positron emission tomography (PET) imaging of gliomas have highlighted characteristic behaviors of time-activity curves (TACs) extracted from the whole tumor w.r.t. the grade, genotype, and outcome. However, gliomas are known to be highly heterogeneous tumors. Here, we aim at highlighting similar dynamic behaviors at the voxel level within the tumor volume in [S-methyl-11C]methionine PET data of 33 glioma patients using principal component analysis (PCA). The PCA model was derived from TACs of 20 patients and subsequently applied to 13 other patients in whom our approach was shown to outperform classical pharmacokinetic modeling to this end. Our parameter-free approach provides additional parametric maps from dynamic methionine PET scans with little modification of the routine protocol and no arterial sampling. This early methodological work paves the way for various clinical studies on glioma heterogeneity with applications for treatment planning and response evaluation. Abstract Recent works have demonstrated the added value of dynamic amino acid positron emission tomography (PET) for glioma grading and genotyping, biopsy targeting, and recurrence diagnosis. However, most of these studies are based on hand-crafted qualitative or semi-quantitative features extracted from the mean time activity curve within predefined volumes. Voxelwise dynamic PET data analysis could instead provide a better insight into intra-tumor heterogeneity of gliomas. In this work, we investigate the ability of principal component analysis (PCA) to extract relevant quantitative features from a large number of motion-corrected [S-methyl-11C]methionine ([11C]MET) PET frames. We first demonstrate the robustness of our methodology to noise by means of numerical simulations. We then build a PCA model from dynamic [11C]MET acquisitions of 20 glioma patients. In a distinct cohort of 13 glioma patients, we compare the parametric maps derived from our PCA model to these provided by the classical one-compartment pharmacokinetic model (1TCM). We show that our PCA model outperforms the 1TCM to distinguish characteristic dynamic uptake behaviors within the tumor while being less computationally expensive and not requiring arterial sampling. Such methodology could be valuable to assess the tumor aggressiveness locally with applications for treatment planning and response evaluation. This work further supports the added value of dynamic over static [11C]MET PET in gliomas.


Introduction
Gliomas are the most common primary brain tumors and are associated with poor prognosis. Glioma diagnosis and follow-up usually rely on magnetic resonance imaging (MRI), although the addition of positron emission tomography (PET) with radio-labeled amino acids such as [S-methyl- 11 C]methionine ([ 11 C]MET) has been shown to provide complementary information for tumor delineation [1] and characterization [2,3], as well as for biopsy [4][5][6] and therapy [7] planning. Whereas clinical amino acid PET imaging of gliomas is almost exclusively based on static acquisitions, the added value of dynamic PET acquisitions has been demonstrated for tumor grading and genotyping [8][9][10][11][12][13], biopsy targeting [9], and recurrence diagnosis [14]. Aside from a longer acquisition time, the main limitation of dynamic PET imaging lies in the difficulty of extracting robust and clinically relevant features from noisy time-activity curves (TACs).
Previous works on dynamic PET imaging of gliomas with O-(2-[ 18 F]fluoroethyl)-Ltyrosine ([ 18 F]FET)-another amino acid PET tracer equivalent to [ 11 C]MET [15]-have highlighted differences in uptake dynamics between high-grade gliomas (HGGs) and low-grade gliomas (LGGs) by visually labeling mean tumor TACs as 'increasing' or 'decreasing'. It has been shown that a fast increasing then progressively decreasing mean TAC is characteristic of HGGs whereas a slowly increasing mean TAC is rather observed in LGGs [8,10], and that foci with a decreasing TAC should be taken into account for surgery guidance [9]. These interesting findings however have some methodological limitations since TAC labeling does not allow continuous quantification of the dynamic behavior and voxelwise extension might become challenging for large amounts of noisy data.
Recently, time-to-peak (TTP) has been investigated as a dynamic feature of interest for quantitative characterization of the mean TAC in gliomas, with promising results for glioma grading [10][11][12][13]16] and recurrence diagnosis [14] in dynamic [ 18 F]FET PET. TTP has the advantage of being easily computed and reflects to some extent the 'increasing' or 'decreasing' behavior of the TACs. However, this parameter is highly sensitive to data noise and depends on the PET reconstruction framing, hence TTP values lie in a discrete range of arbitrarily chosen times.
Pharmacokinetic (PK) modeling is the gold standard method for dynamic PET data analysis. PK modeling relies on compartmental models whose kinetic parameters are estimated from the observed TACs given an arterial input function (AIF) (i.e., the TAC of arterial blood used as an input for the model). This method has the great advantage of providing biologically interpretable kinetic parameters and has been previously used for glioma delineation in 2-deoxy-2-[ 18 F]fluoro-D-glucose [17] as well as for glioma genotyping in [ 18 F]FET PET [11,18]. However, direct kinetic parameters fitting has two major limitations. First, it requires the user to provide an AIF, which can be either measured from arterial blood samples or extracted from large vessels appearing in the image. Nevertheless, arterial sampling is an invasive procedure and is inconvenient in clinical practice. Image-derived input functions (IDIFs), on the other hand, are affected by partial volume effects inherent to PET imaging and by blood tracer metabolites, with direct impact on the kinetic parameters estimation. Second, kinetic parameters fitting is a computationally expensive process and its voxelwise extension may result in substantially long computation times with increased impact of data noise. Alternative methods have been proposed for PK analysis of dynamic PET data such as graphical analysis or reference tissue models, the latter not requiring an AIF. However, these methods generally provide only macro or relative kinetic parameters. A comparison of commonly used PK analysis methods for parametric map extraction from dynamic [ 18 F]FET PET data in diffuse gliomas has been performed recently by Koopman and colleagues [19].
Principal component analysis (PCA) is a commonly used unsupervised multivariate analysis technique aiming at reducing high dimensional data space into a reduced number of components that best explain the observed data variance. The use of PCA for dynamic PET data analysis has been intensively studied by means of simulations, though surprisingly few clinical applications have seemed to emerge from these works. One major limitation of this technique is its limited ability to separate signal from noise for high noise levels or non-Gaussian noise distributions. The impact of noise level and distribution on PCA of dynamic PET images has been previously studied using synthetic imaging datasets by Pedersen and colleagues [20] and Šámal and colleagues [21].
Most of the published dynamic PET studies in gliomas, including those mentioned above, are only concerned with the mean TAC inside one or several pre-delineated volume(s) of interest (VOI(s)). However, gliomas are known to be highly heterogeneous tumors that may comprise multiple subregions with varying genotype, proliferation potential, aggressiveness, hypoxia level, and treatment-resistance abilities. Though computationally more expensive and more prone to noise, voxelwise analysis of dynamic PET data could instead provide a valuable intra-tumor insight. Voxelwise extensions of TTP analysis [16] and PK modeling [18] have been recently proposed but are still prone to the same limitations as their regionwise counterpart, with increased impact of noise and stability issues reported for PK modeling [18].
Besides, most dynamic PET studies rely on a limited number of non-uniformly sampled and non-overlapping frames with variable length and no prior motion correction. However, uniform TAC sampling with a large number of short frames could potentially benefit the dynamic analysis by providing a higher temporal resolution and removing the need to select an arbitrary irregular framing, at the expense of data noise. In this sense, reconstruction of overlapping frames leads to a good comprise between temporal resolution and count statistics and has been previously used for preclinical cardiac PET imaging [22]. Having access to high dimensional data-even noisy-with highly correlated variables is also particularly suitable for PCA.
In this work, we investigate the ability of PCA to extract meaningful dynamic features from a large number of uniformly sampled and motion-corrected [ 11 C]MET PET frames in glioma patients. The emphasis is placed on the ability to quantify at the voxel level of the 'increasing' versus 'decreasing' behavior of TACs previously reported by Pöpperl and colleagues for the whole metabolic tumor [8] since it is expected to be related to the local aggressiveness of the tumor. To this extent, we first assess the robustness of the proposed methodology to noise by means of realistic numerical simulations. We then compare the derived parametric maps to these obtained from classical voxelwise PK analysis. We conclude that PCA outperforms PK modeling in discriminating the characteristic dynamic behaviors reported previously while not requiring arterial sampling. Our results also support the added value of dynamic over static analysis of [ 11 C]MET PET data in gliomas, as previously demonstrated for [ 18 F]FET.

Image Acquisition and Reconstruction
In total, 33 glioma patients (20 males, 28 surgically treated, of median age 55 year) admitted to our institution for a diagnosis or follow-up [ 11 C]MET PET scan were enrolled in this study. Patient's clinical data, 2016 WHO classification, and undergone treatments at imaging time are provided for each lesion in Table A1. The patient cohort was further split into PCA model-building (n = 20, patients 1 to 20 in Table A1) and evaluation (n = 13, patients 21 to 33 in Table A1) sets (see Section 2.8). All patients underwent a 30 min and 30 s PET acquisition started 30 s before the intravenous injection of 287-555 MBq of [ 11 C]MET. All acquisitions were performed on a Vereos digital PET-CT scanner (Philips Healthcare, Best, The Netherlands) with an axial and trans-axial resolution of 4.1 mm at 1 cm from the field-of-view center. In total, 906 overlapping frames of 20 s, spaced by 2 s with a voxel size of 2 mm × 2 mm × 2 mm were reconstructed from the LIST files using time-of-flight ordered subset expectation maximization (TOF-OSEM, 10 subsets and 3 iterations) with computed tomography-based attenuation correction (CTAC). No post-reconstruction filter was applied to the dynamic frames. A routine static PET image (20-27 min post-injection (p.i.)) was also reconstructed for each patient (2 mm × 2 mm × 2 mm-CTAC TOF-OSEM, 15 subsets and 3 iterations-no filter).

Motion Correction
A systematic and progressive drift of the patients' head throughout the acquisition was observed in our dataset, imputed to the low stiffness of the scanner head support. Although motion between two consecutive frames is almost unnoticeable, it turned out that the patients' head had sunk of up to 1.5 mm over 30 min of acquisition. Such movement substantially impacts TACs at the voxel level and had to be corrected prior to the analysis. Since frame-by-frame registration is not suitable for short frames with a low signal-tonoise ratio (SNR), the following approach was used: Frames 36 to 906 (i.e., 40 s to 30 min p.i.) were grouped into consecutive blocks of 31 frames spaced by 30 frames. The mean frame of each block was computed, referred to as a 'long frame.' The last long frame was used as a reference for rigid registration of the other long frames by mutual information maximization [23]. Each computed rigid transform, expressed as a versor, was assigned to the mid-frame time point of the corresponding long frame. For each short frame, a rigid transform at mid-frame time was finally linearly interpolated from the known long frame transforms, which can be trivially performed when expressing transforms as versors. All registrations were performed in Python using SimpleITK [24].

Image-Derived Input Function Extraction
A blood input function was extracted from each registered volume for the pharmacokinetic analyses. A vascular image (0-120 s p.i.) was first computed by averaging registered frames 16 to 66. Two regions of interest (ROIs) covering the petrous (horizontal) segment of each internal carotid artery was then manually drawn on the early frame. The blood input function was finally obtained by averaging the TACs of the 10 brightest voxels on the vascular image within both ROIs. Brighter voxels are indeed considered to be the less impacted by partial volume effect (PVE), as suggested previously [25,26]. An estimation of the spill-out effects introduced was performed both analytically and by means of numerical simulations for realistic internal carotid arteries dimensions and scanner resolution (see Appendixes A.2 and B.1). The estimated spill-out coefficient of 0.51 was used to correct the underestimated image-derived blood input functions.

Time-Activity Curves Extraction and Registration
All TACs within the brain region (excluding skull and peripheral cerebrospinal fluid) were extracted from each registered dynamic volume. For each patient, a brain ROI was first computed as follows: A late frame (20-30 min p.i.) was computed by averaging registered frames 616 to 906. The patient's most recent T2 FLAIR MR image (no more than one month apart) was rigidly registered to the late frame by mutual information maximization [23]. The brain volume was then segmented on the registered T2-FLAIR image using a combination of thresholdings and morphological operations, available in an in-house C++ software based on VTK [27] and ITK [28]. All TACs within the brain ROI were spatially smoothed for noise reduction by averaging within a 3 vox × 3 vox × 3 vox neighborhood. For inter-patient normalization purpose, all smoothed TACs were converted to SUV units and temporally registered to account for inter-patient injection delay. The injection delay was computed for each patient by cross-correlation maximization between the patient's IDIF and the most delayed IDIF observed, used as a reference (patient 17 in Table A1). In total, 6,420,534 TACs were extracted from the 30 min and 30 s [ 11 C]MET PET dynamic scans of our 33 glioma patients.

Biological Tumor Volume Delineation
For each lesion, the biological tumor volume (BTV) was delineated on the average late frame (20-30 min p.i.) using a threshold of 1.6 times the normal brain uptake [29], computed as the mean SUV within a spherical ROI with radius 1 cm placed in the contralateral hemisphere symmetrically to the lesion location with regard to the falx cerebri. Basal nuclei were manually removed when erroneously included in the BTV. All segmentations were reviewed and approved by an experienced nuclear medicine physician.

Pharmacokinetic Modeling
For short (30 min in this work) dynamic brain [ 11 C]MET PET scans, the observed time variations in the voxel activity concentrations are mainly attributed to the tracer transport rather than protein synthesis [30]. The simplest compartmental model describing the transport of the tracer from blood to tissues is the one-tissue compartment model (1TCM) illustrated in Figure 1.
One-tissue compartment model (1TCM) describing the tracer transport from blood to tissues. c p (t) and c 1 (t) are, respectively, the tracer concentration in the blood and tissue compartment and K 1 and k 2 are the associated transport rate constants.
The transport processes involved are classically modeled by first-order kinetics, yielding: where c p (t) and c 1 (t) are, respectively, the tracer activity concentration of the plasma and the tissue compartment at time t, and K 1 and k 2 are the transport rate constants from plasma to tissues and conversely. Calculating the unilateral Laplace transform of Equation (1) and considering zero initial concentration in both compartments leads to the expression of the transfer function H 1 (s) of the 1TCM: where C p (s) and C 1 (s) are the respective Laplace transforms of c p (t) and c 1 (t) for Laplace variable s. In this work, we propose the following change of variables: This transformation is motivated by linear time-invariant system theory where G and τ are commonly referred to as the system static gain and time constant. From a physiological point a view, G is also referred to as the tracer total volume of distribution and τ as the inverse washout coefficient. G thus reflects the steady-state tissue to plasma uptake ratio (i.e., the asymptotic uptake value for a unit step plasma input function) whereas τ describes the uptake dynamics off-steady-state (i.e., the shape of the tissue TAC) and is expected to capture the sought-after 'increasing' or 'decreasing' behavior of tissue TACs. For the sake of clarity, the K 1 /k 2 and 1/k 2 notations will be used hereafter in place of G and τ.
Since a certain fraction of blood vessels is comprised within a single PET voxel (around 5% in healthy brain tissues), the whole voxel TAC is given by: where α is the vascular fraction and c v (t), c b (t), and c t (t) are, respectively, the voxel, blood, and tissue total activity concentration, with c t (t) = c 1 (t) for the 1TCM. Furthermore, c b can be expressed by: where h is the haematocrit (i.e., the erythrocyte volume fraction of blood), c e (t) is the erythrocyte activity concentration, and c p tot (t) is the total plasma activity concentration attributed to the tracer and its metabolites. Since no blood samples were taken from the imaged patients and no population data were available for the erythrocyte uptake of [ 11 C]MET, c e (t) was neglected in first approximation. The haematocrit h was set to its average population value of 0.45. The tracer plasma activity c p (t) was computed from the total plasma activity c p tot (t) using the time-dependent parent fraction f (t): In this work, f (t) was linearly interpolated from population data points reported by Sato and colleagues [31] and depicted in Figure 2.  Time-dependent free [ 11 C]MET plasma fraction used for metabolite correction. Population data points represented by blue triangles with error bars (mean ± standard deviation) were obtained by Sato and colleagues from 18 glioma patients [31]. The least-squares fitted linear approximation is plotted in red.
Taking into account the amount of time required for the tracer to flow from the carotid arteries-where the blood input function is extracted-to the voxel location, the voxel TAC is finally given by: where d is the system delay and h 1 is the inverse Laplace transform of H 1 . The 1TCM kinetic parameter K 1 and k 2 as well as the vascular fraction α and the delay d were individually estimated for each voxel TAC in our dataset (33 patients, 6,420,534 curves) by least-squares transfer function fitting using the SciPy's 'optimize' and 'signal' modules [32] in Python. The added value of overlapping frames over adjacent frames for 1TCM kinetic parameters fitting was investigated by means of numerical simulations (see Appendixes A.4 and B.2).

Numerical Simulations
Numerical simulations were conducted to assess the impact of data noise on the ability of PCA to capture characteristic TAC behaviors. Synthetic time activity curves were generated using the 1TCM and the tri-exponential blood input function model proposed by Feng and colleagues [33]: To be able to compute an analytical solution for the tissue tracer activity, a linear approximation of the time dependent parent fraction f (t) = m t + p was built, leading to the following expression for the plasma input function: The value of p was fixed such that free [ 11 C]MET plasma fraction is equal to 1 at injection time and the value of m was least-squares fitted on the experimental data of Sato and colleagues [31]. The resulting linear approximation is depicted in Figure 2 (red line).
Taking into account the carotid-to-voxel delay d, an analytical expression of the 1TCM tissue tracer activity c t (t) for t ≥ d and the plasma input function in Equation (9) is derived in Equation (A16) (see Appendix A. 3). An analytical expression for the mean voxel activity concentration c t s →t e v for a frame starting at time t s and ending at time t e is then derived in Equations (A17)-(A19).
For each simulated frame c t s →t e v computed using Equation (A17), synthetic noise was generated using the model proposed by Logan and colleagues [34]: where c t s →t e v,noisy is the noisy mean voxel tracer activity concentration for a frame starting at time t s and ending at time t e , β is a scaling factor, λ is the isotope decay constant, and G(0, 1) is a random number generated from a Gaussian distribution with mean 0 and standard deviation 1.
Synthetic TACs were generated using Equations (A17)-(A19) from the 1TCM kinetic parameter values previously fitted on each real TAC in our dataset (33 patients, 6,420,534 TACs). The input function parameter values used in Equation (8) were leastsquares fitted on the corresponding patient's IDIF using SciPy's 'optimize' and 'signal' modules [32] in Python. The generated synthetic dataset is thus similar to our real dataset but has known ground truth signal. Synthetic noise was added to simulated TACs using Equation (10) for increasing noise levels β ranging from 0.25 to 2.0 with step 0.25.

Principal Component Analysis
Principal component analysis is a commonly used dimension reduction technique that aims at finding the linear transformation of the initial data space into the so-called 'components' space which best explains the data variance. More formally, let X be the data matrix of dimensions n × p where n is the number of observations (i.e., the number of TACs) and p the number of variables (i.e., the number of samples per TAC). The method aims at finding the new axis system e 1 , e 2 , . . . , e p in which the variance is maximized under constraints e i e j = δ i,j and s 2 e 1 ≥ s 2 e 2 ≥ · · · ≥ s 2 e p , with δ i,j being the Kronecker delta and S the data covariance matrix given by: where X c is the centered data matrix, i.e., the data matrix from which variable means have been subtracted columnwise. It can be shown that such components are given by the eigenvectors of S ordered by decreasing eigenvalues, corresponding to the respective component variances. As the amount of explained data variance decreases with the component number, it is expected that a sufficiently high amount of data variance can be explained by the first m p components, which is the case for highly correlated initial variables. Our real and synthetic TAC datasets were split into PCA model-building and evaluation sets on a patient basis in order to evaluate the inter-patient generalization performance of the PCA model built. The model-building and evaluation sets comprise TACs from 20 (patients 1 to 20 in Table A1) and 13 (patients 21 to 33 in Table A1) patients, respectively, totaling 3,974,466 and 2,446,068 TACs.
The impact of noise on the true signal reconstruction and on the definition of the first 6 principal components was first investigated using the synthetic TAC dataset. For each level of noise, the principal components were first determined on the noisy model-building TAC set (20 patients, 3,974,466 TACs). The component values were then computed for each noisy evaluation TAC (13 patients, 2,446,068 TACs) at the considered noise level. Estimated denoised TACs were finally reconstructed from the first 4 and 6 component values and compared to the corresponding true unnoisy TAC using the mean squared error. To assess the robustness of our approach in the presence of atypical kinetics, the same numerical experiment was conducted by considering tumor evaluation TACs only (i.e., the 39,780 TACs extracted from the delineated BTVs). The influence of the noise level on the component definition was also assessed by computing the cosine between each of the first 6 principal components determined on the noisy model-building dataset and the corresponding component obtained for a noise level of zero.
The first 6 principal components were then determined on the real model-building TAC set (20 patients, 3,974,466 TACs) and compared to these determined on the synthetic model-building dataset for a noise level of zero by computing their respective cosine. The component values were then computed for each TAC of the 13 evaluation patients and mapped spatially to their respective voxel location. The possible dependence between the minimum, maximum, and mean values of the resulting PC maps within the BTV versus the 1p/19q codeletion and IDH mutation status was investigated by means of Mann-Whitney U tests with a significance level set to 0.05. The PC maps were finally compared to their homologous 1TCM kinetic parameter maps, visually and by means of Spearman's correlation analyses for the 1TCM parametric maps.
Due to the limited amount of memory available, the incremental PCA (IPCA) algorithm was used for all analyses instead of classical matrix decomposition since it allows a batch processing of the data. The implementation used is available in Python as part of the scikit-learn's 'decomposition' module [35]. No noise normalization was performed prior to PCA in this work.

Synthetic Data
The first six principal components (PCs) determined on the synthetic model-building dataset (20 patients, 3,974,466 TACs) for a noise level of zero are depicted in Figure 3. Their corresponding explained variance ratios are, respectively, 90.37%, 8.45%, 0.71%, 0.31%, 0.10%, and 0.03% for PC1 to 6, totaling 99.97% of the explained variance. It should be noted that PC5 and PC6 have very low contribution to the explained variance. A clinical interpretation of the PCs can be proposed based on their respective shape. Aside from its early peak, PC1 ( Figure 3a) assigns a relatively constant weight to the later samples and thus partly reflects the mean tracer uptake. PC2 (Figure 3b) overweights the early TAC samples, then assigns a rapidly decaying weight to the next samples with a slightly negative value for the last samples. This component thus has a high value for TACs with a high early activity as observed in blood after bolus injection of the tracer. PC3 (Figure 3c) assigns a negative weight to the first half of the samples and a positive weight to its second half, with a quasi-linear increase. It thus has a negative value for decreasing TACs, a positive value for increasing TACs and a value around zero for flat TACs. PC4 ( Figure 3d) assigns a strongly positive weight to the very first samples, a strongly negative weight to the next few samples then a small quasi-constant weight to the last samples. PC4 thus assigns a negative value to TACs with delayed arterial peak and a positive value to TACs with non-delayed arterial peak. PC5 and 6 (Figure 3e,f) are less easily interpreted but together only explain 0.13% of the variance.   Figure 4c depicts the impact of noise on the definition of the first 6 principal components, assessed by their cosine with the respective components determined for a noise level of zero (see Figure 3). The PC5 cosine progressively drops from a noise level of 1.0 whereas a rapid drop of the PC6 cosine is observed for noise levels above 0.5. PC1 to 4, on the other hand, remain stable even for high noise levels.  Examples of true, noisy, and PCA-denoised evaluation TACs with 4 components are depicted in Figure 5. Reconstruction with only 4 components is remarkably accurate while efficiently removing noise.

Real Data
The first six principal components determined on the real model-building dataset are depicted in Figure 6. Their corresponding explained variance ratios are, respectively, 72.27%, 7.65%, 0.53%, 0.39%, 0.33%, and 0.30% for PC1 to 6, totaling 81.47% of the explained variance. Strong similarities are, respectively, observed between PC1-4 of the synthetic dataset (see Figure 3) and PC1-3 and 5 of the real dataset (see Figure 6), hence the same clinical interpretations can be made for the latter. These similarities are confirmed by their respective cosines of 0.9987, 0.9891, 0.9383, and 0.9111. The residual differences are suspected to originate from (i) the limited ability of the 1TCM used for data synthesis to fully capture the whole range of observed TAC dynamic behaviors and (ii) unmodeled additional sources of noise and artifacts related to data reconstruction and residual patient motion. It should be noted that PC4-6 of the real dataset have very similar explained variance ratios, hence their order was not considered informative as they could have been reversed for another similar dataset. Since PC4 and 6 of the real dataset are hardly interpretable and do not match any of the first six PCs of the synthetic dataset-as opposed to PC5-they will not be further considered from now on.

Parametric Maps
Hereafter, the PC definitions used are these determined on the real model-building dataset (see Figure 6). The parametric maps obtained by voxelwise mapping of the PC values and of the 1TCM kinetic parameters computed over the voxel TACs of a glioblastoma evaluation patient (patient 21 in Table A1)  Furthermore, the PC3 map ( Figure 7g) clearly highlights regions with characteristic dynamic behaviors, as illustrated by TACs A, B, and C. Voxels with fast increasing then decreasing TACs appear brighter on the PC3 map (red arrow), whereas progressively increasing TACs appear darker (green arrow). Voxels with relatively flat TACs appear in medium gray value on the PC3 map (blue arrow). In contrast, the homologous 1/k 2 map fails to clearly highlight these behaviors, as illustrated by voxels pointed by the green and blue arrows, both appearing darker in Figure 7k.  Static PET images (20-27 min p.i.), PC3, and 1/k 2 maps along with representative smoothed TACs at voxels pointed by the red, and green arrows are depicted in Figure 8 for 4 additional evaluation patients (patients 23, 25, 30, and 32 in Table A1). As opposed to the static PET image, the PC3 map allows to distinguish voxels with 'decreasing' or 'flat' TACs (red) from voxels with 'flat' or 'increasing' TACs (green). The homologous 1/k 2 maps exhibits patterns similar to the PC3 maps within the BTV but with a substantially lower contrast. Similar behaviors were observed for the other eight evaluation patients. Interestingly, the substantial late uptake increase of the green TAC from patient 23 in Figure 8 is not totally captured by PC1-3 and 5 (see black curve) but PC4 and 6 are also required for a more accurate reconstruction.    Table A1) along with the smoothed TACs at voxels pointed by the red (4th row) and green (5th row) arrows, respectively. Red TACs have a 'decreasing' or 'flat' behavior whereas green TACs have a 'flat' or 'increasing' behavior. TAC reconstruction with 4 of the 5 principal components (PC1-3 and 5) is superimposed in black to each TAC.
The biological tumor volume (BTV), mean and maximum tumor-to-background ratio (TBR) evaluated on the static PET image (20-27 min p.i.) as well as tumor contrast (see definition in Equation (A20)) values evaluated on the static PET image (C static ) and on the PC1 map (C PC1 ) are reported for each lesion in Table A3. The corresponding distributions are summarized by means of boxplots in Figure A3a-e for lesions with a non-zero BTV. The Bland-Altman plot of C PC1 versus C static is depicted in Figure A3f, illustrating the systematic higher tumor contrast of the PC1 maps compared to the static PET image.
The minimum, maximum, and mean values of the 1TCM kinetic parameter maps and of the PC maps within the BTV are, respectively, reported for each lesion in Tables A4 and A5. The corresponding distributions are summarized by means of boxplots in Figures A4 and A5 for lesions with a non-zero BTV. None of the performed Mann-Whitney U tests investigating the dependence of the minimum, maximum, and mean PC values within the BTV were found statistically significant. Nevertheless, 1p/19q non-codeleted (p = 0.15) and IDH wildtype (p = 0.32) tumors tend to exhibit lower minimum PC3 values, as suggested by the grouped boxplots in Figure A6. It should be noted that removing the outlier (see asterisk mark in Figure A6a)-identified by a Grubb's test (p < 0.05) performed after verification of the data normality by a Shapiro-Wilk test (p = 0.98)-leads to statistical significance for PC3 min versus the 1p/19q codeletion status (p = 0.04).

Parametric Map Correlations
The similarities between PC and 1TCM kinetic parameter maps observed are confirmed by the pairwise Spearman's correlation coefficients computed voxelwise for each of the 13 evaluation patients taken individually and summarized in Table 1. Feature pairs PC2|α and PC5|d strongly correlate whereas PC1|K 1 /k 2 and PC3|1/k 2 very strongly correlate. Moderate positive correlations are also found between PC1|α-imputed to the overweighting of the very first samples of PC1 (see Figure 6a)-and PC3|K 1 /k 2 . All associated p-values were <0.05-and <0.0001 for most of them due to the large number of analyzed voxels (∼200,000 voxels per patient scan)-except for pair PC5|K 1 /k 2 of evaluation patient 33.

Discussion
We showed the ability of PCA to accurately capture characteristic dynamic behaviors from a broad spectrum of high dimensional [ 11 C]MET PET TACs extracted from the whole brain region of 20 patients, most of whom have been treated multiple times. By means of realistic numerical simulations, we demonstrated the robustness of our approach to noise and validated the accuracy of the unnoisy signal reconstruction on a separate synthetic evaluation dataset not used for PCA model building, without the need for prior noise normalization. A possible explanation for this unexpected robustness to noise is that only TACs within the brain region were considered for this work-as opposed to the works of Pedersen and colleagues [20] and Šámal and colleagues [21]-hence background noise had no influence on the computed components.
Four components (PC1-3 and 5) among the six that best explain variations observed among the real whole brain TACs of 20 patients respectively matched the first four principal components derived from our synthetic model-building dataset. These four components have been found to be of clinical interest. PC1 provides a contrast similar to the routinely acquired static PET image (20-27 min p.i.) but has two advantages over the latter: (i) the PC1 map has a higher tumor contrast (see definition in Equation (A20)) than the static PET image, as shown in Figure A3f and (ii) PC1 does not depend on an arbitrary chosen imaging timing since each of the acquired samples starting from tracer injection contributes to some extent to the component value, making it more robust to inter-protocol variations. PC2 reflects the voxel vascular fraction but its small negative weighting of the late samples (see Figure 6b) makes interpretation less straightforward in tumor regions where high late uptake negatively influences the component value. PC3 provides an interesting dynamic contrast related to the shape of the tissue TAC, distinguishing voxel with fast increasing then decreasing, flat, or progressively increasing TACs. PC5 appears to be related to the relative blood delay and could reflect the effectiveness of micro-vascularization, which is known to be impaired in gliomas [36].
These four PCs strongly to very strongly correlate with the kinetic parameters of the widely used 1TCM PK model, which partly confirms our intuition concerning their clinical interpretation. On the other hand, the correlation between the kinetic parameters and PCs highlights the ability of a simplistic model such as the 1TCM to accurately capture the large variability observed among whole brain TACs of glioma patients. PCA has, however, three advantages over the 1TCM. First, it does not require a blood input function, avoiding invasive arterial sampling procedures or manual extraction of an IDIF potentially affected by PVE and tracer metabolites. Second, PCA model building is a parameter-free procedure, which makes it robust to inter-protocol variations. Third, the PCA parametric maps are much less computationally expensive to compute, with processing times of a few seconds for any new dynamic volume versus more than 3 h of fit for the 1TCM on a high-end computer. Alternatively, a basis function approach could be used for 1TCM parameters fitting to shorten processing times and improve stability [37] but at the expense of the precision on the kinetic parameter values and still longer processing times than the calculation of the component values. Furthermore, the PC3 map has been shown to outperform the homologous 1TCM 1/k 2 map in differentiating voxels with 'increasing' and 'decreasing' TAC, which has been previously linked to the tumor aggressiveness as will be further discussed.
It should however be noted that the PCs must not been seen as a surrogate for the kinetic parameters since each PC may reflect a mixture of kinetic parameters (see off-diagonal correlations in Table 1) and of other factors not captured by the 1TCM. PCs instead capture the largest variability factors observed among TACs resulting from complex underlying biological processes, under an orthogonality constraint. An unequivocal link between one PC and one biologically relevant micro-parameter is not guaranteed eithernor is it the case for the kinetic parameters-which both capture complex amino acid transport processes dependent, among others, on the number of amino acid transporters along the endothelial and tumor cell membrane, the concentrations of all endogenous transporter shared substrates in every model compartment, and the relative sizes of the compartments and thus (over-) cellularity [38].
Our findings may nevertheless have important implications for clinical management of gliomas: PCA applied to dynamic [ 11 C]MET PET data is likely to provide additional quantitative spatial information on glioma heterogeneity. In particular, the PC3 map provides a novel contrast complementary to the static PET image that can help distinguishing voxels with similar late uptake values but different uptake time courses, as illustrated in Figure 8. This component translates the 'increasing' and 'decreasing' behavior of TACs observed of Pöpperl and colleagues at the whole tumor level [8] into a quantitative voxelwise metric. This metric may partly reflect biologically relevant parameters such as local overexpression of LAT1-an important amino acid transporter over-expressed in gliomas at both the blood-brain barrier and tumor cell level [39,40] transporting both [ 11 C]MET [39,41] and [ 18 F]FET [41]. LAT1 being an obligatory exchanger, its over-expression may indeed be associated with an increase in both the influx and efflux of the tracer [41], leading to faster uptake dynamics that more closely follows the arterial input signal, hence characterized by a lower PC3 value. If this hypothesis turned out to be verified, application of the proposed methodology to dynamic [ 18 F]FET data would also be of interest with an even more pronounced effect expected since [ 18 F]FET is more specifically transported by LAT1 [42] and is not incorporated into proteins [15,41]. However, the possible relation between PC3 and the expression of LAT1 still needs to be verified, e.g., by means of targeted biopsies. Although not statistically significant, our preliminary analysis results are encouraging as they suggest that 1p/19q non-codeleted and IDH wildtype tumors may be characterized by lower minimum PC3 values (see Figure A6) and would thus exhibit more 'decreasing' TACs. These findings are in accordance with the previous work of Vetterman and colleagues on dynamic [ 18 F]FET PET in gliomas, which reported a shorter time-topeak in IDH wild type tumors [13]. However, relations between the patient outcome or glioma molecular features such as the IDH mutation or 1p/19q codeletion status and the distribution of the PC values within the BTV should be further investigated as part of a larger scale study with stricter inclusion criteria in terms of undergone treatments and subgroup balance, but it was out of the scope of this preliminary methodological work. To the best of our knowledge, this is the first work to promote the added value of dynamic PET acquisitions with [ 11 C]MET in glioma patients to assess intra-tumor heterogeneity, as previously demonstrated for [ 18 From a clinical point of view, our approach has the advantage of requiring only little modifications of the routine protocol, that is a longer acquisition time on the PET/CT tomograph but identical total duration of the procedure for the patient and no arterial sampling. Instead of injecting the patient outside the scanner and wait for a post injection delay before the acquisition, our preliminary results indeed support that additional information on the local tracer uptake dynamics can be obtained by injecting the patient inside the scanner right after the acquisition start without further modification of the protocol. The main drawback of the approach is related to the longer scanner occupancy, which would nevertheless have little impact if the procedure is reserved to the first patient of the day, selected accordingly. Besides, scanner occupancy will tend to be less problematic in the future with the advent of hybrid PET/MR scanners, which offer more time for PET acquisition as the MR acquisition typically constrains the total duration of the procedure. Once the acquisition is completed, the reconstruction and computation of the PC maps from a pre-built PCA model can be fully automatized as it does not require any manual processing and the results could be sent to the archiving system along with the static PET image. Although beyond the scope of this methodological work, potential clinical applications of the derived PC maps are multiple. The maps could first be introduced as an additional source of information for the definition of the resection margins or dose planning in radiotherapy-areas with an atypical dynamic behavior being potentially preferable to include in the clinical target volume. Additionally, the maps could be used for biopsy targeting in heterogeneous tumors such as gliomas with brighter or darker foci being possibly characterized by more aggressive histological features, reducing the risk of tumor under-grading. The analysis of staggered PC maps could also provide information on the treatment response over time. Finally, the PC maps may potentially provide additional information to distinguish glioma recurrence from radio-necrosis post radiotherapy.
The pharmacokinetic analyses conducted as part of this work were however prone to several limitations. Indeed, the decision was made for this study not to perform arterial sampling which requires a dedicated device and poses risks for the patient and the nursing staff. Instead, the blood input functions used for the kinetic analyses were extracted from the image within the petrous segment of the internal carotid arteries. As a consequence, the input functions used were affected by PVE and tracer metabolite activity, which could not be assessed individually for each patient. Population-based haematocrit and metabolite corrections as well as model-based spill-out coefficient estimation were thus proposed to counteract these effects. Moreover, spill-in correction of the IDIFs could not be trivially performed and is still a challenging open problem in PET imaging [43], which was out of the scope of this work. Nevertheless, it turned out in the course of this study that the introduction of the proposed input function corrections only affected the absolute values of the fitted kinetic parameters but that the contrast of the derived kinetic parametric maps was preserved. Since we were interested in this study in the relative quantification between voxels to highlight intra-tumor heterogeneity of the uptake dynamics, residual uncertainties in the corrected image-derived input functions were not suspected to significantly impact the conclusions of this work regarding the superiority of PCA to this extent. More complex pharmacokinetic models such as the two-tissue compartment model (2TCM) may also be considered to account for the tracer transport from the extra-to the intracellular space [44]. For longer scans, the synthesis of tracer metabolites in tissues-including proteins-would also need to be incorporated into the model as previously proposed for the assessment of the cerebral protein synthesis rate with L-[1-11 C]leucine [45]. Nevertheless, as models grow in complexity, voxelwise fitting becomes even more challenging and prone to severe instabilities, as reported for the 2TCM [18]. In contrast, a model-free approach such as PCA has the advantage of being only concerned with extracting the highest variation factors of TACs, which could be as clinically relevant as model parameters if they turned out to be correlated with clinical or histological data.
For the sake of completeness, it should be noted that other dimension reduction techniques have been proposed in the literature for the analysis of dynamic PET data, such as factor analysis of dynamic structures (FADS) and non-negative matrix factorization (NMF) [46]. These are distinguished by their ability to isolate signal from noise, which strongly depends on the assumptions made on the noise distribution in PET images. Whereas the noise related to the event counting in PET imaging follows a Poisson distribution, the noise distribution in reconstructed PET images is less well characterized due to alterations related to the system hardware and reconstruction algorithm-including scatter and attenuation corrections [46]-hence remains an open problem. Comparison of these methods in the particular case of dynamic [ 11 C]MET PET data of glioma patients would also be of interest as a future work.

Conclusions
We showed the ability of principal component analysis to extract meaningful parametric maps from noisy high dimensional dynamic [S-methyl-11 C]methionine PET scans of glioma patients with little modification of the routine protocol. One of these maps was found to reflect at the voxel level the previously reported 'increasing' or 'decreasing' behavior of TACs within the tumor, which could potentially be linked to the aggressiveness heterogeneity within the tumor. Such maps could be of great interest for tumor characterization as well as for surgery and radiotherapy planning in addition to conventional static PET imaging. This early methodological work paves the way for many possible clinical studies.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on reasonable request from the corresponding author.

Acknowledgments:
The authors would like to thank the technical, medical, and scientific staff in charge of PET tracer synthesis and image acquisition for the needs of this work.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

. Spill-Out Estimation
The spill-out estimation for IDIF correction was performed both analytically and by means of computer simulations. For the sake of simplicity, the scanner resolution was considered spatially constant and isotropic, characterized by a full-width at half-maximum (FWHM) of 4.5 mm consistent with our scanner specification. The scanner PSF g was modeled by a 3D isotropic Gaussian with standard deviation σ = FWHM/ 2 2 ln(2) ≈ 1.91 mm. The petrous segment of the carotid artery from which the whole blood input function is extracted in the image was modeled by a cylinder with radius R = 2.4 mm [47] and height H = 30 mm, determined experimentally based on time-of-flight MR data. For the calculation of the spill-out coefficient, a homogeneous activity distribution f (x, y, z) was constructed for such a cylinder centered at the origin, given by: The convolution c of f with the scanner PSF g at point (x, y, z) is then given by: To solve Equation (A4), the problem was re-expressed in cylindrical coordinates using the following change of variables: leading to: Since no generic analytical solution exists for Equation (A8) at every point (x, y, z), we computed the spill-out coefficient at the central point (0, 0, 0) of the cylinder where it is supposed to be maximum, which gives: Provided that an analytical expression can be found for Equation (A8), the average spill-out coefficient c voxel for a voxel of size 2 mm × 2 mm × 2 mm at the center of the petrous segment of the internal carotid would have been obtained by: Alternatively, a numerical approximation of c voxel was computed by convolution of a 3D image with voxel size 0.02 mm × 0.02 mm × 0.02 mm in the center of which a cylinder of radius R and height H was drawn with a 3D isotropic Gaussian kernel with standard deviation σ/0.02 pixels truncated at 4 σ. An approximated value of c voxel was obtained by averaging the 100 × 100 × 100 central voxels in the convolved image.

Appendix A.3. Numerical Simulations
Taking into account the carotid-to-voxel delay d, the analytical expression of the 1TCM tissue tracer activity c t (t) for t ≥ d and the plasma input function in Equation (9) is given by: The mean voxel activity concentration c t s →t e v for a frame starting at time t s and ending at time t e is then given by: where c b (t) and c t (t) are, respectively, given by Equations (A18) and (A19):

. Overlapping Frames
To quantify the added value of overlapping frames over distinct adjacent frames for the estimation of the 1TCM kinetic parameters, the following numerical experiments were conducted. The Feng's input function parameters in Equation (8) were first leastsquares fitted for each of the 33 patients' IDIF using the SciPy's 'optimize' and 'signal' modules [32] in Python. Frame-averaged blood input functions were then computed for each of the 33 fitted Feng's parameter values using Equation (A18) and for two distinct framing strategies: (i) 906 overlapping frames of length 20 s with a shift of 2 s (the framing used in this work) and (ii) 92 non-overlapping adjacent frames of length 20 s. For the same two framing strategies, frame-averaged voxel TACs were generated using Equations (A17)-(A19) for 10,000 combinations of kinetic parameter values randomly chosen among the whole set of kinetic parameter values previously fitted on the real patients TACs (33 patients, 6,420,534 curves) with the corresponding Feng's parameter values. For each noise level in range 0.0 to 2.0 with step 0.5 and for both framing strategies, synthetic noise was added to the frame-averaged blood input functions and to each of the 10,000 frame-averaged voxel TACs using Equation (10). Kinetic parameter values were finally estimated for each pair of noisy frame-averaged blood input function and voxel TAC using our least-squares transfer function fitting routine in Python and the relative errors between ground truth and estimated kinetic parameter values were computed. To assess the robustness of our approach in the presence of atypical kinetics, the same numerical experiment was conducted by randomly choosing the 10,000 combinations of kinetic parameter values among these extracted from the BTVs only.

Appendix B.1. Spill-Out Estimation
The evaluation of Equation (A14) for respective values of R and H of 2.4 mm and 30 mm gave a maximum spill-out coefficient c(0, 0, 0) ≈ 0.55. Numerical convolution results are depicted in Figure A1. The spill-out coefficient value of the central 0.02 mm × 0.02 mm × 0.02 mm simulation grid voxel was of 0.55, in accordance with the analytical value. The average spill-out coefficient value within the 100 × 100 × 100 central voxels of the simulation grid-corresponding to the central PET image voxel-was of 0.51. This value was used for IDIF correction prior to the PK analyses.

Appendix B.2. Overlapping Frames
The median relative errors on the fitted 1TCM parameter values computed over the 10,000 simulated whole brain TACs are depicted in Figure A2 for each level of noise and both framing strategies. Except for the delay parameter d, overlapping frames systematically provide lower median relative errors. Furthermore, the reduction of the fitting error for overlapping compared to adjacent frames increases with the level of noise. The decrease in median error using overlapping frames for a realistic noise level of 1.5 was of 5.9%, 13.1%, and 12.9% for K 1 , k 2 , and α, respectively. The associated p-values returned by the Wilcoxon's signed-rank test were not considered reliable due to the large number of samples (n = 10, 000), hence the Wilcoxon's effect sizes are reported instead in Table A2. Similar curves were obtained when considering kinetics extracted from the BTVs only, with an overall lower median error for both strategies and all noise levels. The corresponding absolute decrease in the median error was of 2.7%, 7.3%, and 10.5% for K 1 , k 2 , and α at a noise level of 1.5, which confirms the robustness of our approach in the presence of atypical kinetics.   Table A3 reports for each lesion the biological tumor volume (BTV), the mean and maximum tumor-to-background ratio (TBR) evaluated on the static PET image (20-27 min p.i.) as well as the tumor contrast evaluated on the static PET image (C static ) and on the PC1 (C PC1 ) map, given by: where tumor mean and contralateral mean are, respectively, the mean value within the BTV and the contralateral spherical ROI (see Section 2.5). The corresponding distributions are summarized by means of boxplots in Figure A3a-e for lesions with a non-zero BTV. The Bland-Altman plot of C PC1 versus C static is depicted in Figure A3f, illustrating the systematic higher tumor contrast of the PC1 maps compared to the static PET image. The minimum, maximum, and mean values of the 1TCM kinetic parameters and of the PC values within the BTV are, respectively, reported for each lesion in Tables A4 and A5. The corresponding distributions are summarized by means of boxplots in Figures A4 and A5 for lesions with a non-zero BTV. Figure A6 depicts boxplots of the minimum PC3 values within the BTV grouped by 1p/19q codeletion status ( Figure A6a) and IDH mutation status ( Figure A6b).      Table A3 for lesions with a non-zero BTV summarized by means of boxplots (horizontal line: median, box: interquartile range, whiskers: ±1.5 interquartile range, *: outliers) (a-e) and Bland-Altman plot of C PC1 versus C static (f).  Table A4 for lesions with a non-zero BTV summarized by means of boxplots (horizontal line: median, box: interquartile range, whiskers: ±1.5 interquartile range, *: outliers). (d) Figure A5. Distribution of the minimum, maximum, and mean PC values within the BTV in Table A5 for lesions with a non-zero BTV summarized by means of boxplots (horizontal line: median, box: interquartile range, whiskers: ±1.5 interquartile range, *: outliers).

Appendix C. Supplementary Discussion
We showed that the reconstruction of a large number of overlapping frames compared to a reduced number of adjacent frames of the same length benefits the estimation of the 1TCM kinetic parameters increasingly with the noise level, except for the blood delay d. Possible explanations for this gain are that: (i) overlapping frames provide a better sampling of the arterial peak for both blood and tissue TACs and (ii) an increased number of correlated data points helps the fitting algorithm to overcome noise. The latter remark does however not appear to hold for the estimation of the blood delay d, which consists of an arterial peak alignment process that could be negatively impacted by an increased number of early noisy TAC samples. Nevertheless, the blood delay d is rarely considered as a variable of interest in pharmacokinetic studies but is rather used to provide more flexibility to the fitting algorithm. The overlapping framing approach proposed has however two major drawbacks, that are (i) a longer reconstruction time and (ii) a larger memory usage for storage and processing. It should finally be noted that the overlapping framing used in this work may not be optimal for kinetic parameter values estimation. Determining the best dynamic PET framing strategy for pharmacokinetic studies would be of interest but is out of the scope of this paper.
We proposed to average the TACs of the 10 brightest voxels within the petrous segment of both internal carotid arteries on an early PET frame for blood input function extraction (see Section 2.3) since they were expected to be the least affected by PVE. It turned out in practice that these voxels are indeed almost systematically located along the central line of the segment, which is consistent with our simulations since the central line has the largest theoretical spill-out coefficient value. Whereas this method has the great advantage of not requiring an invasive arterial sampling, we showed hereabove that the activity within the central voxels is systematically underestimated by a factor near 2 due to partial volume effects resulting from the limited resolution of the scanner and the restricted diameter of the internal carotid arteries. We proposed to address this issue by scaling the extracted blood input functions by a factor 1/0.51 ≈ 1.96. This correction remains however limited since the geometry of the carotid arteries varies among patients [47]. Additionally, investigating the effects of these variations on the kinetic parameters estimation would also be of interest. Deriving a spill-in coefficient, in contrast, is less trivial since all TACs within the reconstructed volume are expected to impact the input function.