Joint pattern analysis applied to PET DAT and VMAT2 imaging reveals new insights into Parkinson's disease induced presynaptic alterations.

Most neurodegenerative diseases are known to affect several aspects of brain function, including neurotransmitter systems, metabolic and functional connectivity. Diseases are generally characterized by common clinical characteristics across subjects, but there are also significant inter-subject variations. It is thus reasonable to expect that in terms of brain function, such clinical behaviors will be related to a general overall multi-system pattern of disease-induced alterations and additional brain system-specific abnormalities; these additional abnormalities would be indicative of a possible unique system response to disease or subject-specific propensity to a specific clinical progression. Based on the above considerations we introduce and validate the use of a joint pattern analysis approach, canonical correlation analysis and orthogonal signal correction, to analyze multi-tracer PET data to identify common (reflecting functional similarities) and unique (reflecting functional differences) information provided by each tracer/target. We apply the method to [11C]-DTBZ (VMAT2 marker) and [11C]-MP (DAT marker) data from 15 early Parkinson's disease (PD) subjects; the behavior of these two tracers/targets is well characterized providing robust reference information for the method's outcome. Highly significant common subject profiles were identified that decomposed the characteristic dopaminergic changes into three distinct orthogonal spatial patterns: 1) disease-induced asymmetry between the less and more affected dorsal striatum; 2) disease-induced gradient with caudate and ventral striatum being relatively spared compared to putamen; 3) progressive loss in the less affected striatum, which correlated significantly with disease duration (p < 0.01 for DTBZ, p < 0.05 for MP). These common spatial patterns reproduce all known aspects of these two targets/tracers. In addition, orthogonality of the patterns may indicate different mechanisms underlying disease initiation or progression. Information unique to each tracer revealed a residual striatal asymmetry when targeting VMAT2, consistent with the notion that VMAT2 density is highly related to terminal degeneration; and a residual DAT disease-induced gradient in the striatum with relative DAT preservation in the substantia nigra. This finding may be indicative either of a possible DAT specific early disease compensation and/or related to disease origin. These results demonstrate the applicability and relevance of the joint pattern analysis approach to datasets obtained with two PET tracers; this data driven method, while recapitulating known aspects of the PD-induced tracer/target behaviour, was found to be statistically more robust and provided additional information on (i) correlated behaviors of the two systems, identified as orthogonal patterns, possibly reflecting different disease-induced alterations and (ii) system specific effects of disease. It is thus expected that this approach will be very well suited to the analysis of multi-tracer and/or multi-modality data and to relating the outcomes to different aspects of disease. Crown

Introduction

Parkinson's disease (PD) is the second most frequent progressive neurodegenerative disorder (de Lau and Breteler, 2006). The motor deficit of PD is traditionally associated with dysfunction of the nigrostriatal pathway, characterized by progressive loss of dopaminergic neurons in the substantia nigra and loss of their projection fibres to the striatum (Stoessl, 2012). Imaging studies show that the motor deficits start to become clinically relevant when 30–50% of nigral dopaminergic cells are lost (Brooks, 2010). Neurodegeneration of the dopaminergic system tends to follow a fairly well defined spatio-temporal pattern in which the dorsal posterior putamen contralateral to the more affected body side is affected first, followed by degeneration in the ventral and anterior putamen and the caudate, as shown in Fig. 1(Stoessl, 2012). The relatively long preclinical stage of PD, in which subjects remain asymptomatic despite significant dopaminergic neuronal loss, may be due to potential compensatory effects taking place at different stages of dopamine (DA) processing — including DA synthesis, release and turnover (Lee et al., 2000; Nandhagopal et al., 2011; Sossi et al., 2004). Such compensatory mechanisms, likely also involving other neurotransmitter systems (Bezard et al., 2003; Liu et al., 2018; Wile et al., 2016), are deemed responsible for minimizing the effects of dopaminergic deficits on the clinical behaviour prior to onset of motor deficits and are thought to persist in the very early stages of disease (Lee et al., 2000).

Fig. 1

[11C]-dihydrotetrabenazine (DTBZ) PET image (left) and [11C]d-threo-methylphenidate (MP) PET image (right) for a Parkinson's disease (PD) subject. PD subject showed characteristic asymmetric tracer uptake in the less and more affected hemispheres. PD subject also showed spatio-temporal pattern of dopaminergic loss with the posterior putamen (putamen 3) affected before the anterior putamen (putamen 1) and caudate. PET = Positron Emission Tomography.

There is now established recognition that PD is not just a motor disorder; patients often experience non-motor deficits alongside or even before the onset of motor deficits. Non-motor deficits may be more closely related to DA projections outside the nigrostriatal pathway, in addition to alterations in other neurotransmitter systems such as the cholinergic and serotonergic systems (Fu et al., 2018; Liu et al., 2018; Wile et al., 2016). Outside the nigrostriatal pathway, the mesocorticolimbic pathway transmits DA from the ventral tegmental area (VTA) to the ventral striatum (VS) and to the prefrontal cortex (Weingarten et al., 2015); the tuberoinfundibular pathway transmits DA from the hypothalamus to the pituitary gland (Dailly et al., 2004). Study reported that deficient hypothalamic DA transmission may play a role in autonomic and endocrine abnormalities in PD (Javoy-Agid et al., 1984). Studying DA processing in these pathways in addition to the nigrostriatal pathway may provide a more complete picture of dopaminergic denervation in PD, especially in the early disease stage or even before disease onset. Changes in the metabolic and functional connectivity have also been identified at various stages of PD and related to specific clinical manifestations such as cognitive deficit. It is thus reasonable to postulate that PD imprints a general disease-related pattern on several aspects of brain function, with alterations in selected systems reflecting either a specific system disease response or specific clinical manifestation of the disease. The ability to identify such patterns and their evolution as a function of disease progression/specific clinical manifestation may aid in the understanding of disease mechanisms or subjects' propensity towards a particular clinical trajectory.

In this work, we propose the use of a novel joint pattern analysis to study functional similarities and differences between multiple PET targets. Our analysis enhances the more traditional approach where the relationship between two or more sets of imaging data is examined using univariate approaches such as correlation and t-test. Multivariate techniques, such as principal component analysis (PCA) and independent component analysis (ICA) (Eidelberg, 2009; Vo et al., 2017), have been used to decompose individual datasets into functional networks. However, instead of analyzing individual datasets (i.e. comparing functional networks obtained separately from each dataset), we used a data fusion approach to explore and identify common and unique information given by each dataset as functional networks. These common and unique functional networks can provide additional and more direct insights into the interactions between processes observable with different tracers and the differential information provided by each individual tracer. This approach is particularly suitable for neurodegenerative diseases such as PD, where disease affects different stages of DA processing as well as multiple neurotransmitter systems.

Many different joint multimodal analysis techniques have been developed in the neuroimaging field, mainly applied to magnetic resonance imaging (MRI) data. One such data fusion approach commonly used is joint independent component analysis (jICA)(Calhoun and Sui, 2016). More recently, canonical correlation analysis (CCA) has gained popularity. Unlike jICA, CCA provides a relatively less constrained solution to the data fusion problem (Correa et al., 2009, Correa et al., 2010). While jICA assumes that different datasets have exactly the same inter-subject covariations, the CCA models the coherence in the inter-subject covariations to identify associations between datasets (Correa et al., 2009, Correa et al., 2010). CCA has been successfully applied to the analysis of functional MRI (fMRI), electroencephalography (EEG), electromyography (EMG), structural MRI and behavioural data in PD and schizophrenia to explore common inter-subject variations in different datasets (Correa et al., 2009). After extracting the common information among different datasets, the unique information still remaining in each individual dataset can be extracted by using orthogonal signal correction (OSC), which was first introduced as a spectral preprocessing method in spectroscopic calibrations (Fearn, 2000; Wold et al., 1987). OSC was later used together with CCA to draw unique information from EEG and EMG data (Chen et al., 2013).

As this is a first application of such methodology to PET data, we chose to perform the analysis on data obtained from two fairly well characterized presynaptic PET tracers: [11C]-dihydrotetrabenazine (DTBZ) and [11C]d-threo-methylphenidate (MP). The vesicular monoamine transporter type 2 (VMAT2) binding measured by [11C]-DTBZ is proportional to the DA terminal density (Masuo et al., 1990; Stoessl, 2012) and is used to estimate vesicular uptake and storage of DA. DA reuptake, mediated by the membrane DA transporter (DAT), can be targeted in PET by [11C]-MP. Both VMAT2 and DAT are mainly located at the nerve terminals, but can also be found in cell bodies and axons (Fazio et al., 2018; Sulzer et al., 2016). DAT was shown to contribute to maintaining relatively constant synaptic DA levels by removing extracellular DA (Block et al., 2015; Sossi et al., 2007), and may be related to compensatory mechanisms in populations at higher risk of PD (Adams et al., 2005; Nandhagopal et al., 2008) and in early PD (Lee et al., 2000; Nandhagopal et al., 2009). On the other hand, reduction in VMAT2 binding is deemed a more direct measure of dopaminergic degeneration and is less susceptible to disease-related compensation (Vander Borght et al., 1995). Despite the functional differences of VMAT2 and DAT, traditional analyses of PET data showed high correlation between striatal alterations in DTBZ and MP binding in PD (Lee et al., 2000; Nandhagopal et al., 2011) and some reported no differential regulation of the striatal uptake of VMAT2 and DAT (Karimi et al., 2013; Tian et al., 2012). In addition, VMAT2 and DAT distributions in the striatal regions have been studied extensively with traditional univariate analysis, but investigations of their distributions outside the nigrostriatal pathway have been quite limited, especially with imaging studies.

In this study, we first examine the applicability and robustness of the proposed joint pattern analysis approach. We then examine the novel information provided by this approach compared to traditional univariate analysis; specially, we

Compare the common information between VMAT2 and DAT distributions to the results from the traditional univariate analysis to test the method's abilities to capture the characteristic dopaminergic patterns in the striatum. This serves as the main validation of the method.

Explore the decomposition of the common information between VMAT2 and DAT in the orthogonal spatial patterns as reflecting different/independent underlying disease-related mechanisms.

Interpret the unique information specific to VMAT2 or DAT distributions in light of the specific target behavior in the early stages of PD.

Materials and methods

Study participants

The study included 15 early sporadic PD subjects (9 males and 6 females). Exclusion criteria included clinical history of depression, active anti-depressant therapy or medication and a Body Mass Index (BMI) >35. Disease duration was estimated as time from onset of motor symptoms as reported by the subjects. PD subjects were clinically evaluated using the Movement Disorder Society Unified Parkinson's Disease Rating Scale Part III (MDS-UPDRS Part III) and Hoehn and Yahr scale to assess motor dysfunction, Montreal Cognitive Assessment (MoCA) to assess cognitive performance and Beck Depression Inventory (BDI). Detailed clinical characteristics are listed in Table 1. All assessments were performed off medication. The study was approved by the Clinical Research Ethics Board of the University of British Columbia and all subjects provided informed written consent.

Table 1

Clinical characteristics of all subjects. All numbers are reported as mean ± standard deviation.

Number of subjects	Age (years)	Disease duration (symptoms, months)a	Disease duration (diagnosis, months)b	MDS-UPDRS part III	Hoehn and Yahr scale	MoCA	BDI	Levodopa equivalent dose (mg)
15	59 ± 8	56 ± 34	44 ± 29	17 ± 9	1.6 ± 0.5	28.0 ± 1.5	4.4 ± 3.5	380 ± 220

PD = Parkinson's disease subjects; MDS-UPDRS = Movement Disorder Society Unified Parkinson's Disease Rating Scale; MoCA = Montreal Cognitive Assessment; BDI=Beck Depression Inventory.

Disease duration estimated as the time from onset of motor symptoms as reported by the patients.

Disease duration estimated as the time of clinical diagnosis.

Scanning protocols

All study subjects underwent DTBZ and MP PET scans and a T1-weighted MRI scan of the brain. The PET scans were performed on a Siemens High Resolution Research Tomograph (HRRT, Knoxville, TN) with a spatial resolution of 2.5 mm3 (Jong et al., 2007). Subjects were positioned using external lasers aligning the gantry with the inferior orbital-external meatal line, and custom fitted thermoplastic masks were applied to minimize head movement. Prior to PET scans, subjects were withdrawn from all anti-parkinsonian medications for at least 12 h. An average of 300 MBq with average specific activity of 10,194 Ci/mmol of DTBZ and MP were administered by intravenous injection over 60 s using an infusion pump (Harvard Instruments). Tracer injections were separated by at least 2.5 h to allow radioactive decay. Acquired data were binned into 16 time frames (frame durations: 4 × 60 s, 3 × 120 s, 8 × 300 s, 1 × 600 s; image dimension = 256 × 256 × 207; voxel size = 1.22 mm3) with a total duration of 60 min. Transmission scans required for attenuation correction were performed over 10 min with a rotating 137Cs source. PET images were reconstructed using the 3D list-mode ordinary Poisson Ordered Subset Expectation Maximization (OP-OSEM) algorithm (Comtat et al., 2004) with 16 subsets and six iterations, with corrections for decay, dead-time, normalization, attenuation, scattered and random coincidences. After reconstruction, images were smoothed with a 3.0-mm full-width at half maximum (FWHM) Gaussian filter to reduce noise. The frames were spatially realigned with rigid-body transformation to minimize the impact of motion during scans. The structural MRI scans were performed on a Philips Achieva 3.0 T MRI scanner (Phillips Healthcare, Best, NL) using the T1 turbo field echo (TFE) sequence (TR/TE = 7.7/3.6 ms; TFE shots = 218; flip angle = 8°; image dimension: 256 × 256 × 170; voxel size 1 mm3).

Image processing and analysis

The anatomical MRI image of each subject was first coregistered with the subject's mean PET image using statistical parametric mapping (Wellcome Trust Centre for Neuroimaging, University College London). Striatal regions of interest were manually placed on an averaged PET image derived from nine consecutive image slices (slice thickness 1.22 mm) spanning the axial extent of the striatum, using MRI image as guidance. Five elliptical ROIs were placed on the striatum bilaterally – one on the caudate, one on the VS and three covering the full length of the putamen (anterior – putamen 1, middle – putamen 2 and posterior – putamen 3). The same set of image slices was also used to define the occipital cortex reference region for both tracers.

In addition to the manually defined striatal ROIs, we also developed a ROI template in Montreal Neurological Institute (MNI) space using MRI images of healthy controls. The PET-coregistered MRI images were transformed to the MNI space. The ROI template was inverse-transformed to match each subject's PET image for further analysis. This ROI template contained four ROIs placed bilaterally (substantia nigra, thalamus, globus pallidus and hypothalamus), and four individual ROIs (posterior midbrain, pons, raphe nucleus and VTA). Combined with the manually-placed striatal ROIs, this yielded a total set of 22 ROIs. Regional time-activity curves were extracted from each ROI and the Logan graphical method (Logan et al., 1990) was used to calculate the non-displaceable binding potential (BPND) values using time ranges from 17.5 to 60 min. BPND values were then rearranged into the more and less affected hemispheres for all bilateral ROIs based on the average DTBZ BPND values in the three putamen ROIs.

Univariate analysis

To examine the distributions of VMAT2 and DAT density, we first performed one-sample t-test (one-tailed) on the BPND values for DTBZ and MP separately for 22 ROIs. For each tracer, ROIs with BPND values significantly greater than zero (denoted as ‘significant binding’, p < 0.05 after correcting for multiple comparisons) were used for joint pattern analysis. Unpaired two-sample t-test was performed between the more and less affected striatal regions in DTBZ and MP separately to examine disease-induced asymmetry. We then tested the correlation between DTBZ and MP BPND values in all ROIs to check the correlation strength between the two datasets. We also tested the correlation between DTBZ and MP BPND values in the caudate and putamen with disease duration to check the correlation strength of univariate measures for tracking disease progression in the early stages of disease. False positive rates were controlled at p = 0.05 using Bonferroni-Holm's step-down procedure (Holm, 1979).

Joint pattern analysis

We first applied the joint pattern analysis approach to DTBZ and MP BPND values in all ROIs (input data) that had significant DTBZ and MP binding to extract the common (canonical variates from CCA) and unique (from OSC) subject scores and the associated spatial binding patterns in the two datasets. We then compared these patterns with results obtained with the ten striatal ROIs only to indirectly examine the contributions of the extrastriatal regions to the patterns and the robustness of the method. Correlation analysis was performed between the subject scores and clinical measures.

The joint pattern analysis approach applied to the DTBZ (X) and MP (Y) datasets is performed as follow:

Input data matrices X and Y have dimensions [NxM1] and [NxM2] respectively. M1 and M2 are the number of imaging features (in this case, they are either the ten striatal ROIs (M1 = M2 = 10) or the ROIs with significant binding determined by univariate analysis for either DTBZ (M1) or MP dataset (M2)). N is the number of subjects.

Each input data matrix (X and Y) is first demeaned (Xdemean and Ydemean) and whitened (Xwhiten and Ywhiten). The whitening transformation first decorrelates the features in each input data matrix, so that the new data dimensions are linearly independent (orthogonal); it then transforms the covariance matrix into an identity matrix, which ensures the variance of the data along each new dimension is equal to one.where E and Σ are the eigenvector and eigenvalues of Xdemean.

This step serves two important purposes: (i) to reduce feature dimension of a rank-deficient input matrix into fewer components, so that N ≥ max (rank(X), rank(Y)); (ii) to scale all variables to have the same variance so that each variable is assigned equal importance in the subsequent analysis. In this case, M1 and M2 are reduced into top five components to minimize the noise content dominant in later components, while still maintaining at least 90% of the original variance.

CCA (Correa et al., 2009) is then applied to Xwhiten and Ywhiten (both matrices now have dimension [Nx5]). CCA identifies linear relationships between the two datasets to determine the inter-subject covariance. It seeks two mixing matrices (W1 and W2) such that each pair of canonical variates Ui and Vi (i = 1…5) has maximum correlation across the two datasets, while the canonical variates within each dataset are orthogonal (Ui and Uj are uncorrelated)

For i = 1…5,

Thus, the transformed data (canonical variates U and V) contain common (maximally correlated between two datasets) subject profiles, which are composed of subject score of each subject (representing the subject weights for the corresponding mixing matrix). Subject scores are in Z-score form with a mean of zero and a standard deviation of one.

Least absolute shrinkage and selection operator (LASSO) (Klyuzhin et al., 2018a; Tibshirani, 2011) is then applied to regress the canonical variates (Ui and Vi) from the original datasets X and Y to compute regression coefficients (CCA weights A and B) and regression residuals (Xresidual and Yresidual) as shown in Fig. 2. Ten-fold cross validation was used to estimate the best lambda with the cross-validated minimum square error, where lambda is the LASSO penalty coefficient

Fig. 2

Illustration of the decomposition and regression step. X and Y are the whitened input matrices (feature by subject) of non-displaceable binding potential (BPND) values obtained from step 2. The transformed data (canonical variates) U and V are calculated using CCA in step 3, which contains the most highly correlated subject scores along each component (in this case 5). The CCA weights matrices (A and B) are the regression coefficients from least absolute shrinkage and selection operator (LASSO) in step 4. Xresidual and Yresidual are the regression residuals. CCA = canonical correlation analysis.

Since the residuals from step 4 (X and Y) may contain information specific to each dataset besides noise, OSC (Fearn, 2000) is then applied to extract the largest orthogonal component from the LASSO residuals deemed to represent tracer-specific unique information, including unique subject scores (Uunique), unique CCA weights (Aunique), and true noise (X) for each dataset

CCA loadings are defined as the correlation coefficients between each canonical variate (Ui or Vi) and each column of X or Y (feature values for all subjects). CCA loadings represent the feature/region contributions to each pair of canonical variates and are used to construct the spatial patterns.

To determine the significance levels of the correlation between each pair of extracted canonical variates (Ui and Vi), a non-parametric permutation test is performed on the original datasets X and Y with 1000 iterations to construct the empirical null distributions of the correlation coefficients for each pair of canonical variates. The p-value of the original correlation can then be computed as the probability of observing a value at least as extreme as the original correlation in the null distributions. The correlation between the pairs of canonical variates is considered statistically significant if the p-value is <0.05.

To test the stability of the CCA weights and loadings, leave-one-out validation test is performed to compute the error bounds of the feature contributions. CCA loadings are considered statistically significant if the correlation p-value is <0.05 after correcting for multiple comparison.

All codes were written in Matlab and are available upon direct request to the corresponding author, however PET data used in this study are not made available publicly due to patients confidentiality.

Results

DTBZ BPND values were significantly greater than zero (p < 0.05 corrected) in all 22 ROIs, while MP BPND values were not significantly greater than zero in hypothalamus, posterior midbrain, pons, VTA and raphe nucleus (p > 0.05 corrected). Therefore, all 22 ROIs were included for DTBZ and 16 ROIs were included for MP in the joint pattern analysis. Detailed results from univariate analysis are included in the Supplementary Materials.

One subject (S15) appeared as outlier (fell outside the 95% confidence interval) when correlating BPND values with disease duration (Fig. 3). This subject had a disease duration of 23 months, but had the highest BPND values in all striatal regions for both DTBZ and MP (BPND values were more than two standard deviations higher compared to average BPND values in all subjects in most striatal regions). Without this subject, correlations between disease duration and average DTBZ and MP BPND values in the less affected putamen were stronger (R2 = 0.70, p < 0.001 for DTBZ; R2 = 0.45, p < 0.01 for MP). In order to find the best dopaminergic patterns related to disease, we first excluded this subject in the joint pattern analysis, then included this subject in to examine the effect of this outlier on the results.

Fig. 3

Scatter plots for average DTBZ and MP BPND values in the less affected putamen versus disease duration (estimated from the time of symptoms onset) in months. Both DTBZ (left) and MP (right) BPND values correlated significantly with disease duration. S15 fell outside the 95% confidence interval. BPND = non-displaceable binding potential. DTBZ = dihydrotetrabenazine. MP = methylphenidate.

Common information was obtained using BPND values in 22 ROIs for DTBZ and the 16 ROIs for MP that exhibited significant tracer binding. For both DTBZ and MP, the top five whitened components together accounted for 91% of the variance in the original datasets. Each of the five whitened components accounted for at least 7% of the variance. The top three pairs of canonical variates were significantly correlated between DTBZ and MP after permutation test (p < 0.05). The fourth and fifth pairs of canonical variates did not show high correlation across datasets (R2 < 0.5) and were not significant after permutation tests (Table 2), therefore are not discussed in later sections.

Table 2

Correlation strength R2 and significance between each pair of canonical variates.

Pairs of canonical variates	1	2	3	4	5
Correlation R2	0.98	0.90	0.85	0.47	0.28
Permutation p-value	0.048⁎	0.033⁎	0.001⁎	0.123	0.054

Significant correlation at p = 0.05.

The DTBZ pattern along the first pair of canonical variates (R2 = 0.98 between subject scores in the two datasets) showed significant negative loadings in the more affected striatal regions (caudate and putamen), and significant positive loadings in the substantia nigra, hypothalamus and pons (p < 0.01); the MP pattern along this canonical variate showed significant negative loadings in the more affected striatal regions (caudate, anterior putamen (putamen 1) and middle putamen (putamen 2)), and significant positive loadings in the less affected substantia nigra and thalamus (p < 0.05) (Fig. 4A). Along the second pair of canonical variates (R2 = 0.90 between subject scores in the two datasets), the DTBZ pattern showed significant positive loadings in the less affected caudate, VS and VTA, and significant negative loadings in the thalamus and globus pallidus (p < 0.05); the MP pattern showed significant positive loadings in the caudate and VS (Fig. 4B). Subject scores along the first and second pairs of canonical variates did not correlate with any clinical measures.

Fig. 4

Common spatial patterns along the first three pairs of canonical variates for DTBZ and MP. Stars indicate the ROIs with significant CCA loadings. ROI = region of interest; CCA = canonical correlation analysis; GP = globus pallidus; VS = ventral striatum; SN = substantia nigra; VTA = ventral tegmental area. DTBZ = dihydrotetrabenazine. MP = methylphenidate.

The spatial patterns for both DTBZ and MP along the third pair of canonical variates included significant negative loadings in the less affected caudate and putamen (p < 0.01) (Fig. 4C). The subject scores along the third pair of canonical variates (R2 = 0.85 between subject scores in two datasets) correlated significantly with disease duration for both DTBZ (R2 = 0.70, p < 0.001) and MP (R2 = 0.51, p < 0.01) (Fig. 5). Correlations with disease duration remained significant without the subject with longest disease duration (132 months) for both DTBZ (p < 0.001) and MP (p < 0.05).

Fig. 5

Correlation between subject scores and disease duration as estimated from the time of symptoms onset (months) for DTBZ and MP along the third pair of canonical variates. DTBZ = dihydrotetrabenazine. MP = methylphenidate.

The unique DTBZ pattern highlighted the asymmetry between the less and more affected striatal regions with significant negative loadings in the more affected striatal regions (caudate, putamen and VS), globus pallidus and VTA, and significant positive loadings in the pons (Fig. 6 top). The unique MP pattern showed significant positive loadings in the less affected posterior putamen (putamen 3), substantia nigra and thalamus, and significant negative loadings in the more affected anterior (putamen 1) and middle putamen (putamen 2) (Fig. 6 bottom). The unique subject scores for DTBZ and MP patterns did not correlate with any clinical measures.

Fig. 6

Unique spatial patterns along the first three pairs of canonical variates for DTBZ (top) and MP (bottom). Stars indicate the ROIs with significant CCA loadings. ROI = region of interest; CCA = canonical correlation analysis; GP = globus pallidus; SN = substantia nigra; VS = ventral striatum; VTA = ventral tegmental area. DTBZ = dihydrotetrabenazine. MP = methylphenidate.

Including S15, the correlations between disease durations and the subject scores along the third canonical pairs were weaker but still significant for both DTBZ and MP (p < 0.01 for DTBZ and p < 0.05 for MP). Regions with significant contributions to the common and unique spatial patterns remained similar. The correlation strength between the common information in DTBZ and MP along the third canonical pair and disease duration was weaker when only ten striatal ROIs were included in the analysis (with ten striatal ROIs only, R2 = 0.63 and p < 0.001 for DTBZ and R2 = 0.31 and p < 0.05 for MP). The common and unique spatial patterns for both DTBZ and MP however remained the same for the striatal regions as they manifested when all ROIs were included in the analysis.

Discussion

With the univariate analysis, we found all examined ROIs in this study showed DTBZ BPND values significantly greater than zero, while 16 out of 22 ROIs showed MP BPND values significantly greater than zero. The joint pattern analysis decomposed the characteristic gradients of dopaminergic loss in the striatum into orthogonal components ranked by the degree of commonality shared between VMAT2 and DAT distributions: 1) disease-induced asymmetry between the less and more affected dorsal striatum; 2) disease-induced gradient with caudate and ventral striatum being relatively spared compared to putamen; 3) progressive loss in the less affected striatum, which correlated significantly with disease duration. The unique information revealed differences between VMAT2 and DAT distributions.

Our results from univariate analysis, performed primarily to serve as reference for the outcomes of the spatial pattern approach, agree with widely reported findings on the dopaminergic deficit distribution in the striatum (Fazio et al., 2018; Kaasinen and Vahlberg, 2017; Nandhagopal et al., 2009) and are therefore in keeping with existing knowledge about the disease. While not of primary relevance to this study, the results obtained from the extrastriatal regions represent some novel findings of interest.

DTBZ BPND values were significantly greater than zero in regions involved in the nigrostriatal and mesocorticolimbic pathways and brain stem regions in early PD. Previous in-vivo 18F-FP-(+)-DTBZ PET imaging study in healthy controls (Lin et al., 2013) showed highest VMAT2 level in striatal regions and substantia nigra, followed by regions involved in the mesolimbic pathway, brain stem regions and thalamus. In healthy brains, VMAT2 level in the substantia nigra, hypothalamus and raphe nucleus is approximately 40% of those in the anterior putamen, and VMAT2 level in the posterior putamen is approximately the same as in the anterior putamen (Lin et al., 2013). In our case of early PD, DTBZ binding in the substantia nigra and raphe nucleus were also approximately 40% of those in the less affected anterior putamen, but higher than that in the posterior putamen. DTBZ binding in the hypothalamus and thalamus are approximately 45% and 6% of the binding estimated in the anterior putamen in normal brains (Lin et al., 2013), while 60% and 28% was observed in early PD. The hypothalamus and thalamus seem to have better preserved dopaminergic integrity compared to the anterior putamen; however, since VMAT2 is expressed by all monoamine neurons, more preserved DTBZ binding may be also reflective of noradrenergic instead of dopaminergic innervation.

DAT distribution was less widely spread outside the striatal regions compared to VMAT2. We observed asymmetric MP BPND values significantly greater than zero in the striatum and substantia nigra in early PD, which agrees with previous imaging finding of an asymmetric reduction of DAT level in the same regions in early PD (Fazio et al., 2018). We also observed significant MP binding in the thalamus and globus pallidus and insignificant binding in the hypothalamus and brain stem regions; these imaging results are consistent with previous post-mortem immunohistochemical studies which showed significant DAT expression in the thalamus (healthy human and non-human primates) and globus pallidus (healthy human) (Fazio et al., 2018; Prensa et al., 2000; Sanchez-Gonzalez, 2005), and no detectable DAT level in the hypothalamus (Koblinger et al., 2014) and brain stem regions (Sharma et al., 2018).

Common information

DTBZ and MP showed highly correlated subject scores along the first three pairs of canonical variates, corresponding to three distinct orthogonal spatial patterns. Of these three, only the subjects scores along the third pair correlated with disease progression.

Spatial pattern along the first pair of canonical variates showed the familiar early disease-induced asymmetry between the less and more affected striatum in both VMAT2 and DAT distributions as shown in Fig. 4 A. In univariate analysis, differences between the less and more affected striatal regions in either DTBZ or MP binding were not significant after correction for multiple comparison. Pattern analysis accurately captured this characteristic asymmetric tracer reduction independently of the number of regions involved in the analysis, indicating superior robustness of this approach. Higher asymmetry in the dorsal striatum was also associated with more preserved binding in the substantia nigra, hypothalamus and pons for DTBZ and substantia nigra and thalamus for MP, consistent with the fact that asymmetry in the dorsal striatum appears highest at clinical disease onset and decreases over time (Nandhagopal et al., 2009). This particular pattern may be characteristic of disease presence in this range of disease duration rather than progression. In addition, the relatively higher preservation of dopaminergic function in the substantia nigra in early disease appears consistent with recent imaging finding showing greater DAT loss at the axonal terminals compared to cell bodies in early PD (Fazio et al., 2018) and may provide support for the hypothesis of an early involvement of synapses and pre-terminal axons in the neurodegenerative process followed by alterations of cell bodies (Burke, 2013; Kurowska et al., 2016).

Spatial patterns along the second pair of canonical variates reflected the well-known disease-induced rostro-caudal gradient in early disease (Fig. 1). This gradient did not correlate with disease duration, suggesting that it may be predominantly a characteristic evolving in the preclinical stage where there is still relative preservation of the terminal density/dopaminergic function in the caudate and VS. This binding pattern in striatum was further associated with relatively decreased DTBZ binding in the thalamus and globus pallidus compared to the MP pattern where these regions showed no significant contributions.

Spatial patterns along the third pair of canonical variates reflected the progressive loss of dopaminergic function in the less affected striatum for both DTBZ and MP. This is consistent with previous findings (Klyuzhin et al., 2018b; Nandhagopal et al., 2009) which indicate that even though dopaminergic tracer uptake in the more affected striatum may be more sensitive for disease discrimination, tracer uptake in the less affected striatum provides a better marker to track disease progression. The inclusion of extrastriatal regions into the pattern analysis increased the pattern's correlation strengths with disease durations for both tracers compared to using the striatal regions alone; this indicates that dopaminergic denervation in other regions is also affected by disease progression in spite of the fact that the regional loadings in extrastriatal regions were by themselves not significant. The correlation of the MP pattern with disease duration was found to be stronger than what observed with univariate analysis, and the correlation strength of the DTBZ pattern with disease duration was similar to the univariate analysis applied to the averaged less affected putamen without correction for multiple comparison. However, the pattern analysis results did not suffer from multiple comparison problem that may decrease the statistical robustness of the outcomes obtained with the univariate analysis.

In addition, the orthogonality (which loosely implies independency) of the three spatial patterns may imply that disease-induced asymmetry, disease-induced gradients and denervation progression might be underlined by different or independent mechanisms. In a previous study (Nandhagopal et al., 2009) the asymmetry between the less and more affected striatal sides was shown to decrease as disease progresses; the first and third common patterns both showed striatal asymmetry, but with higher loadings on either the more affected (first pattern) or less affected (third pattern) sides. The orthogonality between the two patterns may reflect the fact that different mechanisms may be of most relative relevance to the less and more affected striatal sides related to the fact that degeneration in each side occurs at different stages of disease; however, whether this difference is due to ‘saturated’ dopaminergic loss in the more affected striatum or different underlying mechanisms still need further investigation. The same previous study (Nandhagopal et al., 2009) also showed that there is a marked rostro-caudal gradient of dopaminergic deficit at disease onset. We observed that the striatal gradient is independent of the striatal asymmetry, suggesting the two aspects of dopaminergic denervation may be induced by different mechanisms underlying disease initiation or progression. Interestingly, this decomposition may thus also provide guidance to determine metrics that are either more sensitive to disease discrimination (first two pairs of canonical variates) or are better suited to track disease progression (third pair of canonical variates). Inclusion of data from healthy controls in the analysis may help to further explore this hypothesis.

Overall, the common information in DTBZ and MP binding confirms the applicability and robustness of the proposed joint pattern analysis approach, and was shown to be more sensitive to specific spatio-temporal changes compared to univariate analysis.

Unique information

While it can be assumed that common information mainly reflects characteristic disease-induced alterations related to the integrity of dopaminergic function, unique information reflects patterns in which VMAT2 and DAT are differently affected in early disease. Unique DTBZ pattern showed asymmetry between the less and more affected striatum, similar to the common DTBZ pattern along the first pair of canonical variates, suggesting that VMAT2 may be more sensitive to direct disease effects, i.e. dopaminergic terminal degeneration, compared to DAT. Indeed VMAT2 density is deemed to be least sensitive to disease-induced regulatory changes (Vander Borght et al., 1995). The globus pallidus and VTA appeared relatively more affected and the pons were more preserved.

DAT unique pattern showed more reduced tracer binding in the more affected anterior and middle putamen, with relatively preserved tracer binding in the less affected posterior putamen. The relatively more preserved DAT in the posterior putamen might be a compensatory response to the more severe dopaminergic loss observed with VMAT2. While it is still debatable whether lower DAT contributes to higher synaptic DA levels by reducing DA reuptake, it has also been shown that, in PD, higher DAT is associate with lower DA turnover, i.e. the functional role of DAT may be to maintain relatively constant synaptic DA levels (Sossi et al., 2007). The substantia nigra and thalamus also appeared relatively more preserved with MP, similar to the common MP pattern along the first pair of canonical variates. DAT is found along the projections from the substantia nigra to the thalamus and then to the striatum. The globus pallidus receives input from the thalamus, and the fact that this region appeared relatively more affected for DTBZ compared to MP may again suggest a possible compensatory role effect of DAT in the substantia-thalamus-globus pallidus/striatum pathway in early disease. The downregulation of VMAT2 in VTA may be a reflective of the fact that the mesolimbic pathway is also affected relatively early in PD (Bosboom et al., 2004). However, binding in extrastriatal regions may not be specific to dopaminergic neurons and more accurate interpretations of the functional roles of the extrastriatal regions in these spatial patterns require more detailed studies involving healthy controls.

Limitations

There are several limitations in this study. First of all, in order to unambiguously determine if patterns are related to disease or normal topology differences in the two tracers, the same analysis should be further extended to DTBZ and MP data from healthy controls. However, since all common patterns in both DTBZ and MP highly resembled known characteristic disease-related dopaminergic changes, we believe these patterns are indeed related to disease. As an indirect comparison of spatial patterns in the disease and healthy stages, we applied PCA to DTBZ data in PD and healthy control groups then compared the PCA patterns obtained in the PD group and healthy control group with the CCA patterns presented in the paper. Detailed results of the comparison can be found in the Supplementary Materials. Secondly, both tracers are not entirely selective for dopaminergic neurons (DTBZ is taken up by monoaminergic terminals and MP is not 100% specific for DAT). Although such contribution is very small in the striatum in healthy condition, it may not be entirely negligible in PD or in other regions. While these observations may introduce a potential confound in the interpretation of the data, such confounds are not specific to this analysis method but to any approach comparing MP and DTBZ data. In addition, we applied the proposed method to ROI BPND values instead of parametric BPND maps to reduce the effect of noise. Another important limitation of the study is the relatively small sample size. Results from this study will be further confirmed with larger sample size and parametric images in the future.

Conclusion

Using two extensively used tracers, we showed that the proposed joint pattern analysis approach was able to capture all disease-induced characteristic spatial and temporal distribution patterns with better sensitivity compared to univariate analysis. This approach can be easily extended to the analysis of a larger number of data sets and thus appears very well suited to the analysis of multiple data sets, multi-tracer or multi-modality. It can be further extended to include voxel-level data. The method has several advantages in terms of biologically-relevant information that can be extracted from the data: first, it considers tracer distributions in all ROIs at once, thus providing information not only on localized alterations, but also on spatial patterns of such alterations, emphasizing an network behavior of the targets under investigation. Secondly, the approach decomposes the common information between data sets, in our case DTBZ and MP binding, into distinct orthogonal patterns of characteristic dopaminergic changes that are either more sensitive to disease discrimination or to disease progression and potentially resulting from somewhat independent underlying mechanisms. Thirdly, it allowed to identify unique behavior of each specific target and thus possibly discern the relative target response to disease. While the data considered in this study allowed to validate this approach, application of this method to a larger data set, including healthy controls and/or patients with more advanced disease and/or other tracers, is expected to provide new insights into the effect of disease on multiple targets, their interaction and behavior as a function of disease progression in an entirely data driven manner. Extension of the method to voxel level data and other atypical parkinsonisms might also be of interest in the future.