Participants
Participants were recruited via online advertising. All participants gave informed written consent. Inclusion criteria were age above 18 years and capacity to give written informed consent. Exclusion criteria were (i) any past or current major medical condition, (ii) history of a neurological or psychiatric disorder (including substance abuse/dependence) as determined by the Structured Clinical Interview for Diagnostic and Statistical Manual of Mental Disorders (DSM-IV) and medical review, (iii) history of head injury with a loss of consciousness, (iv) a family history of any psychiatric disorder in first-degree relatives, and (v) contraindications to PET or MRI scanning (significant prior exposure to radiation, pregnancy, or breast feeding) or amphetamine.
Participants attended for a scan on two separate occasions, separated by a minimum of 3 days. On one occasion, participants received an oral placebo (101 mg of lactose/sucrose tablets; number of tablets matched to number administered when receiving dexamphetamine) before scanning, and on the other, they received a dose (0.5 mg/kg) of oral dexamphetamine. The order of drug/placebo scans was randomized. Dosing occurred 3 hours before scanning so that in the case of dexamphetamine, peak drug levels coincided with scan time (
45). The order of scans (i.e., dexamphetamine/placebo order) was randomized, and participants and staff involved in assessment were blinded to this order. The subjective effects of amphetamine scale (
23,
24) was given to participants at baseline for 1.5 and 3 hours to assess the psychological effects of amphetamine. The scale involves ranking each of 10 possible drug effects on a five-point scale ranging from “least” to “most.” Positive effects include “high,” “rush,” “good effects,” “liking,” and “desire for drug”; negative effects included “fidgety,” “anxious,” “dizziness,” “dry mouth,” and “distrust.” Measures obtained at the 1.5-hour time point were used in subsequent analysis, as this was the time point with the most data obtained and at which subjective effects were greatest. The study was approved by the local National Health Service (NHS) Research Ethics Committee (12/LO/1955) and the Administration of Radioactive Substances Advisory Committee.
Image acquisition
Neuroimaging data were acquired using a GE SIGNA simultaneous PET-MRI scanner. [11C]-(+)-PHNO (0.020 to 0.029 μg/kg) was injected immediately before the commencement of scanning as a smooth bolus injection over 30 s. During the initial minutes of data acquisition, the PET signal changes rapidly, and PET modeling will be particularly sensitive to the integrity of the data. The heat generated by the gradient coils during MRI acquisition was noted to have a small effect upon the PET signal, so a 10-min MRI-free period at scan start was used. This initial 10-min period of PET-only data collection was followed by simultaneous PET-MRI acquisition. Mean injected activity was 140 megabecquerel.
Structural MRI
A three-dimensional (3D) BRAVO T1-weighted (T1w) structural scan was obtained with the following parameters: flip angle = 12°, inversion time (TI) = 400 ms, echo time (TE) = 3.2 ms, repetition time (TR) = 8.5 ms, matrix = 256 × 256, number of slices = 188, and 1-mm isotropic voxels.
Arterial spin labeling
Images were acquired using a 3D pseudo-continuous ASL sequence, and 38 slices of a 128 × 128 matrix were obtained, resulting in a spatial resolution of 1.88 mm by 1.88 mm by 4 mm. Parameters were as follows: TE/TR = 10.7/4854 ms, labeling duration = 1450 ms, postlabeling delay = 2025 ms, interslice gap = 0, bandwidth = 62.5 kHz, flip angle = 111°, and number of acquisitions = 4. A proton density image was acquired using the same parameters to compute the CBF map in standard physiological units (milliliters of blood per 100 g of tissue/min).
Automatic CBF map reconstruction used the following quantification algorithm
T1b is the
T1 of blood (1.6 s at 3 T). The partial saturation of the proton density image (PD) is corrected using a typical gray matter value for
T1t of 1.2 s. The partition coefficient λ is set to the whole-brain average of 0.9. Efficiency, ε, is set to 0.6. PLD refers to the postlabeling delay, and labelling time (LT) refers to the labeling duration. PW is the perfusion-weighted image. SF
PW is the scaling factor of the sequence, and NEX
PW is the number of excitations. As only a final label-control difference image is produced by the scanner, motion correction cannot be applied retrospectively. However, this should not markedly affect the data quality because background suppression and the alternation between control and labeled phases reduce motion artifact [see (
46) for full details].
Positron emission tomography
Dynamic emission data were acquired for 90 min following radiotracer administration. A zero echo time (ZTE)–based MR attenuation correction map with the following parameters was obtained: flip angle = 0.8°, matrix = 110 × 110 × 116, voxel size = 2.4-mm isotropic, number of averages = 4, bandwidth = ±62.5 kHz, and acquisition time = 42 s.
Dynamic PET images were reconstructed on the PET/MR scanner using VPFX-S, a fully 3D ordered subset expectation maximization algorithm with time-of-flight information and resolution recovery, six iterations, 16 subsets, no postreconstruction smoothing, matrix of 128 × 128 × 89, and voxel size of 2 mm by 2 mm by 2.78 mm, with corrections applied for detector normalization, randoms, scatter, dead time, and radioactive decay. Attenuation correction was performed using the ZTE map implemented using GE scanner software (
47).
Canonical correlation analysis
CCA is a technique ideally suited to the simultaneous evaluation of two highly multidimensional sets of variables. CCA allows the identification of many-to-many relations in contrast to techniques such as multiple linear regression that requires one set of variables to be condensed to a single variable. For an intuitive understanding, CCA can be understood as an extension of principal components analysis (PCA). PCA decomposes a single set of variables to a smaller number of components that capture the variables’ latent sources of variation, while CCA seeks to simultaneously decompose two sets of variables to maximize the correlation between these two sets (see
Fig. 1) (
14). For a full description of the technique, with an emphasis on applications in neuroimaging, please see the review by Wang
et al. (
14).
For both CBF data (229,007 voxels × N) and PET data (2847 voxels × N), where there are data from N = 51 scans, dimensionality reduction via PCA was undertaken before running CCA (see figs. S1 and S2) so as to avoid the overfitting that occurs if the number of features is significantly greater than the number of subjects. As PCA is affected by the relative scaling of each variable, standardization of each variable was performed by subtracting its mean and dividing by its SD. The optimal number of components is unclear, and PCA was therefore performed to give between K ∈ [2,3,4, …,10] components for both CBF and PET measures. This results in nine pairs of [N, K] matrices C and P, where C contains features derived from the CBF data and P contains features derived from the PET data. For each pair, the number of rows is equal to the number of scans (N = 51), while the number of columns is equal to the number of PCA components. For the PET data, the first 10 components explained 61% of the variance, while for CBF data, 42% were explained.
CCA was implemented using the python library “scikit-learn” (
61). CCA calculates a pair of [1,
K] canonical weight vectors
U and
V (which are each the length of the number of components), which maximize the correlation between
PU and
CV.
PU and
CV are vectors termed canonical variates, and each will be the length of the total number of scans. A pair of canonical variates is termed a mode of covariation. Given that CCA intentionally maximizes the correlation between canonical variates, statistical significance cannot be determined solely on the basis of the magnitude of a correlation coefficient between these canonical variates.
An out-of-sample cross-validation testing approach was used to test the statistical significance of the mode of covariation identified between striatal dopamine receptor availability and CBF. For each choice of PCA dimensionality, K, the data were randomly split 100 times into sets where 80% of the data were used for training and 20% for testing. For each train-test split, we conducted PCA with K components followed by CCA on the training data. The PCA and CCA models that were fit on the training data were then applied to the test data, and the correlation between the resulting test data canonical variates for the first CCA component was calculated. We used the mean correlation coefficient over 100 train-test splits as an estimate of the true out-of-sample performance for that PCA dimensionality. An overall correlation coefficient was then calculated by taking the mean correlation coefficient for the nine (K ∈[2,3,4, …,10]) dimensionalities. We used this method to avoid potential biases that could result from manually setting the number of PCA dimensions. The statistical significance of this correlation coefficient was calculated by comparing it to a null distribution generated by following the same procedure 10,000 times but after permuting CBF images to shuffle the mapping between CBF images and participant identity (i.e., instantiating the null that there is no shared information between striatal D2/3R availability and cortical CBF within subjects).
For most of the participants, there exist two scans (one following placebo, and the other following amphetamine). There, therefore, exists the potential for information leakage between training and test sets, which could erroneously inflate the statistical significance of results. We used two strategies to address this. First, we used a well-established approach in which the same dependency structure was kept when generating a null distribution (
19). As a result, while the observed out-of-sample correlation coefficient may be inflated secondary to dependency structure, an equal inflation will also apply to the null distribution and therefore the
P value will be unaffected. Second, in a complementary analysis, we restricted the shuffling of participants when partitioning the data into train and test sets so that placebo and amphetamine scans for a single subject remain together for each subject and are not separated across test and train sets, meaning that no dependency exists between test and train sets. To ensure that the observed relationship reflected cortex-wide relationships and was not driven by subcortical changes in blood flow, we also performed the above analysis while restricting the analysis to CBF data from cortical regions defined using the Harvard-Oxford cortical atlas (
62).
In a complementary analysis, CCA (and PCA) was performed as on a training sample that randomly left out two pairs of PET-ASL scans. This model was then used to perform CCA on the two pairs of left out test scans. If the CCA scores for CBF and PET were ranked identically, e.g., the first pair of scans having higher scores for both CBF and PET, then this indicates that the CCA has been able to correctly match the pairs of scans and that this was counted as a successful match. If the ranking differed between modalities, then this indicated an unsuccessful match. Accuracy was quantified as the percentage of successful matches out of 100 samples. This was compared to null distributions generated as above by both randomly selecting subjects and selections in which the pair structure was maintained.
After establishing a statistically significant relationship between cortical activity and striatal dopamine, we then sought to investigate the nature of this relationship in greater depth. To do this, we first set the number of PCA components to eight as that had gave the most predictive out of sample relationship between the PET and CBF measure. We then fit CCA on the entire sample (resulting full loading matrices in table S1) and calculated the contribution of each voxel to by correlating CBF and PET measures with their respective canonical variate at each voxel (i.e., CBF voxels to the CBF canonical variate and PET voxels to the PET canonical variate). This is illustrated by the maps of covariation influence in
Fig. 3 (C and D), where the intensity of shading represents the magnitude of correlation coefficient.
Gene expression data
Gene expression data were obtained from six post mortem adult brains from the AHBA transcriptomic dataset (
http://human.brain-map.org; RRID: SCR_007416) (
63). Since the AHBA only includes data relating to two subjects for the right hemisphere, right-sided striatal samples were limited, and in keeping with similar analyses, we only consider left hemisphere samples (
36).
The AHBA microarray expression data need to be both linked to regions of interest and combined across donors. To ensure robust and replicable results, we used the “abagen” (version 0.0.5) python toolbox, an externally developed and validated pipeline, which undertakes the steps described below (
64).
The original MNI coordinates provided with the AHBA do not account for nonlinear deformations in their transformation from native to MNI space. We therefore used “corrected” MNI coordinates to match tissue samples to MNI space; this involves using ANTs to perform a more accurate nonlinear coregistration of each donor’s scan to an MNI template as opposed to the AHBA default based on an affine coregistration of a single donor brain.
Probe-to-gene mappings were reannotated with information from Arnatkevic̆iūtė
et al. (
64). This was followed by intensity-based filtering to remove probes that did not exceed background noise in 50% of all tissue samples. A representative probe was selected for multiple probes indexing the same gene based on which probe had the highest differential stability among donors. Each sample was matched to the voxel in which it lay and, if lying outside the striatal mask, was matched to the closest voxel as long as this was within 3 mm. This resulted in 153 eligible samples. Before aggregation across donors, expression values were normalized for each sample across genes for each donor using a scaled robust sigmoid normalization function.
Gene enrichment analysis
The contribution of each gene to the striatal covariation influence was calculated on the basis of that gene’s VIP score. The error in estimating each gene’s VIP score was assessed by bootstrapping (resampling with replacement of the 153 striatal voxels), z scores were then calculated as the ratio of each gene’s VIP to its bootstrap standard error, and this was used to rank the genes according to their level of influence on the striatal covariation influence map.
GO enrichment analysis was used to characterize the nature of the genes associated with the striatal covariation influence. To avoid artefactual elevation of GO terms for brain expression, we used previously reported methods (
35) to first subset the VIP ranked gene list to include only genes showing above zero expression in cerebral cortex, amygdala, basal ganglia, cerebellum, cerebral cortex, hippocampal formation, and hypothalamus as listed within Human Protein Atlas RNA GTEx brain region gene expression database (
www.proteinatlas.org/download/rna_brain_gtex.tsv.zip; table S3). Enrichment analysis was then performed for this brain refined, and VIP score ordered list was inputted to GOrilla (separately in increasing and decreasing order) (
28). Resulting gene categories were considered significant if meeting a threshold of false discovery rate (FDR)–corrected
P < 0.05 and showed a minimum of three overlapping genes. REViGO (
http://revigo.irb.hr) was then used to cluster categories using the SimREl algorithm (
66).
Using previously reported methods (
35,
36), we next explored whether the VIP ranked gene list was enriched for genes associated with neuropsychiatric illnesses. We used the findings of the psychENCODE consortium to identify genes that are up-regulated or down-regulated in autism spectrum disorder, schizophrenia, and bipolar disorder (
37) and also genes differentially expressed in temporal lobe epilepsy (
40) and Alzheimer’s disease (
41). We first calculated the median rank of the disease-associated gene set within the VIP gene list and took the absolute value of how far this lay from the center of the list. We then compared this value to 10,000 randomly selected gene sets of the same size.