Cocaine dependence and thalamic functional connectivity: a multivariate pattern analysis.

Cocaine dependence is associated with deficits in cognitive control. Previous studies demonstrated that chronic cocaine use affects the activity and functional connectivity of the thalamus, a subcortical structure critical for cognitive functioning. However, the thalamus contains nuclei heterogeneous in functions, and it is not known how thalamic subregions contribute to cognitive dysfunctions in cocaine dependence. To address this issue, we used multivariate pattern analysis (MVPA) to examine how functional connectivity of the thalamus distinguishes 100 cocaine-dependent participants (CD) from 100 demographically matched healthy control individuals (HC). We characterized six task-related networks with independent component analysis of fMRI data of a stop signal task and employed MVPA to distinguish CD from HC on the basis of voxel-wise thalamic connectivity to the six independent components. In an unbiased model of distinct training and testing data, the analysis correctly classified 72% of subjects with leave-one-out cross-validation (p < 0.001), superior to comparison brain regions with similar voxel counts (p < 0.004, two-sample t test). Thalamic voxels that form the basis of classification aggregate in distinct subclusters, suggesting that connectivities of thalamic subnuclei distinguish CD from HC. Further, linear regressions provided suggestive evidence for a correlation of the thalamic connectivities with clinical variables and performance measures on the stop signal task. Together, these findings support thalamic circuit dysfunction in cognitive control as an important neural marker of cocaine dependence.

Introduction

Cocaine dependence is a chronic, relapsing disorder (Yuferov et al. 2005). Previous studies have implicated deficits in cognitive control as a critical psychological factor contributing to continued drug use in dependent individuals (de Wit, 2009, Garavan and Hester, 2007, Li and Sinha, 2008, Porrino et al., 2007). Numerous imaging studies described cortical and subcortical dysfunctions in individuals addicted to cocaine or other stimulants (Aron and Paulus, 2007, Goldstein et al., 2009, Goldstein et al., 2007, Hanlon et al., 2009, Hanlon et al., 2011, Hester and Garavan, 2004, Kaufman et al., 2003, Moeller et al., 2005, Wesley et al., 2011). For instance, chronic cocaine users showed hypo-activation of the thalamus during a visual spatial attention task (Tomasi et al. 2007a). More recently, in a longitudinal study, we showed that decreased error-related activation of the thalamus predicts relapse and time to relapse to drug use in cocaine dependent individuals (Luo et al. 2013).

On the other hand, increased thalamic activations were observed in individuals with cocaine dependence and other addictive disorders when they were exposed to drug cues or situational factors related to drug use (Feldstein Ewing et al., 2010, Filbey et al., 2009, Franklin et al., 2009, Gozzi et al., 2011, Hermann et al., 2006, Jia et al., 2011, McClernon et al., 2009, Rose et al., 2007, Tomasi et al., 2007a, Wang et al., 2007, Weinstein et al., 2010). This seeming inconsistency may relate to functional heterogeneity of the thalamus, with some subregions mediating craving and others subserving control of craving. Understanding how thalamic cortical circuits respond to different psychological constructs and behavioral contexts will help elucidate the multifaceted etiologies of cocaine addiction.

Many studies have demonstrated altered cerebral functional connectivities in cocaine dependence (Cisler et al., 2013, Ding and Lee, 2013b, Gu et al., 2010, Kelly et al., 2011, Konova et al., 2013, Li et al., 2000, Narayanan et al., 2012, Tomasi et al., 2010, Verdejo-Garcia et al., 2012, Wilcox et al., 2011, Zhang et al., 2014). Reduced connectivities were observed in cocaine users between the thalamus and midbrain (Tomasi et al. 2010), putamen and ventral tegmental area (Gu et al. 2010), and anterior cingulate cortex (Verdejo-Garcia et al. 2012). The connectivity between the thalamus and midbrain (Tomasi et al. 2010) and ventral tegmental area (Gu et al. 2010) correlated negatively with the number of years of cocaine use. These studies suggest that cocaine use affects not only the activity but also connectivity of the thalamus. Investigating the functional connectivities of the thalamus may further uncover circuit-level disturbances in cocaine dependence.

Connectivity analysis has delineated functional subclusters of the thalamus in accordance with anatomy (Zhang et al., 2008, Zhang et al., 2010). Multivariate pattern analysis (MVPA), a data-driven technique, can provide new insights into cerebral functional organization (Haxby, 2012, Haynes and Rees, 2006, Norman et al., 2006, O'Toole et al., 2007, Pereira et al., 2009, Yang et al., 2012). Compared to univariate analyses, MVPA fully utilizes activity patterns across multiple variables and increases sensitivity in differential diagnostics (Norman et al., 2006, Pereira et al., 2009). MVPA has been used to predict behavioral variables and outcomes (Dosenbach et al. 2010), decode cognitive states (Brodersen et al., 2012, Freeman et al., 2011, Hassabis et al., 2009, Haxby et al., 2001, Haynes and Rees, 2005, Kamitani and Tong, 2005, Lee et al., 2013, Liu et al., 2011, Rissman et al., 2010), and identify patients with depression (Craddock et al., 2009, Fu et al., 2008, Zeng et al., 2012), functional dyspepsia (Liu et al. 2013), autism (Coutanche et al., 2011, Ecker et al., 2010), binge-eating (Weygandt et al. 2012b), mild cognitive impairment/Alzheimer's disease (Cuingnet et al., 2011, Desikan et al., 2009, Fan et al., 2008, Kloppel et al., 2008, Zhou et al., 2010), attention-deficit/hyperactivity disorder (Zhu et al. 2008), and schizophrenia (Ardekani et al., 2011, Shen et al., 2010), in a variety of behavioral contexts. However, to our knowledge, no studies have used MVPA to examine brain activity and connectivity in cocaine addiction.

Here, we applied MVPA to examine altered thalamic functional connectivity in participants with cocaine dependence (CD). Our earlier work of independent component analysis identified six independent networks of cognitive control from the stop signal task (Zhang and Li 2012). We sought to distinguish CD from healthy controls (HC) with MVPA of voxel-wise thalamic connectivity to these six task-related networks (please see Methods). We posited that, first, thalamic connectivity to these networks should distinguish CD from HC at an accuracy higher than comparison regions with similar volumes (number of voxels). Second, the voxels of which the connectivities to a set of independent components determine the “membership” should aggregate in clusters (as subnuclei) rather than distribute randomly. Thus, we have these specific aims: to investigate whether we could identify CD using only functional connectivity of the thalamus; to examine how the individual voxels of the thalamus distinguish CD from HC with functional connectivities for cognitive control; and to explore the correlation of these connectivities with clinical and performance variables.

Materials and methods

Subjects, informed consent, and assessment

One hundred recently abstinent participants (62 men) with cocaine dependence (CD) and one hundred age- and gender-matched healthy adult (HC) subjects (55 men) participated in this study (Table 1). CD met criteria for current cocaine dependence, as diagnosed by the Structured Clinical Interview for DSM-IV (First et al. 1995). Recent cocaine use was confirmed by urine toxicology screens. They were drug-free while staying in an inpatient unit prior to the current fMRI study. All subjects were physically healthy with no major medical illnesses or current use of prescription medications. None reported having a history of head injury or neurological illness. Other exclusion criteria included dependence on another psychoactive substance (except nicotine) and current or past history of psychotic disorders. Individuals with current depressive or anxiety symptoms requiring treatment or currently being treated for these symptoms were excluded as well. The Human Investigation committee at Yale University School of Medicine approved all study procedures, and all subjects signed an informed consent prior to study participation.

Table 1

Demographics of the subjects.

Subject characteristic	CD (n = 100)	HC (n = 100)	p-Value
Age (years)	40.3 ± 7.4	38.0 ± 10.6	0.08a
Gender (M/F)	62/38	55/45	0.39b
Years of alcohol use	17 ± 9.0	20 ± 10.2	0.02a
Years of Marijuana use	10 ± 4.2	1.0 ± 1.3	0.001a
Amount of average monthly cocaine use (gm) in the prior year	16.9 ± 25.8	N/A	N/A
Amount per use in grams	1.0 ± 1.2	N/A	N/A
Days of cocaine use in the prior month	15.1 ± 8.8	N/A	N/A
Years of cocaine use	17.5 ± 8.3	N/A	N/A
Days abstinent prior to scan	18.2 ± 6.1	N/A	N/A

Note: values are mean ± S.D.

Two-tailed two-sample t test.

χ2 test.

CD's were assessed with the Beck Depression Inventory (Beck et al. 1961) and the State-Trait Anxiety Inventory (Speilberger et al. 1970) at admission. The mean (± SD) BDI (12.2 ± 8.9) and STAI state (36.2 ± 10.7) and trait (40.7 ± 11.2) scores were within the range reported previously for individuals with cocaine dependence (Falck et al., 2002, Karlsgodt et al., 2003, Lopez and Becona, 2007, Rubin et al., 2007). Cocaine craving was assessed with the cocaine craving questionnaire, brief version (CCQ-Brief), for all participants every two to three days (Sussner et al. 2006). The CCQ-Brief is a 10-item questionnaire, abbreviated from the CCQ-Now (Tiffany et al. 1993). It is highly correlated with the CCQ-Now and other cocaine craving measures (Sussner et al. 2006). Each item was rated on a scale from 1 to 7, with a higher total score (ranging from 10 to 70) indicating greater craving. CDs averaged CCQ scores of 20.0 ± 7.8 across all assessments and 18.0 ± 5.5 on the day or within 2–3 days of the scan. The majority of CDs received all of these assessments, and subjects with missing data were not included in the correlation analyses (see below).

Behavioral task and scan procedures

We employed a stop-signal task (SST) as detailed in our previous studies (Duann et al., 2009, Li et al., 2006, Zhang and Li, 2012). Briefly, there were two trial types: “go” and “stop,” randomly intermixed. A small dot appeared on the screen to engage attention at the beginning of a go trial. After a randomized time interval (fore-period) between 1s and 5s, the dot turned into a circle, prompting the subjects to quickly press a button. The circle vanished at button press or after 1s had elapsed, whichever came first, and the trial terminated. A premature button press prior to the appearance of the circle also terminated the trial. Three quarters of all trials were go trials. In a stop trial, an additional “X,” the “stop” signal, appeared after the go signal. The subjects were instructed to withhold button pressing upon seeing the stop signal. Likewise, a trial terminated at button press or when 1s had elapsed since the appearance of the stop signal. The stop trials constituted the remaining one quarter of the trials. There was an inter-trial-interval of 2 s. The stop signal delay (SSD) started at 200 ms and varied from one stop trial to the next according to a staircase procedure, increasing and decreasing by 64 ms each after a successful and failed stop trial (De Jong et al., 1990, Levitt, 1971). Subjects were instructed to respond to the go signal quickly while keeping in mind that a stop signal could appear occasionally. Each subject completed four 10 min runs of the task after a practice session outside the scanner. With the staircase procedure, we anticipated that the subjects would succeed in withholding their response in approximately 50% of stop trials.

Analyses of behavioral data

We computed a critical SSD that represents the time delay between go and stop signals that a subject would need to succeed in 50% of the stop trials (Levitt 1971). Specifically, SSDs across trials were grouped into runs, with each run defined as a monotonically increasing or decreasing series. We derived a mid-run estimate by taking the middle SSD (or average of the two middle SSDs when there was an even number of SSDs) of every second run. The critical SSD was computed by taking the mean of all mid-run SSDs. It was reported that, except for experiments with a small number of trials (less than 30), the mid-run estimate was close to the maximum likelihood estimate of X50 (50% positive response; i.e., 50% stop success or SS in the SST, (Wetherill et al. 1966)). The stop signal reaction time (SSRT) was computed by subtracting the critical SSD from the median go trial reaction time (RT) (Logan 1994).

We computed the fore-period effect as an index of motor preparedness during the SST (Li et al., 2006, Li et al., 2005). Briefly, longer fore-periods are associated with faster RTs (Bertelson and Tisseyre, 1968, Woodrow, 1914). RT was compared between go trials with a fore-period between 3 and 5 s and between 1 and 3 s, and the effect size of RT difference was defined as fore-period effect. It is also known that in a RT task, the RT of a correct response is prolonged following an error, compared with other correct responses, and this prolonged RT is thought to reflect error monitoring (Rabbitt 1966). We thus computed the RT difference between the go trials that followed a stop error (SE) and those that followed another go trial, and termed the effect size of this RT difference “post-error slowing” (PES) (Hu et al., 2015, Ide et al., 2015).

Imaging protocol

Conventional T1-weighted spin echo sagittal anatomical images were acquired for slice localization using a 3T scanner (Siemens Trio). Anatomical images of the functional slice locations were next obtained with spin-echo imaging in the axial plane parallel to the AC–PC line with TR = 300 ms, TE = 2.5 ms, bandwidth = 300 Hz/pixel, flip angle = 60°, field of view = 220 mm × 220 mm, matrix = 256 × 256, 32 slices with slice thickness = 4 mm and no gap. Functional, blood oxygen level-dependent (BOLD) signals were then acquired with a single-shot gradient echo echo planar imaging (EPI) sequence. Thirty-two axial slices parallel to the AC–PC line covering the whole brain were acquired with TR = 2000 ms, TE = 25ms, bandwidth = 2004 Hz/pixel, flip angle = 85°, field of view = 220 mm × 220mm, matrix = 64 × 64, 32 slices with slice thickness = 4 mm and no gap.

Spatial preprocessing

Data were analyzed with Statistical Parametric Mapping (SPM8, Wellcome Department of Imaging Neuroscience, University College London, U.K.). Images from the first five TRs at the beginning of each trial were discarded to enable the signal to achieve steady-state equilibrium between RF pulsing and relaxation. Standard image preprocessing was performed. Images of each individual subject were first realigned (motion corrected) and corrected for slice timing. A mean functional image volume was constructed for each subject per run from the realigned image volumes. These mean images were co-registered with the high-resolution structural image and then segmented for normalization with affine registration followed by nonlinear transformation (Ashburner and Friston, 1999, Friston et al., 1995). The normalization parameters determined for the structure volume were then applied to the corresponding functional image volumes for each subject. Finally, the images were smoothed with a Gaussian kernel of 8 mm at Full Width at Half Maximum.

Independent component analysis

Preprocessed time series were analyzed with a group ICA algorithm (GIFT, http://icatb.sourceforge.net/, version 1.3i) to identify spatially independent and temporally coherent networks (Calhoun et al., 2001, Calhoun et al., 2009). ICA identifies distinct groups of brain regions with the same temporal pattern of hemodynamic signal change.

The data of CD (n = 100) and HC (n = 100) individuals were reduced through principal component analysis (Calhoun et al. 2001) and separated into 30 maximally independent components (ICs) with an infomax algorithm (Bell and Sejnowski 1995). The dimensionality was determined by the modified minimal description length (MDL) criteria as implemented in GIFT and averaged across subjects (Li et al. 2007). A time course for each IC and its corresponding spatial map was obtained. This analysis was repeated 20 times with ICASSO to assess the repeatability of ICs (Himberg et al. 2004). Finally, component time courses and spatial maps, which captured individual differences in the expression of the ICA-derived component, were back reconstructed for each participant (Calhoun et al., 2001, Meda et al., 2009). Following our previous study (Zhang and Li 2012), six ICs, including a motor cortical network for motor preparation and execution, a right fronto-parietal network for attentional monitoring, a left fronto-parietal network for response inhibition, a midline cortico-subcortical network for error processing, a cuneus-precuneus network for behavioral engagement, and a default-mode network (DMN) for self-referential processing, were included for further analyses.

We normalized back-reconstructed spatial maps of each IC into z-scores (Beckmann et al. 2005) and averaged them across runs for each participant. A one-sample t test was applied to the “z maps” across all participants to define significant brain regions associated with each IC.

Multivariate pattern analysis (MVPA) of the thalamic connectivity

Fig. 1 illustrates the analyses step-by-step with details described as follows. A multivariate pattern classifier was trained to identify CD from HC, followed by cross-validation using “leave-one-out”. Briefly, one CD or HC was left out as testing data while the remaining subjects were used to train the classifier. This procedure was repeated until each of the 100 CD and 100 HC was selected once as the validation data. A thalamus mask (698 voxels) was obtained from the Automated Anatomical Labeling atlas (Tzourio-Mazoyer et al. 2002). For each participant, we extracted z values of thalamic voxels for all 6 ICs as a classification feature vector, resulting in a total of 698 × 6 = 4188 features. The feature space was refined by retaining the most discriminating features, as assessed by feature-wise two-sample t tests of CD versus HC using training data only; features that differed between CD and HC at p < 0.05, uncorrected were included in the feature set for classification. Since a distinct set of features were selected from each set of training data, the final feature set differed slightly from iteration to iteration. We thus used the frequency with which each feature was included in the classification sets across all iterations to represent its discriminative power.

Fig. 1

A flow chart of the data analytic procedures. In step 0, the fMRI data of 100 CD and 100 HC were analyzed by independent component analysis (ICA) to generate 30 networks. In step 1, six task-related networks were selected for further analysis. A thalamus mask was applied to generate feature set including only voxels within the thalamus for multivariate pattern analysis (MVPA). In step 2, MVPA was applied using leave-one-out cross-validation to obtain true accuracy rates in group classification. Classification features were selected and support vector machine (SVM) classifier was trained based on the training data set. During testing, the same features from the testing data set were applied to the trained SVM classifier to obtain classification results (CD or HC). This procedure was repeated until each of the 100 CD and 100 HC was selected once as the validation data. The mean accuracy rate was computed to index overall accuracy. Finally, in step 3, the selected thalamic voxels from each ICA network were analyzed across 200 leave-one-out cross validation runs.

We employed a support vector machine (SVM) classifier with a linear kernel function for classification, as in previous work of MVPA (Dosenbach et al., 2010, Liu et al., 2011, Weygandt et al., 2012a, Zeng et al., 2012). Linear SVM helps to weigh down the effect of noisy features that are highly correlated (Pereira et al. 2009). According to the results of cross-validation, we used the generalization rate, sensitivity and specificity to quantify the performance of our classifier. The generalization rate represents the overall prediction rate of CD and HC combined, whereas the sensitivity and specificity each characterizes the prediction rates of CD and HC.

We used permutation tests to examine the statistical significance of the observed classification accuracy (Golland and Fischl, 2003, Zeng et al., 2012). In each permutation, the group label of the training data were randomly permuted prior to training and a generalization rate was computed. The permutation was repeated 1000 times to obtain a distribution of generalization rates, against which the accuracy in classification of the veridically labeled data was tested.

In addition to the thalamus, we examined other AAL masks with a similar volume for comparison. These regions included the paracentral lobule (695 voxels), the caudate (658 voxels), middle part of orbital frontal gyrus (694 voxels), parahippocampal gyrus (728 voxels), and inferior occipital gyrus (638 voxels), all with a volume close to the thalamus (698 voxels).

Results

Behavioral performance

Table 2 shows behavioral measures of the SST. Both CD and HC succeeded in about half of the stop trials, indicating success of the staircase procedure in tracking their performance. Compared to HC, CD showed lower go response rate (p = 0.004; two-tailed two-sample t test) and less PES (p = 0.02). CD also showed prolonged SSRTs, but the difference did not reach statistical significance (p = 0.11). The latter likely resulted from an under-estimation of the SSRT because the “RT” of a larger number of go error trials could not be considered in CD (Verbruggen et al. 2013). The two groups were otherwise not different in behavioral performance.

Table 2

SST performance.

	SSRT (ms)	FP effect (effect size)	Median go RT (ms)	%go	%stop	PES (effect size)
CD (n = 100)	231 ± 51	1.92 ± 1.45	589 ± 113	97.1 ± 6.2	52.6 ± 3.9	1.20 ± 1.87
HC (n = 100)	220 ± 44	2.05 ± 1.55	618 ± 110	99.0 ± 1.8	52.7 ± 3.2	1.79 ± 1.79
p-Value⁎	0.11	0.56	0.07	0.004	0.88	0.02

Note: All values are mean ± standard deviation; CD: individuals with cocaine dependence; HC: healthy controls; SSRT: stop signal reaction time; FP: fore-period; RT: reaction time; %go: percentage of go response trials; %stop: percentage of stop success trials; PES: post-error slowing.

p-Value based on 2-tailed 2-sample t test.

Altered thalamic connectivity in cocaine dependence

We verified the spatial congruency of the independent components identified from the current data set with those identified from our earlier work (Zhang and Li 2012) with spatial linear correlations. These six networks are motor cortical network (IC024, r2 = 0.55), right frontoparietal network (IC009, r2 = 0.46), left frontoparietal network (IC027, r2 = 0.53), midline cortico-subcortical network (IC029, r2 = 0.17), cuneus-precuneus network (IC003, r2 = 0.65), and DMN (IC022, r2 = 0.44), as shown in Fig. 2.

Fig. 2

Six task-related independent component networks identified from ICA (p < 0.000001, FWE corrected), shown in sagittal, coronal and axial views.

A total of 513 ± 10 features were selected at p < 0.05, uncorrected. In Fig. 3, we show the clusters of voxels of which the connectivity with the six ICs form the feature space for classification. With connectivity to the motor cortical network (IC024), a cluster in the area of bilateral ventrolateral nuclei (VL) possibly including mediodorsal nuclei (MD) distinguished CD from HC. Compared to HC, CD showed significantly lower functional connectivity between this cluster and the motor cortical network (p = 0.0002, two-sample t test).

Fig. 3

(A) Identified thalamic clusters, of which the connectivities to each of the six ICA networks were used in the classifier to distinguish CD from HC. Only clusters with at least 10 voxels were shown. (B) Z scores of each cluster for CD and HC with mean ± standard error. p-Value reflects two sample t test result of the Z scores of each cluster. With the exception for the connectivity of the red cluster to DMN (IC022), all clusters showed significant differences between the CD and HC with correction for multiple comparisons. We further tested whether the Z scores of CD or HC were significantly different from zero, with * indicating a p < 0.05, corrected for multiple comparisons, and with ** indicating a p < 0.0001, uncorrected.

With connectivity to the right frontoparietal network (IC009), three clusters, each in the area of the right and left ventroposterior lateral nuclei (VPL) and pulvinar (PUL), and left VL nucleus distinguished CD from HC. Compared to HC, CD showed significantly higher functional connectivity in the right VPL/PUL (p = 0.0001) and lower connectivity in the left VPL/PUL (p = 0.0002) and left VL (p = 0.0002).

With connectivity to the left frontoparietal network (IC027), clusters in the area of bilateral VL and right MD distinguished CD from HC. Compared to HC, CD showed significant lower connectivity in bilateral VL (all p's < 0.00001) and higher connectivity in right MD (p = 0.0009).

With connectivity to the midline cortical-subcortical network (IC029), clusters in the area of bilateral MD distinguished CD from HC. Compared to HC, CD showed higher connectivity in bilateral MD (all p's < 0.0001).

With connectivity to the cuneus-precuneus network (IC003), clusters in the area of bilateral VPL distinguished CD from HC. Compared to HC, CD showed higher connectivity to VPL (all p's < 0.000001).

With connectivity to the DMN (IC022), clusters in the right VL and bilateral MD distinguished CD from HC. Compared to HC, CD showed higher connectivity in these clusters (p = 0.002 and p = 0.001).

Aggregation index analysis

To examine how well the thalamic voxels that distinguished CD and HC aggregate in clusters, we computed an aggregation index (AI) for each independent component for the thalamus and comparison regions (Table 3). The AI was computed as the total number of voxels that appear in a cluster ≥ 10 voxels divided by the total number of voxels. Thus, a higher AI indicated that more voxels aggregate in clusters. Thalamus showed a significantly higher AI than comparison regions, including the caudate (p = 0.03), middle orbital frontal gyrus (p = 0.02), parahippocampal gyrus (p = 0.008), and inferior occipital gyrus (p = 0.03), but not paracentral gyrus (p = 0.14). The results were similar when we redefined the cluster set size to ≥ 20 voxels (the caudate: p = 0.039; middle orbital frontal gyrus: p = 0.055; parahippocampal gyrus: p = 0.020; inferior occipital gyrus: p = 0.119; paracentral gyrus: p = 0.270), or ≥ 40 voxels (the caudate: p = 0.001; middle orbital frontal gyrus: p = 0.025; parahippocampal gyrus: p = 0.059; inferior occipital gyrus: p = 0.050; paracentral gyrus: p = 0.426). In no cases were the AIs of any of the comparison regions higher than that of the thalamus.

Table 3

Aggregation index (AI) of the thalamus and comparison regions.

	IC003	IC009	IC022	IC024	IC027	IC029	Mean	p-Value
Thalamus	0.81	0.41	0.67	0.64	0.80	0.61	0.66	N/A
Comparison regions
Caudate	0.40	0.46	0.45	0.00	0.00	0.67	0.33	0.03
Paracentral	0.00	0.00	0.00	0.95	0.69	0.55	0.37	0.14
Occipital_Inf	0.87	0.55	0.00	0.00	0.00	0.00	0.24	0.03
Frontal_Orb_Mid	0.16	0.40	0.26	0.00	0.60	0.61	0.34	0.02
Parahippocampal G	0.54	0.23	0.15	0.20	0.00	0.60	0.29	0.008

Note: p-value: paired t test of the AI of the thalamus, as compared to each of the other regions, across all six components.

Correlation with performance outcomes and clinical assessments

Thalamic cortical connectivity is critical to cognitive control in the stop signal task (Ide and Li 2011). We previously showed that, compared to HC, CD showed prolonged SSRT and diminished PES (Li et al. 2008). Thus, we examined whether thalamic connectivity to each of the independent component of cognitive control is related to SSRT and PES. We derived the z scores of connectivity (connectivity strength) between each independent network and corresponding cluster for individual participants. Pearson regression across subjects showed positive correlations between SSRT and right MD connectivity to the left frontoparietal network in CD (p = 0.007, r = 0.27) but not HC (p = 0.39, r = 0.09) and between PES and left MD connectivity to the DMN network in HC (p = 0.05, r = 0.20) but not CD (p = 0.68, r = − 0.04). Conversely, a negative correlations was observed between SSRT and left MD connectivity to the DMN in HC (p = 0.01, r = − 0.25) but not CD (p = 0.55, r = 0.06). Further, the linear regressions between SSRT and left MD connectivity to the DMN were significantly different in slope between CD and HC (p < 0.05).

Left VL connectivity to left frontoparietal network positively correlated with CCQ score on the scan day (p = 0.006, r = 0.27) and left VPL/PUL connectivity to right frontoparietal network positively correlated to the amount of cocaine (grams) per use (p = 0.05, r = 0.2). Conversely, right VPL connectivity to the cuneus-precuneus network negatively correlated with the duration of abstinence (days) prior to scan (p = 0.01, r = − 0.26) and left VL connectivity to the right frontoparietal network negatively correlated with the years of cocaine use (p = 0.05, r = − 0.2). None of the connectivity strengths was correlated with days of cocaine use (all p's > 0.1) or daily amount of cocaine use (grams) in the month prior to admission (all p's > 0.07), or years of cocaine use (all p's > 0.05), in CD. The connectivity strengths also did not correlate with years of alcohol use (all p's > 0.06), years of marijuana use (all p's > 0.1), BDI score (all p's > 0.1), STAI state anxiety score (all p's > 0.1), STAI trait anxiety score (all p's > 0.07).

None of the correlations with outcome measures of the SST or with clinical variables were significant, when corrected for multiple comparisons.

Classification results

With leave-one-out cross-validation, the linear SVM classifier achieved a mean accuracy of 72% (76% for CD and 68% for HC, and F1 score = 73%) across 200 iterations (p = 0.001, permutation test; Fig. 4). On the other hand, SVM classification accuracies of the control regions were all less than 60%, and p-values of permutation tests were not significant (Table 4). Further, compared to the control regions, the thalamus showed significantly higher accuracy rates in classification (p = 0.004; two-sample t test; Table 4).

Fig. 4

The permutation distribution of the estimated generalization rate using the linear support vector machine classifier (repetition times = 1000). GR0 represented the real generalization rate as obtained from the current data set. Only 1 sample showed a higher generalization rate than the real data in 1000 runs (p = 0.001).

Table 4

MVPA accuracy rates of the thalamus and comparison regions.

	Total features	Used features (at p < 0.05)	MVPA performance of CD					Permutation test
TPR	TNR	ACC	F1	AUC ROC
Thalamus	4188	513 ± 10	76%	68%	72%	73%	73%	p < 0.001
Comparison regions
Caudate	3948	338 ± 7	60%	53%	56.5%	58%	57%	p ≤ 0.16
Paracentral	4170	334 ± 8	58%	59%	58.5%	58%	58%	p ≤ 0.08
Occipital_Inf	3828	308 ± 6	53%	53%	53%	53%	53%	p ≤ 0.36
Frontal_Orb_Mid	4164	226 ± 7	52%	54%	53%	53%	53%	p ≤ 0.43
Parahippocampal G	4368	298 ± 6	53%	52%	52.5%	53%	53%	p ≤ 0.32

Note: Values under “used features at p < 0.05” are mean ± standard deviation; CD: individuals with cocaine dependence; HC: healthy controls; Paracentral: paracentral gyrus; Occipital_Inf: inferior occipital gyurs; Frontal_Orb_Mid: middle part of orbital frontal gyrus; Parahippocampal G: parahippocampal gyrus. TPR: true positive rate; TNR: true negative rate; ACC: accuracy; F1: F1 score; AUC ROC: Area under the Curve of receiver operating characteristic curve.

Discussion

With MVPA of thalamic connectivity to the six independent functional networks, the current study demonstrated that CD can be distinguished from HC with a higher (72%) accuracy, relative to comparison regions with similar volumes (53% to 58%). Of note, because the testing and training samples are distinct for each validating iteration, this accuracy rate reflects a true estimate of the power of thalamic connectivity in distinguishing CD from HC. Moreover, a higher number of the identified voxels aggregate in clusters in the thalamus, in contrast to the comparison regions, indicating that these connectivities are non-random and reflect the altered functional architecture of the thalamus in cocaine addiction. Together, these results suggest that thalamic connectivities to the functional networks of cognitive control represent a useful circuit maker of cocaine dependence.

A recent study employed a linear classifier to distinguish cocaine-dependent from non-drug-abusing men on the basis of independent component networks of a SST (Elton et al. 2012). The authors identified 15 task-related networks attributed to motor, visual, cognitive and affective processes and used the statistical associations of these components with trial types for group classification. Along with this earlier work, the current results support a critical role of altered cognitive control in characterizing cocaine addiction and, by describing the neural substrates, extend the literature in new directions.

Altered thalamic connectivity for cognitive control

With connectivity to the motor cortical, bilateral frontoparietal, and default-mode networks, clusters in the area of the ventrolateral thalamic nuclei (VL) distinguished CD and HC. Compared to HC, CD showed decreased VL connectivity to the motor cortical and frontoparietal networks and increased connectivity to the DMN. CD and HC were also distinguished by mediodorsal nucleus (MD) connectivity to the left frontoparietal, midline cortico-subcortical, and DMN. Compared to HC, CD showed increased connectivity to all three functional networks. These findings are broadly consistent with cognitive motor dysfunction as a result of cocaine misuse (Elton et al., 2012, Penberthy et al., 2010, Volkow et al., 2001).

The frontoparietal networks play a ubiquitous role in human cognitive functioning (Cole et al., 2012, Corbetta, 1998, Iidaka et al., 2006, Kong et al., 2013, Ptak, 2012, Vincent et al., 2008, Whitman et al., 2013, Zhang and Li, 2012). Frontoparietal dysfunction has been identified in cocaine addicts in multiple experimental domains (Aron and Paulus, 2007, Barros-Loscertales et al., 2011, Garavan et al., 2008, Kaufman et al., 2003, Tomasi et al., 2007b). The current findings showed that connectivity of the anterior and lateral pulvinar (PUL) and the ventroposterior lateral nuclei (VPL) to the right frontoparietal network and connectivity of the MD and left VL to the left frontoparietal network distinguished CD from HC. The right and left frontoparietal networks are each involved in attentional monitoring and response inhibition in the SST (Zhang and Li 2012). Thus, the findings are consistent with the role of the PUL in attentional processing and binding of perceptual information, as shown in imaging (Arend et al., 2008a, Komura et al., 2013, Michael and Desmedt, 2004, Saalmann et al., 2012, Wilke et al., 2010) and lesion (Arend et al., 2008b, Snow et al., 2009) studies. On the other hand, the right VL was associated with response inhibition (Congdon et al., 2010, Schmaal et al., 2013), with activation negatively correlated with the SSRT (Congdon et al. 2010), an outcome measure of response inhibition in the SST (Eagle et al., 2008, Li et al., 2006).

Previous work implicates altered error processing and error-related learning in CD (Franken et al., 2007, Hester et al., 2007, Li et al., 2010, Luo et al., 2013, Madoz-Gurpide et al., 2011, Sokhadze et al., 2008, Vadhan et al., 2008). Consistently, the current findings showed altered connectivity of the dorsal MD to midline cortical-subcortical network, which is involved in error processing in the SST (Hendrick et al., 2010, Ide and Li, 2011, Zhang and Li, 2012). Further, error-related activations elevated in the MD during early compared to late cocaine abstinence, whereas decreased error-related MD activation later during abstinence predicted relapse and time to relapse (Li et al., 2010, Luo et al., 2013). Although it remains premature to link changes in regional activation and connectivity, these findings support thalamic cortical mechanisms as an important etiological process of cocaine addiction.

A more recent study suggested reduced functional connectivity between the MD and the anterior cingulate cortex (ACC) in CD compared to HC, which seems inconsistent with the current findings (Verdejo-Garcia et al. 2012). However, the latter study examined resting state functional connectivity while the current work focused on task-related connectivity. Importantly, Verdejo-Garcia et al. observed reduced connectivity between the MD and rostral ACC, while the ACC of the midline cortical-subcortical network sits at a dorsal location. Dorsal and rostral ACC play different roles in cognitive control and connect with different brain regions (Margulies et al. 2007). Both ACC regions respond to errors, but the dorsal ACC may also play an evaluative role by assessing feedbacks for cognitive control (Bush et al., 2002, Polli et al., 2005, Spunt et al., 2012, Taylor et al., 2006, Van Veen and Carter, 2002).

Our previous studies suggested that activity of the cuneus-precuneus network may reflect behavioral engagement of the participants (Zhang and Li, 2010, Zhang and Li, 2012). Here, the connectivity between the VPL nucleus and cuneus-precuneus network was altered by cocaine use, with less negative connectivity in CD than HC, in accord with previous reports of anatomical connectivity between the precuneus and VPL nucleus (Cavanna and Trimble, 2006, Parvizi et al., 2006, Schmahmann and Pandya, 1990, Yeterian and Pandya, 1993).

The DMN refers to a set of brain regions that decreases activations during cognitive challenges and increases activations during resting, mind wandering, and self-referential processing (Buckner et al., 2008, Raichle and Snyder, 2007). An earlier study (Metzger et al. 2010) with high field fMRI demonstrated coactivation of MD and the pregenual ACC, a key brain region of the DMN (Biswal et al., 2010, Buckner et al., 2008, Walter et al., 2009). Our findings suggested that the MD altered its connectivity to the DMN, consistent with recent findings of DMN dysfunction in drug addiction (Bednarski et al., 2011, Ding and Lee, 2013a, Ding and Lee, 2013b, Li et al., 2011, Ma et al., 2011, Sutherland et al., 2012). For instance, cocaine dependence was associated with altered error-preceding activity of the DMN, a deficit that correlated with the years of cocaine use (Bednarski et al. 2011) and could be remediated with methylphenidate (Matuskey et al. 2013). Here, right VL and bilateral MD connectivity to the DMN distinguished CD from HC.

Broad methodological implications

The current results speak to the utility of thalamic functional connectivity in relation to cognitive control in distinguishing CD from HC. As described earlier, the thalamus is heterogeneous in functions and mediates multiple cognitive and affective processes in addition to cognitive control. Thus, studies using other task-related or even resting state imaging data will provide additional information needed to characterize thalamic processes compromised in cocaine addiction and to fully determine group membership in an analogous data analytic scheme. In these future studies, one would be able to investigate whether the thalamic subnuclei as identified across different behavioral tasks or states would conform to similar anatomical boundaries and yet partake in various cognitive and affective processes by interacting with different functional neural networks.

Leave-one-out versus k-fold cross validation

Besides “leave-one-out” as employed in the current study, k-fold (usually 10-fold) cross validation is another commonly used cross validation approach. On one hand, because the training sets differ in only one observation, leave-one-out approach involves a higher variance compared with 10-fold cross validation. On the other hand, leave-one-out approach yields close to unbiased estimates of the test error, since each training set contains n − 1 observations, which is almost as many as the number of observations in the full data set (James et al. 2013). One previous study compared these two cross validation approaches using 13 different data sets and observed that they have very similar results (Cawley and Talbot 2003).

Limitations of the study and conclusions

A number of limitations need to be considered. First, none of the thalamic connectivities that distinguish CD from HC appeared to correlate with clinical variables or outcome measures of the SST, following correction for multiple comparisons. Therefore, the current work supports thalamic cortical dysfunction as a neural marker of cocaine addiction, but the functional significance is unclear. On the other hand, it is not uncommon that neural phenotypes do not readily associate with clinical characteristics in psychiatric conditions. This discrepancy may suggest the importance in going beyond clinical evaluations in understanding the etiology of mental illness. Second, we did not examine gender differences in the patterns of thalamic connectivity because of the moderate and unbalanced sample size. The enormity of the feature space generated in such analyses can also be extremely challenging. Third, ICA was performed before cross validation, which may result in bias in classification. On the other hand, the differences in group ICA results between a sample of 200 and 199 (leave-one-out) subjects should be very limited. Fourth, many studies have parcellated the thalamus on the basis of diffusion tensor imaging (Behrens et al., 2003, Bisecco et al., 2015, Duan et al., 2007, O'Muircheartaigh et al., 2015, Serra et al., 2014, Ye et al., 2013, Zhang et al., 2010). Future work is warranted to compare the findings from these different parcellation schemes. Fifth, this is a cross-sectional study. Whether and how thalamic functional connectivities may predict trajectories of drug use warrants more investigation.

In summary, MVPA demonstrated thalamic functional connectivities as a potential diagnostic circuit marker of cocaine addiction. The intricate patterns of thalamic subregional connectivities that distinguish cocaine-addicted from non-addicted individuals confirm the utility of this computational approach in psychiatric neuroscience research.

Disclosures

All authors report that they have no conflicts of interest over the past five years to report as related to the subject of the report. Dr. Potenza has consulted for Ironwood, Lundbeck, Boehringer Ingelheim, Shire, INSYS, Lightlake Therapeutics and RiverMend Health; has received research support from the National Institutes of Health, Veteran's Administration, Mohegan Sun Casino, the National Center for Responsible Gaming, Pfizer, and Psyadon Pharmaceuticals; has participated in surveys, mailings or telephone consultations related to drug addiction, impulse control disorders or other health topics; has consulted for gambling, legal and educational entities on issues related to impulse control and addictive disorders; provides clinical care in the Connecticut Department of Mental Health and Addiction Services Problem Gambling Services Program; has performed grant reviews for the National Institutes of Health and other agencies; has edited journals and journal sections; has given academic lectures in grand rounds, CME events and other clinical or scientific venues; and has generated books or book chapters for publishers of mental health texts.