Mapping dopaminergic projections in the human brain with resting-state fMRI

The striatum receives dense dopaminergic projections, making it a key region of the dopaminergic system. Its dysfunction has been implicated in various conditions including Parkinson’s disease (PD) and substance use disorder. However, the investigation of dopamine-specific functioning in humans is problematic as current MRI approaches are unable to differentiate between dopaminergic and other projections. Here, we demonstrate that ‘connectopic mapping’ – a novel approach for characterizing fine-grained, overlapping modes of functional connectivity – can be used to map dopaminergic projections in striatum. We applied connectopic mapping to resting-state functional MRI data of the Human Connectome Project (population cohort; N = 839) and selected the second-order striatal connectivity mode for further analyses. We first validated its specificity to dopaminergic projections by demonstrating a high spatial correlation (r = 0.884) with dopamine transporter availability – a marker of dopaminergic projections – derived from DaT SPECT scans of 209 healthy controls. Next, we obtained the subject-specific second-order modes from 20 controls and 39 PD patients scanned under placebo and under dopamine replacement therapy (L-DOPA), and show that our proposed dopaminergic marker tracks PD diagnosis, symptom severity, and sensitivity to L-DOPA. Finally, across 30 daily alcohol users and 38 daily smokers, we establish strong associations with self-reported alcohol and nicotine use. Our findings provide evidence that the second-order mode of functional connectivity in striatum maps onto dopaminergic projections, tracks inter-individual differences in PD symptom severity and L-DOPA sensitivity, and exhibits strong associations with levels of nicotine and alcohol use, thereby offering a new biomarker for dopamine-related (dys)function in the human brain.


Introduction
The brain's dopamine system plays an important role in a wide range of behavioural and cognitive functions, including movement and reward processing (Joshua et al., 2009;Ruhé et al., 2007). An integral structure of the dopamine system is the striatum, which receives dense dopaminergic projections from the substantia nigra pars compacta (SNc) and ventral tegmental area (VTA) in the midbrain (Steiner and Tseng, 2016). Work in experimental animals has shown that these projections organize along a gradient: dopaminergic neurons in the SNc project preferentially to dorsal caudate and putamen in dorsolateral striatum, while dopaminergic neurons in the VTA project predominantly to the nucleus accumbens (NAcc) in ventromedial striatum (Steiner and Tseng, 2016;Haber, 2014;Björklund and Dunnett, 2007). The projections from the SNc to dorsolateral striatum comprise the nigrostriatal pathway implicated in, for example, the organization of motor planning (Joshua et al., 2009;Faure et al., 2005). The mesolimbic pathway formed by the projections from the VTA to the NAcc has been associated with reward processing (Schultz, 2013;Wise, 2004). In accordance with the partial neuroanatomical overlap in striatum, increasing evidence also suggests partial overlap in the function of both pathways (Haber et al., 2000;Everitt and Robbins, 2005;Wise, 2009). Of note, dopaminergic neurons in the VTA not only project to NAcc but also to prefrontal cortex. These cortical projections form the mesocortical pathway associated with reward-related goal-directed behaviours (Schultz, 2013;Wise, 2004).
In humans, alterations in these dopaminergic projections have been associated with multiple neurological and psychiatric conditions (DeLong and Wichmann, 2007;Money and Stanwood, 2013). A well-known example is Parkinson's disease (PD), a neurodegenerative disorder characterized by a loss of dopaminergic neurons in the SNc (part of the nigrostriatal pathway; Fearnley and Lees, 1991), which frequently causes asymmetric depletion of dopamine in dorsal striatum (first in putamen, later also to a lesser extent in caudate) and leads to impairments in motor as well as a range of nonmotor functions (Brooks and Piccini, 2006;Hornykiewicz, 2008). Dopaminergic dysfunction has also been implicated in substance use disorder given that addictive substances, such as stimulants, alcohol, and nicotine, increase the release of dopamine in ventral striatum (i.e., mesolimbic pathway; Laruelle et al., 1995;Barrett et al., 2004;Nutt et al., 2015).
Despite the important role of the dopamine system in human brain function and its implication in disease, knowledge about this neurotransmitter system is limited and mainly based on experimental work in animals. The investigation of dopaminergic functioning in vivo in the human brain is challenging, although the nuclear imaging techniques position emission tomography (PET) and single photon emission computed tomography (SPECT) can be used for this purpose (Blake et al., 2003;Volkow et al., 1996). Imaging of the density of the dopamine transporter (DaT) using SPECT has become a popular tool to assist in the differential diagnosis of PD as loss of dopaminergic neurons in PD is accompanied by a loss in DaT in striatum, as opposed to lookalike conditions such as dystonic tremor where the DaT signal remains intact (Brooks, 2016). Tracking the loss of DaT signal over time has also been proposed as a progression biomarker for PD (Brooks, 2016). Indeed, DaT reuptakes dopamine from the synaptic cleft after its release and is highly expressed in the terminals of dopaminergic neurons projecting from the midbrain to striatum (Brooks, 2016). Therefore, DaT SPECT imaging can be used to image dopaminergic projections in striatum. However, the radiation exposure and costs of PET/SPECT combined with the low spatial resolution of the scan limit widespread implementation in human brain research and in clinical practice.
In this work, we hypothesize that inter-individual differences in DaT availability induce interindividual variations in the synchronicity of functional activity in the brain, and therefore, that dopaminergic projections in the human striatum can also be mapped using blood-oxygen-level-dependent (BOLD) functional MRI (fMRI) measured at rest. We employ a 'connectopic mapping' data analysis approach to disentangle striatal connectivity into multiple overlapping spatial 'modes' in order to dissect the complex mixture of efferent and afferent connections of the striatum to multiple cortical and subcortical systems (that map onto different neurobiological systems and associated functions; Haak et al., 2018). In previous work, we already showed that the dominant (zeroth-order) mode represents its basic anatomical subdivisions, while the first-order mode maps on to a ventromedialto-dorsolateral gradient associated with goal-directed behaviour in cortex (Marquand et al., 2017) that has been described previously on the basis of tract-tracing work in non-human primates (Haber and Knutson, 2010). Here, we demonstrate -by conducting a series of analyses across different datasets -that the second-order mode of gradual spatial variations in the BOLD connectivity pattern reflects DaT availability in the striatum. We furthermore reveal that this mode tracks inter-individual differences in symptom severity in PD patients, is sensitive to acute dopaminergic modulation (L-DOPA administration), and exhibits strong associations with levels of nicotine and alcohol use in a population-based cohort. Hereby, we provide compelling evidence that this connectivity mode tracks inter-individual differences in dopaminergic projections, and as such, offers a new biomarker for investigating dopamine-related dysfunction across various neurological and psychiatric disorders.

Striatal connection topographies map onto DaT availability
For our first analysis, we applied connectopic mapping (Haak et al., 2018) to resting-state fMRI data from 839 participants of the Human Connectome Project (HCP; Van Essen et al., 2013). Connectopic mapping extracts the dominant modes of functional connectivity change (or connection topographies) within the striatum based on a Laplacian eigenmap decomposition of the similarity matrix derived from functional connectivity (i.e., Pearson correlations) computed between each striatal voxel and the rest of the brain. It provides reproducible and parsimonious representations of overlapping connection topographies at both the group level and at the level of individual subjects. The connectopic mapping approach is detailed in Materials and methods, but a summary of this procedure can be found in Figure 1.
For all analyses described in this paper, connectopic mapping was applied to the left and right putamen and caudate-NAcc subregions separately to increase regional specificity and the secondorder striatal connectivity mode was selected for each of the four striatal regions of interest (ROIs). A spatial statistical model, that is, a trend surface model (TSM; Gelfand et al., 2010), was fitted to both the group-level and the subject-specific connectivity modes to obtain a small set of coefficients summarizing each of the four striatal modes in the X, Y, and Z axes of MNI152 coordinate space, which we used for statistical analyses. A scree test (Cattell, 1966) indicated that a quadratic model (i.e., consisting of six TSM coefficients) provided the best fit for the second-order connectivity mode in putamen and a quartic model (12 TSM coefficients) was found to provide the best fit for the secondorder connectivity mode in caudate-NAcc region.
The subject-specific second-order striatal connectivity modes were highly consistent across the two fMRI sessions (mean ± SD: ρ = 0.98 ± 0.07; averaged across all four subregions) of the HCP dataset, which is in line with what we have demonstrated previously for other brain regions and for the zerothorder and first-order mode of connectivity in striatum. Furthermore, interclass correlation (ICC(2,k)), which indexes measurement consistency for a putative biomarker (Shrout and Fleiss, 1979;Koo and Li, 2016), showed excellent reproducibility of the subject-specific connectivity modes, while still being sensitive to inter-individual differences (see Table 1). Both the variations across subjects and the reproducibility within subjects are illustrated in Figure 2-figure supplement 3.
The group-level second-order connectivity mode across striatum is displayed in Figure 2 (second row). The modes for left and right putamen and caudate-NAcc have been combined in this figure (i.e., the four ROIs were loaded in FslView simultaneously from which the below figures were derived) to aid in visualization and for later comparison to the DaT SPECT scan. The second-order connectivity mode comprises a gradient from the dorsal putamen and dorsal caudate (shown in red) to the ventral putamen and ventral caudate including the NAcc (shown in blue). This coding indicates that the dorsal putamen and dorsal caudate exhibit a connectivity pattern with the rest of the brain that is similar to each other but different from the ventral putamen and ventral caudate and vice versa. This striatal connectivity pattern might thus correspond with the gradient of mesolimbic and nigrostriatal dopaminergic projections to striatum (ventral vs. dorsal striatum) well described by track-tracing studies in rodents and non-human primates (Steiner and Tseng, 2016;Haber, 2014;Björklund and Dunnett, 2007). We therefore investigated its spatial correspondence to DaT SPECT-derived DaT availability in striatum, which is assumed to be an index of dopaminergic projections. To this end, we averaged across DaT SPECT images obtained from 209 healthy control participants from the Parkinson's Progression Markers Initiative (PPMI) dataset (Marek et al., 2011). As can be observed in Figure 2, the group-level second-order striatal connectivity mode indeed displays a remarkably high similarity with the group-level DaT availability in striatum, as quantified by a spatial voxel-wise correlation of r = 0.884 (p<0.001), thereby providing the first evidence for an fMRI-derived striatal connectivity marker strongly associated with dopaminergic projections into striatum. A high correlation (r = 0.925, p<0.001) is also present between the orthogonal TSM coefficients modelling the group-level second-order connectivity mode and the group-level DaT SPECT scan across the striatum, providing more evidence that the second-order connectivity mode maps onto dopaminergic projections. This finding does not strongly depend on the chosen model order, given that repeating this analysis using model order 3 (i.e., a cubic model with nine TSM coefficients for both the putamen and caudate-NAcc regions) resulted in a similar correlation (r = 0.90, p<0.0001). The connectopic mapping pipeline. The functional MRI (fMRI) time-series data from a predefined region of interest (ROI), here the striatum, are rearranged into a time-by-voxels matrix A, as are the time series from all voxels outside the ROI (matrix B). For reasons of computational tractability, the dimensionality of B is losslessly reduced using singular value decomposition (SVD), yielding ∼B. For every voxel within the ROI, its connectivity fingerprint is computed as the Pearson's correlation (CORR) between the voxel-wise time-series and the SVD-transformed data, yielding matrix C. Then similarity between voxels is computed using the η 2 coefficient, resulting in matrix S. Manifold learning using Laplacian eigenmaps is then applied to this matrix, yielding a set of overlapping, but independent, connection topographies or 'connectivity modes' that together describe the functional organization of the striatum. These connection topographies indicate how the connectivity profile with the rest of the brain changes across striatum. Voxels that have similar colours in these connectivity modes have similar connectivity patterns with the rest of the brain. Finally, trend surface modelling is applied to summarize the connectivity modes by fitting a set of trend coefficients (β) that optimally combine a set of spatial polynomial basis functions. See Haak et al., 2018 for further details.
Finally, to further demonstrate the high specificity of the second-order connectivity mode to the DaT SPECT scan, we computed correlations with the TSM coefficients of all PET scans, tapping into various neurotransmitter systems, included in the publicly available JuSpace toolbox (Dukart et al., 2021). Figure 2-figure supplement 1 reveals that the correlation between the TSM coefficients of the second-order connectivity mode with the DaT SPECT scan is not only highly significant but also significantly higher than the correlations with the TSM coefficients of any other PET scan.
In order to demonstrate that also individual variations in this connectivity mode are associated with individual variations in striatal dopaminergic projections, we further aimed to replicate this mapping at the within-subject level in a subsample of PPMI participants (130 datasets from PD patients and 14 from controls) with both DaT SPECT and resting-state fMRI data available. Within a smaller sample of PD patients and controls with good quality connectivity modes (see Appendix 1-Supplementary analyses and Figure 3-figure supplement 1 for further details), we not only replicated the spatial correspondence between the connectivity mode and DaT SPECT scan at the group level (PD group: r = 0.714; control group: r = 0.721) but also observed a within-subject spatial correlation of 0.58 across the four striatal subregions (0.44> r < 0.62; mean = 0.58, 95% CI = [0.56,0.60]) (see Figure 3). These findings were not induced by residual head motion (see

Striatal connection topographies are altered in PD
The strong association of the second-order striatal connectivity pattern with DaT availability suggests that this resting-state fMRI-derived connectivity mode can be used to assess variability (including disease-related alterations) in dopaminergic projections to the striatum. As such, we hypothesized that the second-order connectivity mode would be altered in PD since this disorder is characterized by progressive degeneration of nigrostriatal dopaminergic neurons. In order to validate this hypothesis, we made use of a separate high-resolution PD dataset (Dirkx et al., 2019) including 39 PD patients (19 patients with asymmetric symptoms on the left side of the body, i.e., they were left-dominant; 20 patients were right-dominant) and 20 controls that each underwent two high-resolution restingstate fMRI session (T = 860 ms, 700 time points). Participant characteristics can be found in Table 2. During one session, PD patients received dispersible 200/50 mg levodopa-benserazide (L-DOPA), a precursor of dopamine used for the treatment of PD, during the other session they received placebo (dispersible cellulose). Controls did not receive L-DOPA and placebo but just underwent two typical resting-state fMRI sessions under the same scanning protocol. Thus, this dataset did not only allow us to investigate the effects of clinical diagnosis, but also the effects of acute dopaminergic modulation on the underlying striatal connectivity mode. The associations with diagnosis and L-DOPA were investigated in both groups separately (left-dominant, right-dominant PD) because the side of predominant nigrostriatal dopamine depletion likely influences the pattern of striatal connectivity. As before, the second-order striatal connectivity mode was modelled separately for the putamen and caudate-NAcc subregions to increase regional specificity as PD is known to affect the putamen region of the striatum before the caudate-NAcc region (Kish et al., 1988). Group differences in the TSM coefficients modelling the putamen and caudate-NAcc subregions were subsequently assessed by conducting an omnibus test of all the TSM coefficients, that is, a likelihood ratio test in the context of a logistic Table 1. Interclass correlation coefficients (ICCs) between the two scanning sessions and the session 1 to session 2 within-subject and between-subject spatial correlations.     regression. We applied correction for multiple comparisons (two groups: left-and right-dominant PD × 2 striatal subregions: putamen and caudate-NAcc) using a Bonferroni-corrected α-level of 0.05/4 = 0.0125. In addition, we also investigated associations between the TSM coefficients and symptom severity across PD patients. To this end, we fitted general linear models (GLMs) that included the TSM coefficients modelling the gradient during the placebo session to predict the total score on the motor section (part III) of the Movement Disorder Society Unified Parkinson's Disease Rating Scale (UPDRS; Goetz et al., 2008). This analysis was again conducted separately for the left-and right-dominant . Within-subject correlations between the second-order connectivity mode and the DaT SPECT scan from subjects in the Parkinson's Progression Markers Initiative (PPMI) cohort where both resting-state functional MRI (fMRI) and DaT SPECT data is available. These correlations were obtained in a subsample of the PPMI dataset (6-8 datasets from controls and 73-82 datasets from Parkinson's disease patients depending on the striatal subregion) with good-quality connectivity modes as defined by a high spatial correlation (r > 0.5) with the group-average Human Connectome Project (HCP) connectivity mode. Red dots represent control participants; black dots represent Parkinson's disease patients.
The online version of this article includes the following figure supplement(s) for figure 3:  groups and separately for putamen and caudate-NAcc subregions, and a Bonferroni-corrected α-level of 0.0125 was used for establishing statistical significance. The PD (placebo session) vs. control group analysis (first session) of the TSM coefficients modelling the second-order striatal connectivity mode revealed a significant difference between the rightdominant PD group and the control group in bilateral putamen (but not in the caudate-NAcc region) of the striatum (see Figure 4A; omnibus test of all TSM coefficients: X 2 = 27.17, p=0.007). No significant differences were observed between the left-dominant PD group and the control group. Moreover, within the right-dominant patient group, we observed a trend-level association between UPDRS symptom severity scores and the TSM coefficients modelling the putamen under placebo (GLM omnibus test of the TSM coefficients: X 2 = 22.28, p=0.035). Post-hoc Pearson correlations revealed that this effect was driven by the quadratic TSM coefficients modelling the striatal connectivity mode in the right putamen in the Y (i.e., anterior-posterior) direction (right Y 2 : r = 0.476, p=0.034) and Z (i.e., superior-inferior) direction (right Z 2 : r = −0.460, p=0.041). This association can be observed in Figure 4B, which shows an increase in blue-coded voxels in the second-order striatal connectivity mode as symptom severity increases in PD. This pattern maps very well on the observed decrease in dopaminergic projections as PD becomes clinically more severe, as reflected by higher UPDRS scores. That is, given the spatial similarity of the second-order striatal connectivity mode with the DaT SPECT scan, we can interpret the observed alteration in the connection topography as a decrease in dopaminergic projections to striatum. In Figure 4B, this is evident as an increase in blue-coded voxels and a decrease in red-coded voxels as a function of UPDRS symptom severity. Supplementary analyses showed that the observed group differences and associations with symptom severity were independent of age and sex (Appendix 2-table 1).

Striatal connection topographies are sensitive to acute dopaminergic modulation
After demonstrating that the second-order connectivity mode in putamen is indeed altered in PD, we investigated whether it was also sensitive to the acute effects of the dopamine precursor L-DOPA. We assessed differences in this striatal connectivity mode between the placebo and L-DOPA session in both PD groups by conducting an omnibus test of all the TSM coefficients, that is, a likelihood ratio test in the context of a logistic regression. We applied correction for multiple comparisons (two groups: left-and right-dominant PD × 2 striatal subregions: putamen and caudate-NAcc) using a Bonferroni-corrected α-level of 0.05/4 = 0.0125. These tests did not reveal significant differences between the placebo and L-DOPA session in the putamen or the caudate-NAcc region. However, For the FAB, lower scores indicate worse functioning; for the UPDRS, higher scores indicate worse functioning. The FAB score was evaluated off medication. FAB = frontal assessment battery (score 0-18); UPDRS = Unified Parkinson's Disease Rating Scale part III (score 0-132); NS = not significant; NA = not applicable; PD = Parkinson's disease; L-DOPA = levodopa-benserazide.
The online version of this article includes the following source data for table 2: Source data 1. Source data for participant characterists listed in Table 2.
treatment response to L-DOPA is known to differ among PD patients. To take this variability across patients into account, we conducted GLM analyses relating differences in the L-DOPA-induced change (difference between L-DOPA and placebo session) in the second-order striatal connectivity mode to differences in treatment response. These analyses showed that in both patient groups significant associations were present between the L-DOPA-induced change in UPDRS symptom severity scores and the L-DOPA-induced change in TSM coefficients in putamen (GLM omnibus test of all TSM coefficients in the right-dominant patient group: X 2 = 25.48, p=0.012; in the left-dominant patient group: Post-hoc Pearson correlations revealed that this effect was driven by the quadratic TSM coefficients modelling the striatal connectivity mode in the right putamen in the Y (i.e., anterior-posterior) direction (right Y 2 : r = 0.476, p=0.034) and Z (i.e., superior-inferior) direction (right Z 2 : r = −0.460, p=0.041). The correlation between UPDRS symptom severity scores is displayed for the right Y 2 coefficient. To visualize this association, the reconstructed second-order connectivity mode in the right putamen is shown for five Parkinson's disease patients (data points circled) with increasing UPDRS symptom severity scores.
X 2 = 34.07, p=0.001). Post-hoc Pearson correlations revealed that these effects were driven by the linear TSM coefficient (Y 1 ) modelling the striatal connectivity mode in putamen in the Y (i.e., anteriorposterior) direction (right-dominant PD: left Y 1 : r = −0.548, p=0.012; left-dominant PD: right Y 1 : r = 0.345, p=0.15). As can be seen in Figure 5, a larger L-DOPA-induced reduction in UPDRS scores is associated with a larger positive change (i.e., an increase in red-coded voxels) in the superior-anterior part of the putamen, which we hypothesize maps onto an increase in dopamine-related connectivity. Supplementary analyses demonstrated that the observed effects of L-DOPA were independent of age and sex (see Appendix 2-table 1).

Striatal connection topographies are associated with the amount of substance use
Finally, dopaminergic signalling is also implicated in reward processing, alterations of which have been associated with substance use and compulsive behaviours (Laruelle et al., 1995;Barrett et al., 2004;Nutt et al., 2015). As such, we investigated whether the second-order striatal connectivity mode was also associated with tobacco and alcohol use. These quantities are amongst the set of demographic variates available from the HCP. In order to increase specificity and also in order to transcend analysis from a categorical comparison to a continuous characterization predicting the relative amount of substance usage, we consider here a subset of high tobacco and alcohol users. Specifically, we selected otherwise drug-naïve HCP participants who reported to have consumed ≥3 light and/or ≥ 1 heavy alcoholic units per day during the week preceding the scan (N = 30) and participants reporting to have smoked ≥5 cigarettes every day during the week preceding the scan (N = 38). GLM analyses investigating associations between the TSM coefficients modelling the second-order striatal connectivity mode and the amount of use over the past 7 days were conducted separately for the alcohol users and smokers and separately for putamen and caudate-NAcc subregions. We applied correction To visualize this association, the difference in the reconstructed second-order connectivity modes between the placebo and L-DOPA session is shown for the left putamen (at slice X = -24) for five PD patients (data points circled). Red-coded voxels are hypothesized to map onto an increase in dopaminergic connectivity, blue-coded voxels onto a decrease. A significant effect was also observed for the left-dominant PD group (GLM omnibus test of all TSM coefficients: X 2 = 34.07, p=0.001), but as post-hoc Pearson correlations did not reveal significant associations with one of the individual TSM coefficients in this group, this association is not shown. ant, anterior; post, posterior putamen.
for multiple comparisons using a Bonferroni-corrected α-level of 0.0125 (see Materials and methods for more details). In smokers, we observed a significant association with the total number of cigarettes smoked over the past week for the TSM coefficients modelling the second-order connectivity mode in the caudate-NAcc region (X 2 = 49.55, p=0.002), but not for the putamen region. Subsequent computation of the Pearson correlations between the individual TSM coefficients and the amount of use revealed that this association was driven by multiple TSM coefficients in both the left and right caudate-NAcc (left X 3 : r = −0.409, p=0.011; right Z 1 : r = 0.408, p=0.011; right Y 3 : r = 0.367, p=0.024; right Z 3 : r = −0.451, Figure 6. The second-order mode of connectivity in striatum is associated with the amount of tobacco use (top) and alcohol use (bottom). Strong associations were observed between the trend surface model (TSM) coefficients modelling the connectivity mode in the caudate-nucleus accumbens (caudate-NAcc) region and the total amount of tobacco use as well as alcohol use over the past week (general linear model [GLM] omnibus test tobacco use: X 2 = 49.55, p=0.002; alcohol use: X 2 = 64.45, p<0.001). To visualize these relationships, Pearson correlations between one of the significant TSM coefficients and the amount of use are shown as well as the reconstructed second-order connectivity mode in the right caudate-NAcc (at slice X = 14) for four tobacco users and in the left caudate-NAcc (at slice X = -14) for four alcohol users with increasing of amounts of use (data points circled). Circles and arrows indicate where in the connectivity mode tobacco and alcohol use-related changes can be observed. Correlation plots for the other TSM coefficients can be found in  p=0.005; right Y 4 : r = 0.362, p=0.026; see Figure 6-figure supplement 1). As can be observed in Figure 6 (top panel), alterations in the second-order striatal connectivity mode are subtle but consist of an increase in inferior blue-coded voxels in the caudate-NAcc as tobacco use in this population cohort increases.
In the 'heavy' drinkers, we also observed a strong association with the total number of alcoholic drinks consumed over the past week for the TSM coefficients modelling the connectivity mode in the caudate-NAcc region (X 2 = 64.45, p<0.001), but again not for the putamen region. Subsequent computation of the Pearson correlations between the individual TSM coefficients and the amount of use revealed that this association was driven by multiple TSM coefficients in both the left and right caudate-NAcc (left X 1 : r = −0.378, p=0.039; left Y 1 : r = 0.399, p=0.029; left Y 2 : r = 0.488, p=0.006; left Z 2 : r = −0.386, p=0.035; left X 3 : r = 0.519, p=0.003; left Y 4 : r = −0.417, p=0.022; right Z 2 : r = −0.418, p=0.022; right Z 4 : r = 0.488, p=0.006; see Figure 6-figure supplement 2). Similar to tobacco use, Figure 6 (bottom panel) shows that higher levels of alcohol use are accompanied by a subtle increase in blue-coded voxels as well as a decrease in red-coded voxels in the caudate-NAcc region of the second-order connectivity mode. We argue that these subtle increases in blue-coded voxels (and decrease in red-coded voxels) in high nicotine and alcohol users map onto decreases in dopaminergic connectivity, which corresponds with reported reductions in dopamine release in striatum in patients with nicotine and alcohol dependence (Nutt et al., 2015;Balfour, 2015). Supplementary analyses revealed that the associations with the amount of tobacco and alcohol use persisted under different usage thresholds and were independent of age and sex (see Appendix 2-tables 1-3). Finally, five subjects were included in both the tobacco and alcohol use analyses, but the associations with tobacco use (X 2 = 39.40, p=0.025) and alcohol use (X 2 = 62.01, p<0.001) also remained significant after excluding these subjects.

Discussion
In this work, we provide evidence for a resting-state fMRI-derived biomarker of dopamine function in the human striatum. Specifically, we demonstrated that one particular mode of functional connectivity in the striatum showed a high spatial correspondence to DaT availability, a marker of dopaminergic projections derived from DaT SPECT imaging. This observation generated multiple hypotheses that we validated using both data from PD patients and healthy controls. We showed that this secondorder striatal connectivity mode is associated with symptom severity and sensitive to acute dopaminergic modulation by L-DOPA in persons with PD, a disorder characterized by a degeneration of dopaminergic neurons projecting to striatum (Fearnley and Lees, 1991;Brooks and Piccini, 2006;Hornykiewicz, 2008). We also demonstrated that this mode is associated with the amount of tobacco and alcohol use, both of which have been related to alterations in dopaminergic signalling (Laruelle et al., 1995;Barrett et al., 2004;Nutt et al., 2015). As such, our results provide evidence that the second-order mode of functional connectivity in striatum maps onto dopaminergic projections and can be used as a non-invasive biomarker for investigating dopaminergic (dys)function in PD and substance use. While our results still need to be replicated out of sample to warrant immediate application in clinical practice, formal quantification of test-retest reliability already suggests that this gradient approach has very high measurement consistency and therefore lends itself for further investigation into the clinical utility across the various neurological and psychiatric disorders associated with dopaminergic functioning.
By applying connectopic mapping, we shift away from the vast majority of resting-state fMRI studies that employ hard parcellations to investigate functional brain connectivity. Gradient-based approaches such as connectopic mapping were only developed recently, but have already been successfully employed in several studies to investigate functional connectivity in cortical (Haak et al., 2018;Margulies et al., 2016;Saadon-Grosman et al., 2020) and subcortical regions (Marquand et al., 2017;Yang et al., 2020;Tian et al., 2020). Recent work (Hong et al., 2020) has furthermore demonstrated that connectivity gradients can generally be obtained with high reproducibility and reliability, and can predict phenotypic variations with higher accuracy than connectivity measures derived from traditional parcellation-based approaches, making them of interest for potential biomarker development. Indeed, hard parcellations only allow investigating the average functional connectivity signal in one or more regions of interest and thereby ignore both the topographic organization of and functional multiplicity in the brain. In contrast, gradient-based approaches do not only enable characterization of smooth, gradual changes in functional connectivity, but also enable the detection of multiple, overlapping modes of functional connectivity in a region that might exist at the same time (Haak and Beckmann, 2020). The work presented here signifies the importance of both features by not only showing that the second-order striatal connectivity mode comprises a smooth gradient from the dorsal putamen and dorsal caudate to the ventral putamen and ventral caudate including the NAcc, but also that this second-order mode -but not the zeroth-order or first-order mode -maps onto DaT availability. In doing so, we are the first to demonstrate a direct mapping between a functional connectivity-derived marker and dopaminergic projections.
DaT is highly expressed in the terminals of dopaminergic neurons projecting from the midbrain to striatum (Brooks, 2016). The high spatial correlation (r = 0.884) between the group-average secondorder connectivity mode in the HCP dataset and the group-average DaT SPECT image in the PPMI dataset as well as the significant within-subject spatial correlation between the connectivity mode and DaT SPECT scan (r = 0.58) in PPMI subjects therefore suggests that this connectivity mode maps onto these dopaminergic projections to striatum. We further demonstrated that the association of the second-order connectivity mode with the DaT SPECT scan was stronger than that of all the other investigated PET markers indexing various neurotransmitter systems (see Figure 2-figure supplement 1). This figure also shows that correlations of the second-order connectivity mode with dopamine receptors D1 and D2 in striatum, which are present on postsynaptic dopaminergic neurons, were substantially lower and not significant (r = −0.290, p=0.086 and r = 0.241, p=0.156), suggesting that the second-order connectivity mode is specific to presynaptic dopaminergic projections. Animal work has furthermore shown that dopaminergic projections form a gradient with nigrostriatal neurons from SNc projecting predominantly to the dorsolateral striatum (putamen and caudate) and mesolimbic neurons from the VTA projecting predominantly to the ventromedial striatum (NAcc; Steiner and Tseng, 2016;Haber, 2014;Björklund and Dunnett, 2007) representing a functional connectivity gradient formed by the SNc projections to the dorsolateral (putamen/caudate) and VTA projections to the ventromedial striatum (NAcc). Studies demonstrating that the average striatal DaT binding as obtained by DaT SPECT or PET imaging is highly correlated with averaged post-mortem SN cell counts in humans are in support of this view (Snow et al., 1993;Colloby et al., 2012;Kraemmer et al., 2014). However, to our knowledge the relationship between DaT SPECT/PET and VTA cell counts in humans has not been investigated, and future work will thus be necessary to determine the exact relationship between this striatal connectivity mode, DaT availability assessed by DaT SPECT, and dopaminergic projections.
While we were able to replicate the spatial correlation between the second-order connectivity mode and the DaT SPECT scan at the within-subject level, this spatial correlation (r = 0.58) is not as high as the spatial correlations observed at the group level (i.e., r = 0.721 and r = 0.714 for PPMI controls and PD patients respectively, and r = 0.884 between the DaT SPECT scan in PPMI controls and the connectivity mode in HCP participants). This is not surprising given the relatively low temporal resolution of the resting-state fMRI scan of the PPMI dataset (TR = 2400 ms, 260 volumes). While this resolution is sufficient for typical resting-state fMRI analyses at the group level, the precise delineation of the very fine-grained and overlapping connectivity modes using connectopic mapping at the single-subject level calls for high spatial and temporal resolution data (Haak et al., 2018). However, to our knowledge, there is currently no dataset (publicly) available that includes both a high-resolution resting-state fMRI scan and a DaT SPECT scan from the same participants. With respect to Figure 3, we further note the difference in the within-subject correlation for the putamen (r = 0.61/0.62) compared to caudate-NAcc region (r = 0.51/0.44). We tentatively speculate that this difference might relate to a stronger and more stable dopamine-related resting-state fMRI signal in putamen compared to caudate-NAcc resulting from more dopaminergic projections to putamen (Hörtnagl et al., 2020), and the putamen being larger in size and spatially further away from the ventricles and therefore less susceptible to motion-related artefact than the caudate-NAcc region.
Nevertheless, adding to its association with dopaminergic projections are the alterations of this connectivity mode observed in PD. This disorder is characterized by a loss of nigrostriatal dopaminergic neurons projecting from SNc to the striatum, which is most prominent in the putamen (Fearnley and Lees, 1991). Corresponding to the pathology of the condition, we observed a significant difference of this connectivity mode between right-dominant PD patients and control participants in bilateral putamen. Moreover, within this patient group the second-order striatal connectivity mode was also sensitive to inter-subject variability as revealed by the association with symptom severity. This association is visualized in Figure 4B, which shows that portions of the gradient in right putamen that map on low DaT availability (blue) increase as symptom severity in PD increases. We argue that the second-order striatal connectivity mode hereby follows the expected pattern of a reduction in dopaminergic projections to the putamen (as indexed by decreased DaT availability) as symptom severity increases in PD.
While group differences in the connectivity mode were thus present in bilateral putamen, the association with symptom severity was driven by the TSM coefficients modelling the right putamen. This latter finding might appear counterintuitive as a tremor dominant to the right side of the body in PD (i.e., right-dominant PD) is assumed to correspond with a dopamine depletion that is dominant to the contralateral striatum, that is, left striatum. However, post-mortem studies have reported so-called flooring effects by demonstrating a complete absence of dopaminergic fibres in the dorsal putamen in the most affected hemisphere in PD patients ≥ 4 years after disease onset (Kordower et al., 2013). Furthermore, resting-state fMRI studies have reported larger PD vs. control group differences in the anterior putamen of the lesser affected compared to the more affected hemisphere (Helmich et al., 2010). Taken together, these findings might suggest that in many patients with right-dominant PD in our sample (mean disease duration is 3.96 years) dopamine-related connectivity in the contralateral left putamen is showing flooring effects, making associations with symptoms only detectable in the right putamen. Future research is necessary to confirm this hypothesis. It should also be noted that we did not find a significant difference between controls and the left-dominant PD group. Since there is no evidence that different mechanisms underlie left-dominant and right-dominant PD (apart from the difference in the most-affected hemisphere), this might be a power issue that requires further investigation.
Not only was the second-order mode of connectivity in striatum sensitive to variability in symptom severity in a clinical cohort, but also to behavioural variability associated with self-reported alcohol and tobacco use in a healthy, non-clinical population. Substance use has frequently been associated with alterations in dopamine release in the ventromedial striatum (NAcc), part of the mesolimbic dopaminergic pathway (Laruelle et al., 1995;Barrett et al., 2004;Nutt et al., 2015). Corresponding with these findings, we observed significant associations with the amount of alcohol and tobacco use over the past week in the caudate-NAcc region, but not in the putamen region of the second-order striatal connectivity mode. More specifically, we observed subtle increases in blue-coded voxels and decreases in red-coded voxels as substance use increased, suggesting decreased DaT availability or more generally decreased dopaminergic signalling in the caudate-NAcc region in high nicotine and alcohol users. These results are consistent with findings from previous DaT SPECT and PET studies reporting reductions in striatal DaT availability in patients with alcohol dependence (Grover et al., 2020;Yen et al., 2016;Laine et al., 1999;Repo et al., 1999) and nicotine dependence (Yang et al., 2008). However, a limitation that should be mentioned is that the resting-state fMRI sequence of the HCP dataset has not optimized for subcortical brain regions.
As such, not only are these alterations in PD and high alcohol and tobacco users of the HCP dataset consistent with the hypothesis that the second-order striatal connectivity mode reflects dopaminergic projections, the alterations are also specific to the hypothesized striatal subregions and dopaminergic pathways. That is, we found that PD -a disorder characterized by death of nigrostriatal dopaminergic neurons leading to motor impairments-was associated with connectivity alterations in the putamen, which is a key region of the nigrostriatal pathway that has predominantly been implicated in motor function (Joshua et al., 2009;Faure et al., 2005). On the other hand, tobacco and alcohol use were associated with connectivity alterations in the caudate-NAcc region, which is part of the mesolimbic pathway that has repeatedly been implicated in reward processing and substance use (Schultz, 2013;Wise, 2004).
Finally, we observed that the change in the second-order mode of connectivity in striatum induced by L-DOPA administration was associated with the change in symptom severity in PD patients. L-DOPA is used as a drug for the treatment of PD, yet not all patients are equally responsive to L-DOPA treatment. When L-DOPA crosses the blood-brain barrier, it is converted into dopamine and is assumed to increase dopaminergic signalling (Lewitt, 2008). However, there are differences between PD patients in treatment response, which can be explained by a variety of factors, including differences in the level of systemic L-DOPA uptake from the gut (Nonnekes et al., 2016). Our finding thus indicates that the second-order striatal connectivity mode is differentially sensitive to acute dopaminergic modulation across PD patients and that the amount of this change is associated with the amount of change in symptom severity. This adds to our hypothesis that this mode is associated with dopamine-related functional connectivity and furthermore indicates that studying the dopaminergic system by applying connectopic mapping to resting-state fMRI offers advantages over PET and SPECT scans: PET and SPECT are not only invasive and limited by their low spatial resolution but also depend on indirect measures of dopaminergic signalling such as availability of DaTs and receptors and are therefore not very sensitive to acute, temporal alterations in dopaminergic signalling. In contrast, here we show that connectopic mapping does allow for the investigation of both fine-grained spatial and short-term temporal changes in dopamine-related functional connectivity.
In conclusion, our results provide evidence that the second-order mode of resting-state functional connectivity in striatum is associated with dopaminergic projections and can be developed into a noninvasive biomarker for investigating dopaminergic (dys)function. This may have wide-ranging clinical and scientific applications across disorders associated with dopaminergic functioning. For example, in the diagnostic work-up of movement disorders where DaT SPECT is currently used to distinguish between PD and essential tremor or dystonic tremor, the resting-state fMRI derived second-order connectivity mode might be used instead. The correlation with symptom severity suggests that this mode might also be used as a progression biomarker, for example, to track differences in rate of progression in future intervention studies of new experimental medications aimed at modifying the course of PD. Our results furthermore suggest that this striatal connectivity mode is associated with functions of both the nigrostriatal and mesolimbic pathway, and that it might be possible to differentiate between the two dopaminergic pathways by considering in which striatal subregion that gradient is altered: connectivity alterations seem to occur in putamen for functions associated with the nigrostriatal pathway and in ventral caudate/NAcc for functions associated with the mesolimbic pathway. However, the exact mapping of this striatal connectivity mode on both pathways as well as its relation with the first-order, ventromedial-to-dorsolateral striatal gradient, which we previously linked to goal-directed behaviours, is subject for further investigation.

Materials and methods
Resting-state fMRI data of the HCP dataset For our first analysis, we used resting-state fMRI data from the HCP, an exceptionally high-quality, publicly available neuroimaging dataset . HCP participants were scanned on a customized 3 T Siemens Skyra scanner (Siemens AG, Erlanger, Germany) and underwent two sessions of two 14.4 min multiband accelerated (TR = 0.72 s) resting-state fMRI scans with an isotropic spatial resolution of 2 mm. Here, we included participants from the S1200 release who completed at least one resting-state fMRI session (2 × 14.4 min) and for whom data was reconstructed with the r227 reconstruction algorithm. (The reconstruction algorithm was upgraded in late April 2013 from the original 177 ICE version to the 227 upgraded ICE version. As the reconstruction version has been shown to make a notable signature on the data that can make a large difference in fMRI data analysis [for details, see https://wiki.humanconnectome.org/display/PublicData/Ramifications+of+Image+ Reconstruction+Version+Differences], we only included participants with r227 reconstructions.) This resulted in the inclusion of 839 participants (aged 22-37 years; 458 females). Resting-state fMRI data were preprocessed according to the HCP minimal processing pipeline (Glasser et al., 2013), which included corrections for spatial distortions and head motion, registration to the T1w structural image, resampling to 2 mm MNI152 space, global intensity normalization, and high-pass filtering with a cutoff at 2000s. The data were subsequently denoised using ICA-FIX -an advanced independent component analysis-based artefact removal procedure (Salimi-Khorshidi et al., 2014) -and smoothed with a 6 mm kernel.

Connectopic mapping of the striatum in the HCP dataset
We estimated connection topographies from the HCP resting-state fMRI data using the first session (2 × 14.4 min) for each subject. To this end, we used connectopic mapping (Haak et al., 2018), a novel method that enables the dominant modes of functional connectivity change within the striatum to be traced on the basis of the connectivity between each striatal voxel and the rest of the brain (see Figure 1). In previous work, we showed that the dominant mode (zerothorder mode) of connectivity in the striatum obtained with connectopic mapping represented its anatomical subdivision into putamen, caudate, and NAcc. Since higher-order modes are restricted by lower-order modes, we decided to take the anatomical subdivision in the striatum into account by applying connectopic mapping in the current work to the left and right putamen and caudate-NAcc striatal subregions separately, thereby also increasing regional specificity. When referring to the second-order mode of connectivity in striatum, we thus refer to the combination of the second-order connectivity modes of putamen and caudate-NAcc. We did not apply connectopic mapping to the NAcc and caudate separately as the left NAcc and right NAcc only include 136 voxels and 127 voxels, respectively. We expect that this very small region is too homogenous in terms of connectivity with cortex to estimate reliable overlapping connectivity modes. Masks for the striatal regions were obtained by thresholding the respective regions from the Harvard-Oxford atlas at 25% probability.
In brief, we rearranged the fMRI time-series data from each striatal subregion and all grey-matter voxels outside the striatum into two time-by-voxels matrices. Since the latter is relatively large, we reduced its dimensionality using a lossless singular value decomposition (SVD). We then computed the correlation between the voxel-wise striatal time-series data and the SVD-transformed data from outside the striatum, and subsequently used the η 2 coefficient to quantify the similarities among the voxel-wise fingerprints (Haak et al., 2018). Next, we applied the Laplacian eigenmaps non-linear manifold learning algorithm (Belkin and Niyogi, 2002) to the acquired similarity matrix, which resulted in a set of overlapping, but independent, vectors representing the dominant modes of functional connectivity change across striatum (i.e., connection topographies). Note that this can be done at the group level by using the average of the individual similarity matrices or individually for each subject (as used for statistical analysis). For each subject, modes were aligned to the group-level connectivity mode (by inversion if negatively correlated) to enable visual and statistical comparisons across subjects. We selected the second-order striatal connectivity mode (both the group-average and subject-specific modes) for further analyses.
Finally, to enable statistical analysis over these connection topographies, we fitted spatial statistical models to obtain a small number of coefficients summarizing the second-order connectivity mode of each striatal subregion in the X, Y, and Z axes of MNI152 coordinate space. For this, we use 'trend surface modelling' (Gelfand et al., 2010), an approach originally developed in the field of geostatistics, but that has wide-ranging applications due to its ability to model the overall distribution of properties throughout space as a simplified surface. Here, we use the TSM approach to predict each individual subject's connection topography by fitting a set of polynomial basis functions defined by the coordinates of each striatal location. We fit these models using Bayesian linear regression (Bishop, 2006), where we employed an empirical Bayes approach to set model hyperparameters. Full details are provided elsewhere (Bishop, 2006), but this essentially consists of finding the model hyperparameters (controlling the noise-and the data variance) by maximizing the model evidence or marginal likelihood. This was achieved using conjugate gradient optimization. For fixed hyperparameters, the posterior distribution over the trend coefficients can be computed in closed form. This, in turn, enables predictions for unseen data points to be computed. We used the maximum a posteriori estimate of the weight distribution as an indication of the importance of each trend coefficient in further analyses. To select the degree of the interpolating polynomial basis set, we fit these models across polynomials of degree 2-5 and then compared the different model orders using a scree plot analysis (Cattell, 1966). This criterion strongly favoured a polynomial of degree 2 for the putamen subregion and a polynomial of degree 4 for the caudate-NAcc subregion. This means that the connectivity mode in putamen was modelled with linear and quadratic functions in the X, Y, and Z directions of MNI152 coordinate space (six TSM coefficients) and the connectivity mode in the caudate-NAcc region with linear, quadratic, cubic, and quartic functions in the X, Y, and Z directions of MNI152 coordinate space (12 TSM coefficients). The TSM coefficients of the fitted polynomial basis functions describe the rate at which the connectivity mode changes along a given spatial dimension and can be used for statistical analysis. The polynomials summarized the connectivity modes well, explaining the following mean ± SD of the variance: left putamen: 90.5% ± 4.16%; right putamen: 90.2% ± 4.64%; left caudate-NAcc: 88.6% ± 2.54%; right caudate-NAcc: 89.4% ± 2.15%.

DaT SPECT imaging in the PPMI dataset
Mapping the second-order striatal connectivity mode onto DaT availability Next, we quantified the similarity between the second-order mode of connectivity in striatum and the DaT SPECT image. To this end, we combined the average (i.e., group level) second-order connectivity modes of putamen and caudate-NAcc obtained in the high-resolution HCP dataset and computed the voxel-wise spatial correlation of this mode with the average (i.e., group level) DaT SPECT image of striatum obtained in the PPMI dataset. Given that the voxel-wise spatial correlation between both images might be inflated due to potential spatial autocorrelation effects (i.e., the images represent smooth spatial functions), we additionally computed the correlation between the TSM coefficients modelling the group-average connectivity mode and the group-average DaT SPECT scan since TSM coefficients are orthogonal. More specifically, the TSM coefficients modelling the left and right putamen (2 × 6) and caudate-NAcc regions (2 × 12) were combined and the correlation was computed across all these 36 TSM coefficients. In addition, to show that our results do not heavily depend on the chosen model order, we repeated this analysis using model order 3 (i.e., a cubic model with nine TSM coefficients) for both the putamen and caudate-NAcc regions (i.e., 4 × 9 = 36 TSM coefficients).
For a subsample of PPMI participants with a DaT SPECT scan, there was also a low-resolution resting-state fMRI scan available (130 datasets from PD patients and 14 from controls). We therefore also investigated the within-subject spatial correspondence between the DaT SPECT scan and the second-order connectivity mode for these subjects in the PPMI dataset. This procedure is detailed in Appendix 1-Supplementary analyses.

Resting-state fMRI data of the PD dataset
Given that PD is characterized by a loss of dopaminergic neurons (Fearnley and Lees, 1991), we investigated whether the second-order striatal connectivity mode was altered in PD. For this analysis, we used high-resolution resting-state fMRI data from a cohort consisting of 39 patients with PD (aged 38-81, 23 females) and 20 controls (aged 42-80, 9 females), recruited at the Centre of Expertise for Parkinson & Movement Disorders at the Radboud University Medical Center (Radboudumc) in Nijmegen and scanned at the Donders Institute in Nijmegen, the Netherlands (Dirkx et al., 2019). All patients were diagnosed with idiopathic PD (according to the UK Brain Bank criteria), and all patients had a mild to severe resting tremor besides bradykinesia. In 20 patients, the motor symptoms were right-dominant, in 19 patients left-dominant (dominance here refers to the side of the body displaying the most prominent motor symptoms [including tremor]), which is believed to correspond with a dopamine depletion dominant to contralateral hemisphere in the brain. The study was approved by the local ethics committee, and written informed consent according to the Declaration of Helsinki was obtained from all participants. Detailed sample characteristics can be found in Table 2.
Patients with PD underwent two 10 min resting-state fMRI sessions, that is, a placebo session and L-DOPA session, separated by at least a day on a 3 T Siemens Magnetom Prisma fit scanner. Resting-state fMRI scans were obtained with an interleaved high-resolution multiband sequence (TR = 0.860 s, voxel size = 2.2 mm isotropic, TE = 34 ms, flip angle = 20°, 44 axial slices, multiband acceleration factor = 4, volumes = 700). Under both conditions, patients were scanned after overnight fasting in a practically defined off state, that is, more than 12 hr after intake of their last dose of dopaminergic medication. During one session, patients were scanned after administration of L-DOPA, that is, they received a standardized dose of 200/50 mg dispersible levodopa/benserazide. During the other session, patients received placebo (cellulose powder). The cellulose powder and L-DOPA/benserazide were dissolved in water and therefore undistinguishable (visually and olfactory) for the participants as confirmed by a pilot study. Patients also received 10 mg domperidone to improve gastrointestinal absorption of levodopa and reduce side effects. The order of sessions was counterbalanced and the resting-state fMRI scan started on average 48 min (range: 25-70 min) after taking L-DOPA or placebo. Symptom severity was assessed during both sessions with part III (assessment of motor function by a clinician) of the Movement Disorder Society UPDRS (Goetz et al., 2008), and an electromyogram (EMG) of the hand was recorded to monitor tremorrelated activity. In light of ethical considerations, control participants did not receive L-DOPA and placebo, they just underwent two typical resting-state fMRI sessions during which the UPDRS was not administered.
Preprocessing of the resting-state fMRI data included removal of the first five volumes to allow for signal equilibration, primary head motion correction via realignment to the middle volume MCFLIRT (Jenkinson et al., 2002), grand mean scaling, and spatial smoothing with a 6 mm FWHM Gaussian kernel. The preprocessing pipeline was furthermore designed to rigorously correct for potential tremor-induced head motion-related artefacts. To this end, we used ICA-AROMA (Pruim et al., 2015), an advanced ICA-based motion correction procedure to identify and remove secondary head motionrelated artefacts with high accuracy while preserving signal of interest (Pruim et al., 2015;Parkes et al., 2018). Next, any remaining motion artefacts were removed from the data by regressing out the EMG parameters in addition to the white matter and CSF signal (Helmich et al., 2012). Finally, the data were temporally filtered with a high-pass filter of 0.01 Hz before being resampled to 2 mm MNI152 space.
We conducted four different analyses. All these analyses were conducted separately for the leftand right-dominant PD groups -given that the dopamine depletion is dominant to different hemispheres in these two groups -and separately for the putamen and caudate-NAcc subregions (left + right hemisphere combined) to increase regional specificity as PD is known to affect the putamen region of the striatum before the caudate-NAcc region. In all these analyses, we therefore corrected for multiple comparisons using a Bonferroni-corrected α-level of 0.0125 (0.05/4 (2 patient groups * 2 subregions)).
In our first analysis, we compared the second-order striatal connectivity mode between the control and PD groups. To this end, we conducted omnibus tests comparing the TSM coefficients modelling the second-order striatal gradient of the placebo session in the PD group with the TSM coefficients modelling the striatal gradient of the first session of the control group. More specifically, group differences in the TSM coefficients were assessed by using a likelihood ratio test in the context of a logistic regression. We report the X 2 (likelihood test) and corresponding p-value of tests that revealed significant group differences. Second, since PD is a heterogeneous disease, we also conducted an analysis taking this variability into account by investigating associations between the TSM coefficients modelling the second-order striatal gradient and symptom severity in PD patients. To this end, we conducted GLMs that included the TSM coefficients modelling the gradient during the placebo session to predict UPDRS symptom severity scores. For all identified associations with UPDRS symptom severity, we report the X 2 (likelihood test) and the corresponding p-values, and post-hoc compute Pearson correlations between UPDRS symptom severity and the individual TSM coefficients to determine which coefficients most strongly contributed to the effect.
Third, we assessed differences in the second-order striatal connectivity mode between the placebo and L-DOPA session in both PD groups. More specifically, session differences in the TSM coefficients were assessed by using a likelihood ratio test in the context of a logistic regression. We report the X 2 (likelihood test) and corresponding p-value of tests that revealed significant differences between the placebo and L-DOPA session. Finally, treatment response to L-DOPA is known to differ among patients with PD. To take this variability across patients into account, we also investigated whether the L-DOPA-induced change in the second-order striatal connectivity mode was associated with L-DOPA-induced changes in UPDRS symptom severity. More specifically, we calculated the difference in UPDRS symptom scores and TSM coefficients between the placebo and L-DOPA session and investigated associations within the GLM framework. For all identified associations, we post-hoc computed Pearson correlations between the change in UPDRS symptom severity and the change in individual TSM coefficients to determine which coefficients most strongly contributed to the effect.
Investigating the second-order striatal mode in relation to tobacco and alcohol use Given that alterations in dopaminergic functioning have also been implicated in substance use, we also investigated the association between the second-order striatal connectivity mode and tobacco use as well as alcohol use across high users within the HCP dataset. To this end, we selected HCP participants testing negative for acute drug and alcohol use but who reported to have consumed ≥3 light and/ or ≥1 heavy alcoholic drinks per day over the past week (N = 30), and participants reporting to have smoked ≥5 cigarettes every day over the past week (N = 38). Effects of smoking and drinking were analysed separately, and we again modelled the second-order striatal connectivity mode separately for left and right putamen and left and right caudate-NAcc. We conducted GLMs that included next to the amount of use over the past week (i.e., the total number of alcohol drinks or the total number of times tobacco was smoked), the TSM coefficients modelling the second-order striatal connectivity mode (i.e., the TSM coefficients modelling the putamen or caudate-NAcc). Multiple comparison correction was applied using a Bonferroni-corrected α-level of 0.0125 (2 substances * 2 striatal subregions). For all identified associations with the amount of alcohol use or tobacco use, we report the X 2 (likelihood test) and the corresponding p-values and post-hoc compute Pearson correlations between the amount of use and the individual TSM coefficients to determine which coefficients most strongly contributed to the association.

Supplementary analyses
To further demonstrate the high specificity of the second-order connectivity mode to the DaT SPECT scan over and above other PET scans (i.e., other neurotransmitter systems), we computed correlations with the TSM coefficients of all PET scans, tapping into various neurotransmitter systems, included in the publicly available JuSpace toolbox (Dukart et al., 2021). This analysis is described in Appendix 1-Supplementary analyses. We also investigated whether the mapping of the second-order connectivity mode onto DaT SPECT scan was influenced by residual head motion. Furthermore, for the other analyses described above (effects of diagnosis and L-DOPA in the PD dataset and associations with tobacco and alcohol use in the HCP dataset), we conducted post-hoc sensitivity analyses to rule out that the group differences and associations revealed by our analyses were dependent on age and sex. In addition, we investigated whether the associations with the amount of tobacco and alcohol use persisted under different usage thresholds. All these analyses are also described in Appendix 1-Supplementary analyses. We investigated the spatial correspondence of the second-order connectivity mode to multiple PET markers indexing various neuromodulatory systems. To this end, we made use of various PET scans tapping into different neurotransmitter systems (group-averages of 11-36 controls) implemented in the publicly available JuSpace toolbox (https://github.com/juryxy/JuSpace). The included PET scans for the serotonin system are 5HT1a receptor (5HT1a_WAY), 5HT1b receptor (5HT1b_P943), 5HT2a receptor (5HT2a_ALT), and the serotonin transporter (SERT_DASB and SERT_MADAM). The included PET scans for the dopamine system (other than DaT SPECT) are dopamine type 1 receptor (D1_SCH23390), dopamine type 2 receptor (D2_RACLOPRIDE), and FDOPA (FDOPA_f18). In addition, the following neurotransmitter systems were also included: for GABA the GABAa receptor (GABAa_FLUMAZENIL), for noradrenalin the noradrenalin receptor (NAT_MRB), and for the opiod system the mu opiod receptor (MU_CARFENTANIL). We applied TSM to each of these PET scans in the striatum (as in the main analysis, the rest of the brain was masked and not included) and computed correlations between TSM coefficients modelling the second-order connectivity mode and the TSM coefficients modelling these PETderived markers. For each PET scan, the correlations with the TSM coefficients were normalized using the Fisher r-to-z transformation and the absolute correlation was taken. These normalized correlations are visualized in the top panel of Figure 2-figure supplement 1. As can be observed, the correlation of the second-order connectivity mode with the DaT SPECT scan is substantially higher than that of any other PET marker. Though it is noticeable that some of the correlations with the PET markers for the serotonin system (SERT_DASB transporter and 5HT1b receptor) that are also known to have a high density in the striatum are also relatively high. Nevertheless, these markers only reach about half of the correlation value of DaT SPECT. To support the robustness and significance of the correlation of our second-order connectivity mode with the DaT SPECT scan over and above the correlation with the markers of the serotonergic system, we tested the correlation between the TSM coefficients obtained for the second-order connectivity mode and the DaT SPECT scan in striatum for significance using permutation testing (N = 10,000). More specifically, we permuted corresponding TSM coefficients obtained for each of the PET markers and thereby generated a null distribution by computing the absolute (Fisher r-to-z normalized) correlations between the connectivity mode TSM coefficients and the permuted TSM coefficients of the other PET markers. Permutations were conducted separately for each coefficient, not permuting across coefficients to ensure interchangeability under the null assumption of no differentiation across different PET markers.

Competing interests
As can be observed in the bottom panel of Figure 2-figure supplement 1, all permuted correlations are lower than the correlation observed between the DaT SPECT scan and the connectivity mode, indicating that the observed correlation between the connectivity mode and DaT SPECT scan is highly significant and unlikely to be obtained by chance. Furthermore, using this null distribution, we defined the Bonferroni-corrected threshold for significance corresponding to p=0.0008 (i.e., p=0.01/12 PET and SPECT scans), which we added to the top figure displaying the correlations with the other PET tracers. This not only confirms that our results are highly significant, but also that the correlations obtained for the other PET markers, including those of the serotonin system, are not only substantially lower, but also do not pass the threshold for significance based on the null distribution. Of note, the correlation with other markers of the dopamine system, such as the D1 and D2 receptor, as opposed to DaT SPECT is not particularly high. However, this is not surprising since these receptors are present on postsynaptic neurons and are likely representative of postsynaptic dopaminergic projections from striatum to cortex, rather than the presynaptic dopaminergic projections from the midbrain to the striatum reflected by the DaT SPECT scan.
Within-subject correspondence between the second-order striatal connectivity mode and the DaT SPECT scan In the article, we demonstrated that the second-order striatal connectivity mode at the group level (obtained by averaging this mode across all 839 HCP subjects) showed a very high spatial correlation (r = 0.884) with the group-level DaT SPECT image of striatum (obtained by averaging the DaT SPECT images across all 209 PPMI controls). We also aimed to demonstrate that this mapping can be replicated at the within-subject level by investigating the within-subject spatial correspondence between this connectivity mode and the DaT SPECT scan acquired in the PPMI dataset. However, while the PPMI dataset has resting-state fMRI data available for a small subsample of its participants (14 controls with one resting-state fMRI dataset each and 82 PD patients with 130 resting-state fMRI datasets combined; in case of multiple assessments per subject they were separated by at least 1 year), it is of a relatively low temporal and spatial resolution (TR = 2400 ms, 210 time points, 3.3 mm isotropic resolution) compared to the HCP data (TR = 720 ms, 2,400 time points, 2.0 mm isotropic resolution). While this resolution is sufficient for typical resting-state fMRI analyses, the precise delineation of the very fine-grained and overlapping connectivity modes using connectopic mapping calls for high-resolution data. The single-subject connectivity modes in the PPMI dataset (as opposed to the HCP single-subject modes and grouplevel modes) might therefore not be of sufficient quality and reliable for every subject. To address this issue, we first computed the spatial correlation of each subject's individual connectivity mode with that of the group-average HCP connectivity mode as well as with the DaT SPECT scan of each subject (see Figure 3-figure supplement 1). In this analysis, the second-order striatal connectivity mode was modelled separately (and correlations were calculated separately) for the left and right putamen and caudate-NAcc subregions. This revealed highly significant positive correlations (0.68 > r < 0.91, all p<4.0e21) across both controls and patients, suggesting that if the connectivity mode of a subject resembles the HCP group-average connectivity mode -assumed to be an index of good quality -a high spatial similarity can be observed between the connectivity mode and the DaT SPECT scan of that subject. Next, we selected those subjects with good-quality connectivity modes as determined by a spatial correlation of r > 0.5 with the group-average connectivity mode in the HCP dataset. Within this sample of 73-86 datasets from PD patients and 6-8 datasets from controls (dependent on the striatal subregion), we not only replicated the spatial correspondence between the connectivity mode and DaT SPECT scan at the group level (patients: r = 0.714; control group: r = 0.721) but also observed significant within-subject spatial correlations (0.44 > r < 0.63; mean = 0.58, 95% CI = [0.56,0.60]) between the connectivity mode and DaT SPECT scan (see Figure 3). While we were able to replicate the spatial correlation between the second connectivity mode and the DaT SPECT scan at the within-subject level, this correlation (r = 0.58) is not as high as the spatial correlations observed in the group level (i.e., r = 0.721 and r = 0.714 for PPMI controls and PD patients, respectively, and r = 0.884 between the DaT SPECT scan in PPMI controls and the connectivity mode in HCP participants). This is, however, not surprising given the relatively low temporal and spatial resolution of the resting-state fMRI scan of the PPMI dataset. However, to our knowledge, there is currently no dataset available that includes both a highresolution resting-state fMRI scan and a DaT SPECT scan from the same participants.

Post-hoc analyses of head motion
To demonstrate that the mapping of the group-average second-order connectivity mode onto the group-average DaT SPECT scan (r = 0.884) was not influenced by residual head motion, we generated connectivity modes for the 10% highest movers (N = 84, meanFD range: 0.0376-0.0538) and 10% lowest movers (N = 84, meanFD range: 0.1354-0.3155) of the HCP dataset and computed the spatial correlation to the group-average DaT SPECT scan (N = 209, no head motion metrics were available, so we used the entire sample). This analysis revealed very similar connectivity modes (see Figure 3-figure supplement 2) and a very similar spatial correlation for the low FD group (r = 0.883) and the high FD group (r = 0.886), indicating that this mapping was not induced by residual head motion.
Next, we aimed to demonstrate that the mapping between the connectivity mode and the DaT SPECT scan at the within-subject level was also not influenced by residual head motion. To this end, we divided the subsample of the PPMI dataset where both a resting-state fMRI scan and DaT SPECT are available (used for mapping the second-order connectivity mode onto the DaT SPECT at the within-subject level) in half based on the meanFD (meanFD cutoff = 0.126) and computed the within-subject correlation between the DaT SPECT and connectivity mode for the 'low motion' and 'high motion' halves of the sample separately. This analysis revealed that the withinsubject spatial correlation between the DaT SPECT and second-order connectivity mode was virtually identical between the low and high meanFD samples (see Figure 3-figure supplement  2). Although it does appear that for the caudate-NAcc region the correlations between the connectivity mode and the DaT SPECT scan are slightly lower for the high motion half than for the low motion half. This indicates that the high spatial correlation between the connectivity mode with the DaT SPECT scan is not artificially induced by head motion, but that residual head motion instead weakens this correlation. All these additional analyses thus suggest that the proposed biomarker is not picking up on residual motion.

Post-hoc analyses of age and sex
For all the analyses described in the article (effects of diagnosis and L-DOPA in the PD dataset and associations with smoking and drinking in the HCP dataset), we conducted post-hoc sensitivity analyses to rule out that the group differences and behavioural associations revealed by our analyses were dependent on age and sex. To this end, we conducted two types of analyses. First, we repeated our main analyses by including covariates for age and sex in our statistical models in addition to the TSM coefficients to verify that effects remained (close to) significant when including these demographic variables. Next, we only included age and sex in our statistical models (without the TSM coefficients) to verify that effects were not explained by age and/or sex only. The outcomes of these analyses (X 2 and p-value) are listed in Appendix 2-table 1 and demonstrate that none of the significant effects observed in our main analyses were dependent on age or sex. However, adding age and sex (age in particular) did increase the significance of findings substantially for the analyses investigating the L-DOPA-induced changes. This might be explained by the fact that patients who are older often have more severe PD and do not benefit as much anymore from L-DOPA treatment.

Post-hoc analyses using different usage thresholds for tobacco and alcohol use
We also investigated whether the associations of the second-order mode of connectivity in striatum with the amount of tobacco use and alcohol use persisted under different usage thresholds. For both tobacco and alcohol use, we chose a daily usage threshold lower (≥2× tobacco/≥1× alcoholic drink) and a daily usage threshold higher (≥8× tobacco/≥3× alcoholic drink) than the one used in the main analysis (≥5× tobacco/≥3× light alcoholic and/or ≥1× hard liquor drinks a day). Please note that the aim of these analyses is not necessarily to show that effects remain significant as under different usage thresholds the sample size and statistical power will change, but rather that the explained variance remains high. Nevertheless, apart from the low-usage threshold for alcohol use, all effects also remained significant, as can be observed in Appendix 2-table 2 and Appendix 2-table 3, indicating that the associations with tobacco and alcohol use were not specific to the chosen usage threshold. However, a pattern that is visible is that associations become stronger when only including the highest users in this population-based sample in the analysis.