Data-driven, voxel-based analysis of brain PET images: Application of PCA and LASSO methods to visualize and quantify patterns of neurodegeneration

Spatial patterns of radiotracer binding in positron emission tomography (PET) images may convey information related to the disease topology. However, this information is not captured by the standard PET image analysis that quantifies the mean radiotracer uptake within a region of interest (ROI). On the other hand, spatial analyses that use more advanced radiomic features may be difficult to interpret. Here we propose an alternative data-driven, voxel-based approach to spatial pattern analysis in brain PET, which can be easily interpreted. We apply principal component analysis (PCA) to identify voxel covariance patterns, and optimally combine several patterns using the least absolute shrinkage and selection operator (LASSO). The resulting models predict clinical disease metrics from raw voxel values, allowing for inclusion of clinical covariates. The analysis is performed on high-resolution PET images from healthy controls and subjects affected by Parkinson’s disease (PD), acquired with a pre-synaptic and a post-synaptic dopaminergic PET tracer. We demonstrate that PCA identifies robust and tracer-specific binding patterns in sub-cortical brain structures; the patterns evolve as a function of disease progression. Principal component LASSO (PC-LASSO) models of clinical disease metrics achieve higher predictive accuracy compared to the mean tracer binding ratio (BR) alone: the cross-validated test mean squared error of adjusted disease duration (motor impairment score) was 16.3 ± 0.17 years2 (9.7 ± 0.15) with mean BR, versus 14.4 ± 0.18 years2 (8.9 ± 0.16) with PC-LASSO. We interpret the best-performing PC-LASSO models in the spatial sense and discuss them with reference to the PD pathology and somatotopic organization of the striatum. PC-LASSO is thus shown to be a useful method to analyze clinically-relevant tracer binding patterns, and to construct interpretable, imaging-based predictive models of clinical metrics.

Introduction Pathological processes associated with neurological disorders often develop in distinct spatiotemporal patterns. These patterns can be imaged with positron emission tomography (PET) or single-photon emission computed tomography (SPECT) using appropriate radiotracers. However, traditional quantitative PET and SPECT image analysis metrics, such as the standardized uptake value (SUV) and non-displaceable binding potential (BP ND ), are often evaluated as averages over a specific region of interest (ROI). This approach neglects the spatial distribution of radiotracer binding, which may be affected by disease within the ROIs. There is thus a growing realization that better methods of spatial image analysis are needed to achieve a more complete disease characterization and to improve prediction and tracking of disease progression [1][2][3]. The primary approach to spatial analysis in the context of these objectives has so far relied on the use of texture-and shape-based image features [4][5][6][7][8], leaning heavily on the emerging field of radiomics [9][10][11]. Following a general radiomics approach, a large number (between tens and hundreds) of different features are computed from the images within predefined ROIs. The features are then used as inputs to one of the established machine learning models (neural nets, decision forests, etc.) to predict clinical measures of interest [12][13][14].
However, investigations that use radiomic features are limited in two aspects. First, PET and SPECT brain imaging studies often have a relatively small sample size (10 to 50 subjects) compared to the number of tested features ("large p, small n" problem) [14][15][16][17]. The values of features may likewise depend on internal variables required to define the features, image reconstruction algorithm and ROI definition criteria [7,[18][19][20][21]. This aspect may result in a very large model-parameter search space and substantially increase the risk of model overfitting, particularly if the employed statistical testing strategy is not robust or valid [22][23][24][25]. Second, a relatively high complexity and locality of the radiomic features hinders visualization and biological interpretation of the disease-relevant spatial patterns [26].
In the present work we explore a data-driven approach to the analysis of PET brain images that aims to overcome the shortcomings described above, i.e. overfitting, optimal feature search, and interpretability. In contrast to relying on generic radiomic features, our method identifies disease-and structure-specific spatial patterns of radiotracer binding on a voxel level, and constructs image metrics optimized for prediction of clinical measures. The method uses individual-subject voxel values as inputs, and consists of two steps. In the first step, the high dimensionality of the voxel data is reduced in a blind (i.e. independent of the clinical metrics) manner using principal component analysis (PCA). In the second step, several principal components (PCs) with the highest percent of explained total variance are selected and linearly combined to predict clinical metrics of interest. The combinations are determined using regularized fitting by the least absolute shrinkage and selection operator (LASSO) [27]. In comparison to ordinary multivariate fitting, LASSO selects only those input variables that improve prediction accuracy; inputs that do not contribute to the accuracy are set to zero. We refer to this method as the principal component LASSO (PC-LASSO), and to the best of our knowledge, it has not been previously used in the analysis of PET images. The use of sequential dimensionality reduction and variable selection techniques, wherein PCA is blind and LASSO is used in conjunction with model cross-validation, reduces the risk of overfitting, and enables the use of raw data from 10 3 -10 4 voxels to build predictive models from a much lower number (10-100) of subjects. The fitted PC-LASSO models can be used to compute PC-LASSO estimators, which capture information in a non-local manner by assigning a weight to each voxel that quantifies the voxel's relevance to the predicted clinical metric. The weights can be visualized in the context of the topology of the imaged structure, providing direct information on the spatial characteristics and evolution of disease.
We test the PC-LASSO method by applying it to high-resolution PET images of subjects suffering from Parkinson's disease (PD), the second most common progressive neurodegenerative disorder after Alzheimer's disease [28]. PD is associated with degeneration of presynaptic dopaminergic terminals in the basal ganglia, while initially inducing a compensatory response in the postsynaptic neurons [29,30]. Accordingly, the method is applied to striatal PET images obtained with a pre-synaptic and a post-synaptic dopaminergic tracer: 11 Cdihydrotetrabenazine (DTBZ), a marker of the vesicular monoamine transporter type 2 (VMAT2) that reflects the density of dopaminergic terminals, and 11 C-raclopride (RAC), a dopamine D2 receptor antagonist. We tested whether PC-LASSO can a) identify the spatial progression patterns for DTBZ and RAC, which are expected to be different, and b) predict clinical metrics from the imaging data better than the mean tracer binding ratio (BR), thus demonstrating the clinical relevance of the captured patterns.
We first examine and compare the DTBZ and RAC binding patterns captured by PCA. The robustness and temporal evolution of the patterns are assessed by comparing them between the less-and more-affected sides of the brain, as well as between images obtained from subjects at different disease stages. We then combine several PCs using LASSO to construct crossvalidated models of clinical disease metrics: disease duration (DD) and motor impairment scores. We compare the fit quality and prediction accuracy between the PC-LASSO models and the mean BR within an ROI. Finally, we analyze the distribution of voxel weights in the obtained PC-LASSO estimators in the context of PD pathology.

Subjects and scans
The study included 10 healthy control (HC) subjects and 41 PD subjects combined from different clinical studies, but imaged using the same protocol on the same scanner. All except two PD subjects had both a DTBZ and a RAC scan; one PD subject lacked a DTBZ scan, and another lacked a RAC scan. Among the HC group, 6 subjects were scanned with both DTBZ and RAC, while another 4 subjects only had a DTBZ scan. Although the number of HC subjects was relatively low, their images were only used as part of a larger group that also included early PD subjects, as elaborated below. All HC and PD subjects underwent a T1-weighted structural magnetic resonance (MR) imaging scan.
Clinical and population characteristics of the HC and PD groups are listed in Table 1. To analyze the radiotracer binding patterns in different disease stages, the PD subjects were split into two sub-groups of approximately equal size according to the DD. The threshold that produced similarly-sized sub-groups was found to be 3 years. The sub-group with 0 � DD � 3 years (N = 22) is referred to as early PD, and the sub-group with DD � 4 years (N = 19) is referred to as moderate PD. Since chronic intake of dopamine receptor agonist treatment, but not levodopa, has been shown to affect RAC binding [31], we recorded the agonist levodopaequivalent dose (AgLED) for each PD subject.
The severity of movement impairment in the PD subjects was assessed off-medication (withdrawal period 12 hours) according to the motor subscale of the Movement Disorder Society's Unified Parkinson's Disease Rating Scale (MDS-UPDRS Part III) [32]. The better (less affected) and worse (more affected) sides were identified for each subject using the total MDS-UPDRS III scores for the left and right sides. As one of the clinical metrics for imagebased prediction, lateral motor score (LMS) was computed by adding the individual items of MDS-UPDRS III for the leg (leg rigidity + toe tapping + leg agility) and arm (arm rigidity + finger tapping + hand movements + pronation/supination) assessments. The LMS was computed separately for the better and worse sides ( Table 2). The tremor MDS-UPDRS III items were not included, as tremor is known to have worse correlation with dopaminergic deficit [33,34]. The study was approved by the University of British Columbia Ethics Board and all subjects gave informed written consent.

PET and MR image acquisition
PET scans of the PD subjects were performed with anti-Parkinson medication withdrawn at least 12 hours prior to scanning. Before tracer injection, transmission scans for attenuation correction were performed over 10 minutes with a rotating 137 Cs source. Subsequently, the subjects were administered 320 ± 34 MBq of DTBZ and 322 ± 34 MBq of RAC via intravenous injection. PET data were acquired in list mode using a High Resolution Research Tomograph (HRRT, Siemens, spatial resolution (2.5 mm) 3 ). Acquired coincidence data were binned into 16 temporal frames with frame durations of 1 min (×4 frames), 2 min (×3 frames), 5 min (×8 frames), 10 min (1 frame) (60 minutes in total). Dynamic images were reconstructed using the 3D list-mode ordinary Poisson ordered-subset expectation maximization (OP-OSEM) algorithm [35] with 16 subsets and 6 iterations, and corrections for attenuation, scattered and random coincidences. After reconstruction, the images were smoothed using a 2.0 mm (full width at half-maximum) Gaussian filter to reduce noise, and rigidly frame-to-frame realigned to correct for possible motion during the scans. The image dimensions were 256×256×207 with voxel size (1.219 mm) 3 . MR images were acquired using a T1-weighted Turbo Field Echo sequence (TR 7.7 milliseconds) on a Philips Achieva 3T scanner. The acquired MR image dimensions were 256×256× 170 with voxel size (1.0 mm) 3 . Post-acquisition, the MR images were resampled using trilinear interpolation to match the PET voxel size (yielding new dimensions 211×211×140 voxels). All subcortical structures included in the analysis were automatically segmented from the MR images using Freesurfer 6.0 [36].
Segmentations of the left and right putamen and caudate were used for the main part of the analysis. The ventral striatum (VS) regions, also segmented by Freesurfer (labeled as the accumbens area), are more closely related to the executive rather than motor function. Therefore, they were combined with the caudate segmentations to reduce the number of analyzed regions.

Computation of parametric PET images
DTBZ and RAC images averaged over 30-60 min post-injection were rigidly co-registered to the corresponding subject's MR images, using the SPM 12 software (www.fil.ion.ucl.ac.uk/ spm/); the quality of the registration was visually inspected for each subject. Following coregistration, parametric BR images of DTBZ and RAC were computed for each subject, by dividing the voxel values in the respective activity images by the mean activity in a reference region. Freesurfer segmentation of the cerebellum (gray and white matter) was used to define the reference region for RAC. Segmentation of the occipital cortex, produced by masking the Freesurfer segmentation of the cortex, was used to define the reference region for DTBZ.
We used parametric images of BR rather than more commonly used BP ND , since the BR can be readily computed from static as well as dynamic scans. In a preliminary study, the proposed method was applied to parametric BR as well as BP ND images, obtained using the simplified reference tissue model (SRTM2) [37,38]; the results were found to be very similar. This was expected, since the BR (measured after the tracer equilibration) has been previously shown to be highly correlated with BP ND values [39]. In addition, BR is often used as the outcome in imaging studies with AV-133, the 18 F-labelled version of 11 C-DTBZ [40]. In principle, PC-LASSO can be applied to different types of parametric images, since it primarily works with covariance patterns rather than absolute voxel values.

Adjustment of disease duration
We found that for DTBZ, the mean BR in the putamen was better fitted against DD by an exponential function, rather than by a linear function, in agreement with previously reported results reported with BP ND [41]. Thus, to enable the use of linear fits produced by LASSO, we linearized DD to obtain an adjusted DD (aDD) computed as aDD = exp(−0.17 × DD), followed by re-normalization between min(DD) = 0 years and max(DD) = 16 years. The coefficient −0.17 years -1 was found from the exponential fit of the mean DTBZ BR against DD. The re-normalization was applied to make aDD increase with disease progression, within the same range as DD. This facilitated the interpretation of the results (scatter plots and prediction errors).

Processing and PCA of PET images
The pipeline for processing of DTBZ and RAC parametric BR images is schematically illustrated in Fig 1. For each subject, the segmented (and co-registered) MR and BR images were separated into left and right sides, and each side was processed separately. Examples of MR images, striatum segmentations, and PET images for one of the PD subjects are shown in Fig 2A. On each side, the MR image-defined segmentations of the putamen and caudate were combined to generate a single labeled volume that was warped to a common striatal template ( Fig 2B) using 3D diffeomorphic mapping [42]. The resulting transformation was saved and applied to the respective BR image. Thus, the left and right striatal BR images of all HC and PD subjects were warped to a common labeled reference space (Fig 2C). The PD images were further sorted into better and worse brain sides, contralateral to the less and more clinicallyaffected body sides.
PCA was applied to the voxel BR data in the common space, separately for voxels in the putamen and caudate, and separately for the better and the worse brain sides. Each voxel in the common space was treated as a variable, and voxel BRs corresponding to different subjects were treated as observations. PCA finds orthogonal axes in the original variable space, termed PC loadings, that are aligned with the directions of greatest variance in the data. The following decomposition of the BR data into K PCs was performed: where ! f n is the vector of voxel BR values of the n-th subject, ! o k is the vector of voxel weights in the k-th PC, and s nk are the PC scores, quantifying the projections of the original data onto the new basis ! o k . In other words, PC scores s nk quantify the relative contribution of the k-th PC to the BR values of the n-th subject. Vectors ! o k are the orthogonal PC loadings that define the contribution of each voxel to the PCA-defined basis functions.

Spatial analysis of PC loadings
PC loadings obtained from PCA of all PD subjects were analyzed. The spatial patterns in voxel weights, captured by the PC loadings, were visually examined in the putamen and caudate. The better-side patterns were compared to the worse-side patterns. PC scores for each subject corresponding to these loadings/patterns were used in the multivariate PC-LASSO fits to obtain PC-LASSO models of aDD and LMS. Parametric PET images of the left and right striatum were warped to a common space using an MR-derived striatal template. In the common space, the images were re-sorted according to the better (less clinically affected) and worse sides, and used in PCA to obtain side-specific PC loadings and PC scores. The PC loadings and scores in the putamen and caudate were analyzed separately. LASSO was used to obtain the fitting coefficients between the PC scores and clinical PD metrics-aDD and LMS. Using the fitting coefficients, the PC loadings were linearly combined to obtain the PC-LASSO estimators. https://doi.org/10.1371/journal.pone.0206607.g001 Voxel-based analysis of brain PET images: Visualization and quantification of neurodegeneration patterns PLOS ONE | https://doi.org/10.1371/journal.pone.0206607 November 5, 2018 To assess the consistency of the patterns in different disease stages, PCA was additionally applied to the DTBZ and RAC caudate and putamen images of three subject sub-groups: The mixed HC/PD sub-group was constructed with the intention to capture the patterns of tracer binding associated with the onset of clinical disease. It consisted of the HC subjects, and Voxel-based analysis of brain PET images: Visualization and quantification of neurodegeneration patterns a subset of subjects from the early PD sub-group (with DD up to two years). Only one brain side (left or right) was used from each HC subject, chosen randomly. The age difference between the HC and PD subjects in the sub-group was not statistically significant (p > 0.1). One PD subject in the mixed sub-group was on agonist therapy-that subject was excluded from the PCA of RAC images. The early and moderate PD sub-groups (Table 1) were intended to capture the patterns of tracer binding at later stages of the disease.

PC-LASSO models for clinical metric prediction
Following the analysis of PC loadings, we used LASSO-based fitting with cross-validation to combine multiple PCs in a model optimized for prediction of clinical disease metrics. For a general set of input (predictor) and output (predicted) variables, LASSO finds a vector of fit- where N is the number of samples, x nj are the input (independent) variables, J is the number of input variables, y n is the outcome (dependent) variable, b 0 is the intercept, and λ is a non-negative regularization parameter. Compared to the regular least squares fitting, regularization that partially restricts overfitting is achieved through the inclusion of the term l P J j¼1 jb j j. The choice of λ controls the degree of regularization. Large λ forces the terms with small β j to go to zero. With λ ! 0, LASSO becomes identical to regular multivariate least-squares fitting, and with λ ! 1 all but the constant (β 0 ) terms are eliminated. LASSO is particularly well-suited for problems with a small sample size and a large number of independent variables. The models that were fitted to the data using LASSO had the general form where the PC scores s nj obtained from PCA of all available PD subjects are used instead of the variables x nj ; n represents the subject index, and j represents the PC number. The PCs were numbered (ranked) in a descending order according to their variance-accounted-for (VAF). Based on the initial analysis of PC loadings, it was determined that large clusters of voxels with similar weights were only present in the top 5 PCs; loadings of PC6 and onwards resembled random noise. Thus, to make the number of pre-optimization independent variables constant in different tested models, only the scores of the top 5 PCs were used as LASSO inputs (J = 5). The variable y n represents the predicted clinical metric; we fitted separate models to predict aDD, better-side LMS and worse-side LMS ( Table 2). With aDD, we tested better-side and worse-side PC scores as the input variables. In the models of better and worse LMS, PC scores from the corresponding contralateral side of the brain were used as inputs. All models with RAC data included age and AgLED as covariates. The total number (train+test) of samples in the cross-validated fitting was equal to 40 for both tracers.
Model optimization and fitting was performed using search over λ ranging from 0 to λ max in 100 steps, where λ max produced an intercept-only model (β j � 0). For each tested λ, the data were randomly divided 500 times into training and test sets using 0.3 holdout ratio: 70% of the data were used for training, and 30% for testing. The mean cross-validated mean squared error (MSE) in the test sets was measured (MSE test ). For λ that produced the lowest mean MSE test (λ min ), MSE on the entire data was measured (MSE all ). MSE test was used to estimate the degree of overfitting and predictive accuracy for the unseen (future) data. MSE all was used as a metric of model fit to the available data. Robustness of the LASSO fits corresponding to λ min was assessed from the means and standard deviations of the term coefficients β j .
The MSEs of the PC-LASSO models of aDD and LMS given by Eq 3 were compared to those of 1) a linear model that only included the mean BR values, and 2) a constant model that included only the intercept β 0 . The constant model was included in the analysis to obtain reference MSE values that indicate a lack of meaningful outcome prediction. To assess the robustness of the PC-LASSO models, we analyze the LASSO trace plots of the fitting coefficients and MSE test (i.e. their values as functions of λ). We compare the predicted-versus-actual (PVA) plots between the mean BR and PC-LASSO models, and provide the values and significance levels of the fitting coefficients β j for aDD and LMS.

Computation of PC-LASSO estimators
Since both PCA and LASSO are linear in nature, the LASSO-fitted combinations of PC scores that predict clinical metrics can be transformed back to the space of the original variables, i.e. voxels. Indeed, using the orthogonality of PCs and Eq 1, a LASSO-fitted linear model given by Eq 3 can be expressed as where the vectorṽ ¼

Analysis of PC loadings computed using all PD subjects
DTBZ PC loadings in the putamen and caudate are visualized in Fig 3A in the order of decreasing VAF. Loadings of PCs with relatively high VAF contained large clusters of voxels with similar weights, while loadings with low VAF had more heterogeneous weight distributions. PC1 had much greater VAF compared to other PCs, and the respective loadings contained only positive voxel weights; thus, mathematically it represents a weighted mean, and describes a global variance in the tracer binding. PC2 loadings contained both positive and negative weights, and captured antero-posterior gradient in the putamen and infero-superior gradient in the caudate. Patterns reflected by PC3 loadings combined antero-posterior and infero-superior gradients. Patterns in the loadings of PC4 and PC5 were more intricate; PC5 in the caudate had prominent voxel weights in the ventral striatum region. In the putamen, the difference between the better and worse side PC loadings (Fig 3A) reflected the disease asymmetry: the worse-side patterns can be seen as more progressed better-side patterns. For example, in the PC2 loadings, on the worse side the positive/negative weight boundary was shifted towards the anterior putamen (inferior caudate) compared to the better side. There were no similar side-to-side differences in the caudate PC loadings, possibly indicating that the asymmetry was not present or was less pronounced in the caudate.
Patterns captured by the RAC PC loadings were different from the DTBZ patterns, with the possible exception of PC1 (Fig 3B). For example, patterns captured by DTBZ PC2 and 3 in the putamen could not be found among the RAC PC loadings; the gradient direction in the RAC PC2 loadings was different compared to DTBZ. Also in contrast to DTBZ, the disease asymmetry was not clearly manifested in the RAC PC loadings on the better and worse sides.

PC loadings computed in different disease stages
Better-side DTBZ PC loadings computed in different subject sub-groups are visualized in  Voxel-based analysis of brain PET images: Visualization and quantification of neurodegeneration patterns the sub-groups consistently decreased with disease progression, while VAF by PC2-5 increased. The voxel weight distributions in the loadings of PC1, 2 and 3 were found to progressively change as a function of disease stage. For example, in the mixed HC/PD sub-group, PC1 loading had greater weights in the posterior putamen (superior caudate), while in early and moderate PD, the region with higher weights shifted towards the anterior putamen (inferior caudate); in the putamen PC3 loadings, the ventral region became increasingly more involved, while the magnitude of the antero-posterior gradient decreased.
With RAC, voxel weight patterns in the PC loadings were also similar between different sub-groups (S1 Fig). However, in contrast to DTBZ we did not observe a gradual change in RAC PC loadings with disease progression.

Performance and analysis of PC-LASSO models
The measured values of MSE test and MSE all for the DTBZ putamen models (constant, mean BR and PC-LASSO models) are summarized in Table 3. When putamen DTBZ data were used as the model inputs, PC-LASSO models of aDD and better-side LMS outperformed the mean BR and constant models. The largest improvement provided by PC-LASSO, compared to the mean BR, was in MSE all on the worse side of putamen (mean BR MSE all = 18.3, PC-LASSO MSE all = 13.4). Although it is known that MSE all has a tendency to diminish with greater number of independent variables, the fact that MSE test of PC-LASSO was also lower implies that the optimized model was not overfitted. In the caudate, the MSE test of all models were similar (Table 4).
When RAC data were used as the model inputs, the prediction errors were similar between the PC-LASSO, mean BR and constant models, indicating that very little aDD and LMSrelated information, if any, is captured by this tracer in general.
We examined in detail three fitted DTBZ PC-LASSO models with lower MSE test compared to that of the mean BR and constant models: 1. aDD predicted from the better-side PC 1, 3, and 4 scores in the putamen (scores that were reatined as inputs by LASSO) 2. aDD predicted from the worse-side PC 1-5 scores in the putamen 3. better-side LMS predicted from the contralateral PC 1, 2, and 4 scores in the putamen These three models had the general form given by Eq 3 with different input variables, coeffients and predicted variables. Trace plots of the fit coefficients β j and MSE test against λ for the three considered DTBZ PC-LASSO models are shown in Fig 5. The MSE test of the mean BR and constant models are plotted for reference. Distinct global minima were observed in the MSE test of all three PC-LASSO models; these minima were located at non-zero values of λ, and were considerably lower than the MSE test measured with the mean BR, PC1 scores alone, or constant models. With large values of λ, only the PC1 scores were retained.
The PVA data plots for the PC-LASSO models and the respective mean DTBZ BR models are shown in Fig 6. The plots demonstrate that the PC-LASSO models indeed fit the data better than the mean BR. Improvement of the aDD fit with PC-LASSO was observed over the entire range of actual aDD values. Improvement of the LMS fit was primarily observed in the intermediate range (2-6) of LMS; the high actual LMS values (LMS > 6) were severely under-estimated. The first two columns specify the predicted clinical metric (dependent variables) and the side of the brain from which the imaging data (independent variables) were taken. The last column specifies which terms were retained after optimization. The mean values and standard errors of MSE test from 500 test sets are shown. A square root of the MSEs provides the absolute values of prediction errors. b.s. = better side; w.s. = worse side; x-lat. = contralateral.
https://doi.org/10.1371/journal.pone.0206607.t003 The cross-validated mean values, standard deviations, and significance levels (p-values) of the fit coefficients β j , obtained with λ = λ min , are plotted in Fig 7. The mean values of the coefficients obtained from the training sets agree well with the values that were determined using the entire data. As indicated by the relative magnitudes of the means and standard deviations in Fig 7, the fit coefficients were consistent between different training sets. Interestingly, the contribution of PC4 in the model 3 was more significant than that of PC1. This may imply that the severity of motor impairment in PD is associated with neurodegeneration in specific subregions of the putamen, rather than in the entire structure.

DTBZ PC-LASSO estimators of aDD and LMS
The PC-LASSO estimators for aDD and LMS, corresponding to the three considered DTBZ models fitted on all data, are visualized in Fig 8. The better putamen estimator for aDD contained positive weights primarily in the anterior region, while negative weights were located in the interior and posterior regions. The worse putamen estimator for aDD featured negative weights in the anterior putamen, and positive weights in the posterior putamen-opposite to the corresponding better-side estimator for aDD. The estimator for the better-side LMS featured a positive-to-negative gradient in the ventral-to-dorsal direction.

Robustness and clinical relevance of patterns captured by the PC loadings
One of the main findings obtained with PCA was that tracer binding patterns captured by the PC loadings were generally similar between different sub-groups of subjects. This provides evidence that the patterns were robust and inherent to the disease, and are thus very likely to be Voxel-based analysis of brain PET images: Visualization and quantification of neurodegeneration patterns found in other PD cohorts. Comparison of DTBZ PC loadings between the sub-groups on a finer level uncovered the spatial evolution of patterns as a function of disease stage.
The pattern captured by PC1 describes a global, but not uniform, decrease in tracer binding. From the analysis of VAF it follows that in early disease, the majority of the variance between the subjects is explained by this pattern. The pattern of PC2 captured the antero-posterior gradient in DTBZ binding that develops in PD. The gradients captured by PC3, PC4 and PC5 describe variances of tracer binding in various directions that may be related to individual differences in disease progression. Although the gradient patterns of PC2 and higher order PCs reflected a much smaller component of variance compared to PC1, according to the VAF values they become more prominent with disease progression.
The difference between the better-and worse-side DTBZ PC loadings in the putamen was a reflection of the well-known asymmetrical nature of PD. Surprisingly, no corresponding side difference was found in the caudate loadings (Fig 3A). This suggests that the spatiotemporal progression of neurodegeneration may be symmetrical in the better and worse sides of the caudate, at least during the symptomatic stages of the disease.
Although it is very well known that pre-and post-synaptic processes respond differently to the disease, PCA was able to identify more accurately the spatial differences between the DTBZ and RAC binding alterations. Voxel-wise comparison of DTBZ and RAC second PCs in the putamen shows that pre-and post-synaptic functional gradients differ in terms of spatial localization. The predominant DTBZ gradient was in the antero-posterior direction, while the predominant RAC gradient was in the infero-superior direction (Fig 3).
In contrast to DTBZ, the binding patterns of RAC did not change with disease progression. These observations indicate that compensatory upregulation of the D2 receptor density may not exactly match the spatial localization of the dopaminergic terminal loss. The superior putamen is known to be primarily involved in the leg motor control [43]; it is thus possible to speculate that the limited spatial extent of the D2 receptor density upregulation, localized mostly in the superior putamen as captured by the putamen RAC PC2, is partially responsible for the fact that most new PD diagnoses are made based on the impairment of arm (but not leg) movement.

Performance and biological interpretation of the PC-LASSO estimators
In terms of the clinical metric prediction error, PC-LASSO models provided an improvement over the mean BR models when using DTBZ imaging data, but not when using RAC. The fact that PC-LASSO did not always outperform the mean BR or constant models demonstrates that a greater number of independent variables did not automatically produce better prediction accuracy, and validates the method as being resistant to overfitting.
The DTBZ PC1 score appeared to be the strongest single predictor of aDD and LMS as it was the only term retained at higher values of λ. However, the final optimized models fitted using λ min (with lowest MSE test and MSE all ) retained non-zero terms for PC2-5, i.e. terms that were different in meaning (and orthogonal) to PC1. This implies that the scores of PC2-5 had independent explanatory value while capturing clinically-relevant tracer distribution patterns.
Similarly to individual PC loadings, the obtained DTBZ PC-LASSO estimators for aDD and LMS (Fig 8) can be spatially interpreted. Taking into consideration the signs of the fit coefficients (Fig 7), it follows that on the better side of putamen, greater aDD is associated with increasing antero-posterior gradient, whereas on the worse side, it is associated with Voxel-based analysis of brain PET images: Visualization and quantification of neurodegeneration patterns decreasing gradient. The latter is likely a reflection of the tracer binding floor effect in moderate-to-advanced disease.
Another observation of interest is that according to the putamen PC-LASSO estimator for the LMS, motor impairment in limbs is associated with the magnitude of the ventral-to-dorsal gradient in DTBZ binding. The found association between the ventral-to-dorsal gradient and limb motor scores appears to be consistent with the known striatal somatotopy in healthy state, wherein the forelimb and hindlimb cortical regions project to the ventral-to-dorsal putamen [43]. The lack of accurate PC-LASSO predictions for higher LMS values likely implies that DTBZ (and RAC) imaging data only provide a partial characterization of the neural circuitry that is responsible for motor function. Notwithstanding this limitation, the PC-LASSO estimators for aDD and LMS visually demonstrate how different properties of the same PET image (in this case gradients along different directions) may be related to different clinical metrics.

Method limitations and future work
There are several limitations to the proposed method. First, LASSO was used to fit linear models, while the DTBZ binding is known to change non-linearly as a function of disease progression [44]. To overcome this limitation, we computed aDD for model fitting; however, the appropriate model for linearizing other disease metrics, or binding of other tracers, may not always be known a priori. Second, PCA does not uniquely determine the sign of ! o k in Eq 1. This property may complicate identification of similar PC loadings obtained in different datasets. In future work, it would be beneficial to establish a metric of similarity between PC loadings that could account for possible sign differences, as well as small differences in voxel weights (e.g. between better and worse sides). We envision the use of some form of clustering to group similar PC loadings together. Third, the study was performed on cross-sectional data, which are not optimal when investigating disease progression. Other limitations of the method are the facts that it requires non-rigid image registration, and may only be applicable to tracers that have consistent/reproducible binding patterns.
Various modifications of the PC-LASSO pipeline can be envisioned. Our rationale for using PCA here was that, among similar techniques, it is most general, and is well-suited to capture intensity gradients through the inclusion of both positive and negative voxel weights. PCA loadings are orthogonal by construction, and the corresponding scores are expected to be uncorrelated, which makes LASSO appropriate to perform regularized fitting. However, depending on the application, PCA in principle can be replaced with another linear dimensionality reduction technique, such as sparse and/or generalized PCA [45,46] or nonnegative matrix factorization [47]. While the linearity of PCA did not pose an obstacle in our work, a non-linear dimensionality reduction method (e.g. kernel PCA) may also be used where it is required. In turn, LASSO can be replaced with logistic LASSO, group LASSO, elastic net, ridge regression, or a more sophisticated non-linear machine learning algorithm.
It would be of interest to compare the predictive accuracy of PC-LASSO models to that of frequently-used local radiomic features in different imaging scenarios; given that we only had a limited number of subjects available, such comparison was beyond the scope of this work. Future work could also include the use of PC-LASSO or similar method to predict disease progression in large longitudinal imaging study with multiple tracers. A larger number of subjects will also allow for a direct comparison of predictive performance between the PC-LASSO models and most frequently used radiomic features.

Conclusion
We propose a novel data-driven method to construct models that predict clinical disease metrics from imaging data. The method is comprised of voxel-wise PCA followed by LASSO, and readily allows incorporation of clinical covariates at the same level as the voxel data. We applied PC-LASSO to the analysis of dopaminergic PET tracer binding in the striatum of PD subjects. The PC loadings obtained in different groups of subjects revealed predominant voxel-level binding patterns associated with the initial symptom onset and disease progression. The constructed DTBZ PC-LASSO models had lower cross-validated error of clinical metric prediction than the mean BR, confirming that patterns captured by the PC loadings are disease-related and offer additional explanatory and predictive power. The PC-LASSO estimators captured information in a non-local manner, and hence enabled data-driven visualization and interpretation of spatial patterns manifested in the images. The method is applicable to a variety of PET and SPECT imaging studies that focus on basal ganglia as well as other brain structures, including combination of data acquired with different tracers.