Mesolimbic Dopamine Function Is Related to Salience Network Connectivity: An Integrative Positron Emission Tomography and Magnetic Resonance Study

Background A wide range of neuropsychiatric disorders, from schizophrenia to drug addiction, involve abnormalities in both the mesolimbic dopamine system and the cortical salience network. Both systems play a key role in the detection of behaviorally relevant environmental stimuli. Although anatomical overlap exists, the functional relationship between these systems remains unknown. Preclinical research has suggested that the firing of mesolimbic dopamine neurons may activate nodes of the salience network, but in vivo human research is required given the species-specific nature of this network. Methods We employed positron emission tomography to measure both dopamine release capacity (using the D2/3 receptor ligand 11C-PHNO, n = 23) and dopamine synthesis capacity (using 18F-DOPA, n = 21) within the ventral striatum. Resting-state functional magnetic resonance imaging was also undertaken in the same individuals to investigate salience network functional connectivity. A graph theoretical approach was used to characterize the relationship between dopamine measures and network connectivity. Results Dopamine synthesis capacity was associated with greater salience network connectivity, and this relationship was particularly apparent for brain regions that act as information-processing hubs. In contrast, dopamine release capacity was associated with weaker salience network connectivity. There was no relationship between dopamine measures and visual and sensorimotor networks, indicating specificity of the findings. Conclusions Our findings demonstrate a close relationship between the salience network and mesolimbic dopamine system, and they are relevant to neuropsychiatric illnesses in which aberrant functioning of both systems has been observed.


Supplemental Methods Participants
Ethical permission was obtained from the local ethics committee, and all participants provided informed written consent. Healthy controls were recruited via advertisements online and in local media. Subjects had no history of psychiatric or neurological disorders, and had a urine drug screen and pregnancy test (where appropriate) prior to scanning.

PET Data Acquisition and Analysis
Participants were not permitted to smoke or consume caffeine for four hours preceding the scan. After acquiring a CT scan for attenuation correction, PET images were acquired using a Siemens Biograph HiRez XVI PET scanner (Siemens Healthcare, Erlangen, Germany) at Imanova Centre for Imaging Sciences.

Experiment 1: FDOPA Study
One hour prior to scanning, participants received 400mg entacapone and 150mg carbidopa, to prevent formation of radiolabelled metabolites and reduce peripheral metabolism.
Approximately 160 MBq of 18 F-DOPA was administered by bolus intravenous injection. The quantification pipeline was consistent with previous works. 1 Correction for head movement during the scan was performed by denoising the non-attenuation-corrected dynamic images using a level 2, order 64 Battle-Lemarie wavelet filter. Frames were realigned to a single reference frame, acquired 20 minutes post-injection, employing a mutual information algorithm. 2,3 The transformation parameters were then applied to the corresponding attenuated-corrected dynamic images, creating a movement-corrected dynamic image, which was used in the analysis. Realigned frames were then summated to create an individual motion-corrected reference map for the brain tissue segmentation. The cerebellum was used as a reference region, and Ki cer was calculated with the Patlak-Gjedde graphical approach adapted for reference tissue input function 4 . Image processing and quantification was done using in-house code with MATLAB 2012b.

Experiment 2: PHNO Study
Approximately 170 MBq of 11 C-(+)-PHNO was administered by bolus injection. After the administration of the radiotracer, dynamic emission data were acquired continuously for 90 minutes. The dynamic images were reconstructed using a filtered back-projection algorithm into 31 frames (8 x 15 seconds, 3 x 60 seconds, 5 x 120 seconds, 15 x 300 seconds) with a 128 matrix, a zoom of 2.6 and a transaxial Gaussian filter of 5mm.
An individual parcellation of the brain was implemented in MIAKAT release 4.2.6 (http://www.miakat.org), 5 SPM12 and FSL (version 5.0.9). Cerebellar grey matter was used as the reference region, and the simplified reference tissue model (SRTM) was used to derive BPND from the regional time activity curves. 6,7 The magnitude of dexamphetamine-induced dopamine release within the limbic striatum was quantified as the percentage change in BPND in the dexamphetamine condition vs. baseline (no dexamphetamine) condition.

MRI Data Acquisition
Participants were instructed to remain still, keep awake, and keep their eyes closed. A structural image was obtained using a gradient-echo scan (TR=7.0s, TE=2.8s, flip angle=11°, in plane resolution=1mm x 1mm, slice thickness=1.2mm, 196 slices).

Experiment 2: PHNO Study
MRI data was obtained using a Siemens MAGNETOM Verio 3-T magnetic resonance imaging scanner. Functional imaging involved a multiband sequence based on the multiband EPI WIP v012b provided by the University of Minnesota, [8][9][10][11] using a multiband acceleration factor of 2. 238 volumes were acquired, consisting of 72 interleaved slices (2mm thickness, and inplane resolution of 3 x 3 mm), with a TR of 2000 ms, TE of 30 ms and a scan time of 7 minutes 56 seconds.

fMRI Preprocessing
Image pre-processing was performed via the CONN toolbox (version 17.b) 13 for Statistical Parametric Mapping software (SPM 12 (6906)). A standard preprocessing pipeline was used consisting of slice timing correction, realignment, and normalisation to MNI space. Images were smoothed with a Gaussian kernel of 8mm full-width-half-maximum. The ART toolbox was used to account for motion and artefact detection using anatomical component based correction (aCompCor) of temporal confounds relating to head movement and physiological noise. This method models noise effects at a voxel level based on estimates derived from principal components of noise regions of interest (white matter and CSF, eroded by one voxel to minimise partial volume effects), and then removes these from the BOLD timeseries using linear regression. Six residual head motion parameters and their first order temporal derivatives were also entered as regressors into the first level model. A confounding effect accounting for magnetisation stabilisation, and its first order derivative was entered.
Artifact/outlier scans (average intensity deviated more than 5 standard deviations from the mean intensity in the session, or composite head movement exceeded 0.9 mm from the previous image) were also regressed out. Preprocessed data were temporally bandpass filtered (0.008-0.09 Hz) Time-series were extracted from N=333 predefined nodes of interests of the Gordon cortical atlas. The salience and default mode network nodes of the Gordon atlas are displayed in Figure S1. For each participant, a graph representing a functional connectivity network was constructed, each edge representing the level of functional connectivity between a pair of nodes, which was computed as the z-transformed Pearson's correlation coefficient between their mean time-series.

MRI-analysis: Atlas Selection
The Gordon parcellation is based upon resting state boundary maps observed in a sample of 120 healthy young adults, and shows superior within parcel homogeneity when compared to other parcellations, making it an ideal choice for the analysis of resting state data. 14 In order to demonstrate robustness of our findings, we also undertook all analyses using two alternative atlases -the Power atlas (a collection of 264 10 mm diameter spheres derived from connectivity data in over 300 healthy volunteers performing various tasks), 15 and the CONN network atlas (a 32 node atlas in which nodes are defined on the basis of an independent components analysis of 497 subjects from the Human Connectome Project). 13

MRI-analysis: Network Strength
For a network formed of N nodes, the average network strength ̅ can be computed similarly to the link density ρ of unweighted networks 16

MRI-analysis: Community Detection
On the basis of the original Gordon atlas labels 41 nodes were a priori defined as belonging to the default mode network, and 44 to the cinguloopercular/salience network (referred to in the current paper as the salience network). 14 The Power atlas assigns 32 and 58, while the CONN atlas assigns 7 and 4 nodes to the salience and default mode networks respectively. As a result, the networks of interest were defined across a wide range of scales both in terms of node volume, and maximal network size.
In addition to the apriori network labels, however, we also ran a whole brain community detection algorithm for each atlas, 17 to generate definitions of the salience and default mode networks based on the connectivity patterns present in the current datasets. To accomplish this each subject's fully weighted functional connectome was subjected to the Louvain community detection algorithm, and the results of this were used to generate community assignments at the group level (individual level community assignments are not appropriate for subsequent analyses). 17 Due to the non-deterministic nature of the Louvain algorithm, a previously described consensus clustering approach was employed, 18 negative weights were treated symmetrically, and the gamma parameter was set to 1.7 as this produced community sizes in relative agreement with existing parcellation schemes.

Identifying Dopamine Associated Nodes-Network Based Statistic
In order to identify whether specific subnetworks show a significant relationship with limbic dopamine synthesis capacity we used the Network-Based Statistic (NBS) to investigate salience, default mode, sensorimotor, and visual networks separately (the method is summarised in Figure 2A in the main text). 19  We next averaged across individuals to create a group level graph. Proportional thresholding was performed on this group averaged matrix by assigning a value of 1 to all edges with connection strength above a set threshold, and setting all remaining edges to 0. We used 100 thresholds, retaining 20% of edges at the most lenient threshold, and 7% at the most stringent. There is no 'correct' set of thresholds but at more lenient thresholds one risks including a high degree of spurious connections, while at more stringent thresholds the graph became overly fragmented. The fact that this is a more lenient range than reported elsewhere is appropriate given we are investigating intranetwork connectivity, where there will be a lower proportion of spurious edges. 20 Graph metrics were computed using the Brain Connectivity Toolbox. 21 Node degree refers to the number of neighbours a node has, and is thus a measure of the local, direct importance of a node: ∑ , . 22 While this intuitively captures the relative importance of a node within a network, in correlation based graphs, it may also reflect membership of a larger community, as opposed to the importance of the node in information processing. 23 We therefore also calculated for each node betweenness centrality , 22 and the nodeparticipation coefficient . 24 The node betweenness centrality measures the proportion of shortest paths between all pairs of nodes that pass through it, and reflects its position as a potential information broker in the network. It is formally defined as: , with , the number of shortest paths between nodes m and n. 22 The node participation coefficient was calculated after first assigning each node to a community using the Louvain community detection algorithm. 17 A participation coefficient of zero means that all the edges of a node are restricted to its own community, indicating a rather local role, whereas a value approaching 1 means that its edges are evenly distributed among all the communities of the graph -indicating that the node plays a role in integrating different clusters of the graph. The participation coefficient is defined as 1 ∑ , with the number of communities and , the degree of a node restricted to community C.
At each MRI threshold, every node was ranked on each of these metrics, and the mean rank of each node across MRI thresholds was then calculated. We then set a rank threshold, and for each metric selected only the nodes ranking above it. If any node ranked above this threshold for all three metrics it was termed a combination hub (main text Figure 2B steps A-C), highlighting its importance as an all-round information processing node. By next lowering the rank threshold we gradually increased the number of nodes meeting combination hub criteria, and so defined sets of combination hubs comprising between 10 and 40% of the total number of nodes. In some cases it is possible a specific hub threshold might have no eligible nodes (e.g. in the main paper figure 4B -the 18 F-DOPA salience network does not have combination until the threshold reaches 15%).

Identifying Overlap Between Dopamine Associated Nodes and Network Hubs
After identifying nodes within the default mode and salience networks that showed an association with measures of limbic dopamine function, we sought to identify whether these dopamine associated nodes tended to overlap with nodes classified as combination hubs (as defined above using the rfMRI data).
The overlap of dopamine associated nodes and combination hubs was quantified using the Dice Similarity Coefficient: 25,26 A is the set of nodes in the dopamine associated subnetwork and B is the set of combination hub nodes. The Dice Coefficient was calculated for each of the 100 NBS thresholds (t=1.3-3.1) and then averaged to give a single 'true' score (main text Figure 2B part D). We then randomly selected an assortment of nodes, equal in number to the number of nodes present in the most leniently thresholded original network-based statistic subnetwork (main text Figure 2B part E). Next, we randomly deleted a node from this original assortment whenever the number of nodes in the 'true' subnetwork dropped as NBS threshold stringency increased (main text Figure 2B part F). This gave us 100 thresholds for this randomly generated subnetwork, and for each we calculated the Dice Coefficient with the same combination hubs, and then calculated a single mean 'random' Dice Coefficient as before. We repeated this procedure 10,000 times yielding 10,000 random Dice Coefficients (main text Figure 2B parts G-H), which allowed us to test the significance of the true Dice Coefficient (main text Figure   2B part I). This procedure was then repeated for each of the combination hub thresholds (10-40%), thereby giving a p-value for each hub threshold.
At some more lenient network-based statistic thresholds the network-based statistic defined dopamine associated networks contained all nodes of the SAL/DMN networks. In these cases, all nodes will overlap with the hub nodes, and so it is not valid to test if overlap is statistically significant. In these cases, we increased NBS stringency until the network no longer contained all the nodes in question.
We finally examined overlap between the dopamine associated nodes identified in Experiment 1 and those identified in Experiment 2. In this case we only compared overlap at NBS thresholds where both experiments showed the same number of dopamine associated nodes. We calculated the dice coefficient between the two networks, and compare it to a null distribution generated as before.

Software
Statistical analysis was undertaken in MATLAB 2016b and R 3.3.2.