Disrupted topological organization in whole-brain functional networks of heroin-dependent individuals: a resting-state FMRI study.

Neuroimaging studies have shown that heroin addiction is related to abnormalities in widespread local regions and in the functional connectivity of the brain. However, little is known about whether heroin addiction changes the topological organization of whole-brain functional networks. Seventeen heroin-dependent individuals (HDIs) and 15 age-, gender-matched normal controls (NCs) were enrolled, and the resting-state functional magnetic resonance images (RS-fMRI) were acquired from these subjects. We constructed the brain functional networks of HDIs and NCs, and compared the between-group differences in network topological properties using graph theory method. We found that the HDIs showed decreases in the normalized clustering coefficient and in small-worldness compared to the NCs. Furthermore, the HDIs exhibited significantly decreased nodal centralities primarily in regions of cognitive control network, including the bilateral middle cingulate gyrus, left middle frontal gyrus, and right precuneus, but significantly increased nodal centralities primarily in the left hippocampus. The between-group differences in nodal centralities were not corrected by multiple comparisons suggesting these should be considered as an exploratory analysis. Moreover, nodal centralities in the left hippocampus were positively correlated with the duration of heroin addiction. Overall, our results indicated that disruptions occur in the whole-brain functional networks of HDIs, findings which may be helpful in further understanding the mechanisms underlying heroin addiction.

Introduction

Drug addiction, a major social problem, appears to be a chronic brain disease that involves complex interactions between biological and environmental variables and is characterized by a compulsive drive to take drugs despite serious negative consequences [1]. Heroin users are a major proportion of drug addicts, especially in China [2]. Since the development of neuroimaging technologies, many studies have been concerned with the mechanisms underlying drug addiction [3], [4], [5].

Resting-state functional magnetic resonance imaging (RS-fMRI), a non-invasive imaging technique, has been widely used to explore the intrinsic functional organization of the human brain [6], [7], [8]. Several studies that used this technique investigated heroin-related changes in spontaneous brain activity [9], [10], [11], [12], [13] and suggested that heroin addiction is related to widespread functional abnormalities in many brain regions. These regions include the amygdala [13], anterior cingulate cortex (ACC) [10], hippocampus [13], insula [13], lingual gyrus [9], orbitofrontal cortex (OFC) [9], and temporal cortex [10]. In addition, functional connectivity alterations have been found in heroin-dependent individuals (HDIs). Ma et al. [12] indicated that heroin users showed increases in functional connectivity between the nucleus accumbens and the ventral/rostral ACC and decreases in connectivity between the prefrontal cortex and the OFC as well as between the prefrontal cortex and the ACC. Liu et al. [11] detected abnormal connectivity between the prefrontal cortex, ACC, ventral striatum, insula, amygdala and hippocampus in heroin users. However, to date no study has considered heroin-related whole-brain functional networks during the resting-state.

Graph theory analysis provides a powerful tool for characterizing topological organization, including identifying global and nodal properties in whole brain functional networks. It has been applied to the study of normal brains [14], [15] and of various brain-related diseases, such as Alzheimer's disease [16], epilepsy [17], [18], depression [19], and schizophrenia [20], [21]. Although two previous studies [11], [22] explored brain functional networks in heroin addiction patients using graph theory analysis, these studies focused on regional functional connectivity [11] or on four specific circuits (control, reward, motivation/drive and memory) [22]. However, what heroin-related alterations occur in the whole-brain functional networks remains unknown. Because previous studies [9], [10], [11], [12], [13] have indicated that the brain connectivity alterations in HDIs are widespread, in this research we attempted to analyze the topological properties of whole-brain functional networks in HDIs based on graph theory.

In this study, we constructed brain functional networks with RS-fMRI data for a group of HDIs and a group of controls, and compared the topological organization of their brain networks using graph theory. In addition, considering the duration of heroin addiction as a vital clinical variable in understanding the effects of heroin on functional abnormalities [9], [23], [24], we analyzed the correlations between the altered network parameters and the duration of heroin addiction.

Materials and Methods

Subjects

We recruited seventeen heroin-dependent individuals (HDIs: 15 M/2 F, aged 26–50 years, mean ± SD = 36.29±6.86 years, right-handed) from the Addiction Medicine Division of Guangdong No. 2 Provincial People's Hospital. The HDIs were screened using the Structured Clinical Interview (SCID-IV) for the Statistical Manual of Mental Disorders, Fourth Edition (DSM-IV), to confirm the diagnosis of heroin dependence. Urine tests with a positive finding for heroin use were acquired before enrolling in the treatment program. None of the HDIs had used any other types of drugs according to a laboratory report and an interview conducted by a clinical psychologist (30 years of clinical experience). They were hospitalized for 6–7 days before RS-fMRI scanning took place, none of the HDIs used heroin, as confirmed by the medical personnel responsible for their care. All of HDIs except two were under daily methadone maintenance treatment at the time of study. In addition, we recruited fifteen age- and gender-matched normal controls (NCs: 12 M/3 F, aged 20–46 years, mean ± SD = 31.27±8.10 years, right-handed) as the normal controls. Table 1 lists the demographic details of all the volunteers in this study. The detailed clinical descriptions for each of HDIs are listed in Table S1.

Table 1

Demographic information for the heroin-dependent individuals (HDIs) and the normal controls (NCs) in the present study.

Characteristics	HDIs (n = 17)	NCs (n = 15)	p-value
Female/Male	2/15	3/12	0.645a
Age (years)	36.29±6.86	31.27±8.10	0.067b
Range (years)	26–50	20–46
Education (years) Range (years)	10.24±3.21 2–15	11.07±4.11 5–17	0.526b
Head motion
Translation (mm)	0.126±0.051	0.106±0.045	0.251b
Rotation (mm)	0.149±0.120	0.126±0.063	0.512b
Nicotine (Median, No. cigarette/day)	20 (0–40)	20 (0–40)	0.116b
Heroin use (years) Range (years)	9.21±5.28 1–19	N/A
Heroin dosage (g/day) Range (g/day)	0.80±0.54 0.1–2.0	N/A
Dosage of methadone (g/day)	39.41±20.45	N/A

The duration of heroin usage means the period from the time of their initial heroin use to the time of their seeking medical attention.

a Fisher's exact test.

t-test.b Two sample

Neither the HDIs nor NCs had any history of neurological illness or head injury or had been diagnosed with schizophrenia or an affective disorder according to their past medical history. This study was approved by the Research Ethics Review Board of the Southern Medical University in Guangzhou of China. Informed written consent was obtained from each subject prior to the MRI scanning.

Data acquisition

MRI data were obtained on a 1.5T Philips Achieva Nova Dual MR scanner in the Department of Medical Imaging, Guangdong No. 2 Provincial People's Hospital. The RS-fMRI data were obtained using a T2*-weighted gradient-echo echo-planar imaging (EPI) sequence with the following parameters, TR = 2000 ms, TE = 50 ms, flip angle = 90°, matrix = 64×64, FOV = 230×230 mm2, thickness/gap = 4.5/0 mm, 22 axial slices covering the whole brain, 240 volumes obtained in about 8 min. During the RS-fMRI scanning, all lights in the scanner room were switched off, and the subjects were instructed to close their eyes, to keep still, not to think systematically about anything, and not to fall asleep. In addition, we acquired 3D high resolution brain structural images using a T1-weighted 3D turbo-gradient-echo sequence (TR = 25 ms, TE = 4.1 ms, flip angle = 30°, matrix = 256×256, FOV = 230×230 mm2, thickness = 1.0 mm, and 160 sagittal slices).

Data preprocessing

All the MRI data were processed using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/) and DPARSF_V2.0 (http://www.restfmri.net/forum/index.php) [25]. For each subject, we first removed the first 10 volume images from the RS-fMRI data for scanner stabilization and for the subject's adaptation to the environment, leaving 230 volumes for further analysis. Then we performed slice timing to correct for the acquisition time delay between slices within the same TR, realignment to the first volume to correct the inter-TR head motions, spatial normalization to a standard MNI template and resampling to a voxel size of 3×3×3 mm3. No spatial smoothing was applied by following previous studies [26], [27], [28]. Finally, we performed band-pass filtering for each voxel in the frequency of 0.01–0.08 Hz to reduce low-frequency drift and high-frequency physiological noise. The RS-fMRI data for each subject were checked for head motion. No subject was excluded according to the criteria that the translation and rotation of head motion in any direction were not more than 1.5 mm or 1.5°.

Network analysis

Network construction

In order to construct brain functional networks for each subject, we applied an automated anatomical labeling (AAL) atlas [29] to parcellate the brain into 90 regions of interest (ROIs) (45 in each hemisphere). The names of the ROIs and their corresponding abbreviations are listed in Table S2. The time series for each ROI was calculated by averaging the signals of all voxels within that region and by linearly regressing out the following nine nuisance covariates: three translation and three rotation head motion parameters and the white matter, cerebrospinal fluid (CSF), and global mean signals. For each subject, we obtained a 90×90 correlation matrix by calculating the Pearson's correlation coefficient in the residual time courses between all ROI-pairs. This matrix contained both negative and positive values, we used the absolute value of each element as the inter-regional functional connectivity by following previous studies [19], [26], [30], [31]. Finally, this correlation matrix was thresholded into a binarized matrix with a sparsity value (the ratio between total number of edges and the maximum possible number of edges in a network). By taking each ROI as a node and the functional connectivity as an edge, we obtained a 90×90 connectivity matrix for each subject and analyzed the topological organization of the whole-brain functional networks according to graph theory.

Clearly, the choice of a sparsity value has a major effect on the topological organization of networks [32], [33]. By setting a specific sparsity as the threshold, we were able to ensure that the brain functional networks corresponding to each subject contained the same number of edges. In order to balance the prominence of the small-world attribute with an appropriate level of sparseness in the networks for all subjects, we determined the range of sparsity according to the following criteria: 1) the averaged degree (total number of edges divided by N/2, with N = 90 here, denoting the number of nodes) over all nodes of each network was larger than log(N) [34], [35]; and 2) the small-worldness of the network for each subject was larger than 1.1 [33], [34]. Thus, we determined the range of sparsity (0.05≤s≤0.36) in which the network for each subject holds the small-worldness property. Using different threshold values over the range of 0.05≤s≤0.36 and intervals of 0.01, we set the connectivity matrix into a series of binaried connectivity matrices for each subject and calculated the topological properties. The subsequent network analysis was based on the series of binarized connectivity matrices for each subject.

Network parameters

We described the global topological properties of the brain functional networks by using the following seven global network parameters: the clustering coefficient (), characteristic path length (), normalized clustering coefficient (), normalized characteristic path length (), small-worldness (), global efficiency (), and local efficiency (). Their expression and detailed descriptions are listed in Table S3.

Two nodal centrality metrics, nodal degree () and nodal efficiency (), were used to describe the nodal properties of brain functional networks. Their expressions and descriptions are also presented in Table S3.

Instead of selecting a single sparsity threshold, we used the integrated network parameters over the range of sparsity to detect the between-group differences in the topological parameters of the brain functional networks. The integrated global parameters were given by [33]:where the sparsity interval equals 0.01 and refers to any of the global parameters (, , , , , , and ) at a sparsity of . Similarly, the integrated nodal parameters can be calculated by [33]:where represents either nodal parameters (, ) of node at a sparsity of kΔs.

Hub identification

Hubs refer to highly connected nodes in the network [32]. Here, following the method used in previous studies [15], [36], we used nodal betweenness centrality () to determine the hub regions of the brain functional networks (for a detailed description, see Table S3. For each node, we first calculated its normalized nodal betweenness centrality as follows:where is the integrated nodal betweenness centrality of node in the network of subject , is the number of subjects in each group and is the number of nodes (here  = 90). Nodes satisfying the criterion of were considered to be the hubs of the brain functional networks [33]. Based on this criterion, we then identified the hubs of the brain functional networks separately for HDIs and NCs.

Statistical analysis

Between-group differences

Two sample t-tests were performed to assess differences in age, duration of education, cigarette smoking, and head motions between the heroin addict group and the control group using SPSS (version 17.0). We used Fisher's exact test to estimate the difference in gender between the two groups (SPSS, version 17.0). Significant between-group differences were determined at p<0.05 (two-tailed).

A nonparametric permutation test [37] was performed to determine significant differences in each integrated network metric (five global parameters and two nodal centrality metrics) between the two groups. Briefly, for each network metric, we first calculated the between-group difference in the mean values. To obtain an empirical distribution of the difference, we then randomly reallocated all the values into two groups and recomputed the mean differences between the two randomized groups (10,000 permutations). The limits of the 95th percentile for each empirical distribution were used as the critical values for a two-tailed test of whether the observed group differences could occur by chance. To check the statistical power for the between-group comparisons in nodal metrics, we also estimated the effect sizes (Cohen d) according to Cohen's definition [38].

We used a network-based statistic (NBS) approach [39] to detect differences in the inter-nodal functional connections between the HDIs and NCs. In brief, a primary cluster-defining threshold was first used to identify suprathreshold connections for which the size (i.e., number of edges) of any connected components was then determined. A corrected p-value was calculated for each component using the null distribution of the maximally connected component size, which was derived empirically using a nonparametric permutation approach. The detailed descriptions are provided in Text S1.

Correlations between network parameters and duration of heroin addiction

We analyzed the correlation between each of the network parameters and the duration of heroin addiction in HDIs using a multiple linear regression. The significance levels were set at p<0.05 (two-tailed).

Although the two groups were statistically matched for age, the heroin group was an average of 5 years older than the control group. To control for any potential age-related effect, all of the above analyses were repeated after removing the confounding effect of age using a multiple linear regression.

Results

Demographic information

Statistical comparisons showed no significant differences in gender, age, duration of education, cigarette smoking, and head motions between the heroin group and the control group (Table 1).

Global parameters

Fig. 1 shows the plots of the global parameters (, , , , , , and ) of the whole-brain functional networks changing with sparsity in both the HDIs and NCs. Fig. 1 also shows the comparisons for the values of , , , , and to be lower, but the values of and to be higher, in HDIs compared to NCs.

Figure 1

Global parameters of the brain functional networks for the heroin-dependent individuals (HDIs) and the normal controls (NCs) changing with the sparsity threshold.

The error bar represents the standard deviation of a parameter at a given sparsity across all subjects. The symbol (*) means that significant between-group difference in the given parameter was detected (p<0.05). Except for the sparsity range of 0.05≤sparsity≤0.09, no statistically significant between-group differences were detected for other values of sparsity. , clustering coefficient; , characteristic path length; , normalized clustering coefficient; , normalized shortest path length; , small-worldness; , global efficiency; , local efficiency.

Fig. 2 shows statistical comparisons of the integrated global parameters between HDIs and NCs. The HDIs exhibited significantly lower values for integrated small-worldness (p = 0.035) and the integrated normalized clustering coefficient (p = 0.049) compared to the controls. However, we found no significant between-group difference in any of the integrated parameters , , , , and .

Figure 2

Bar plots of the differences in the integrated global topological parameters of brain functional networks between the heroin-dependent individuals (HDIs) and the normal controls (NCs).

The symbol (*) indicates significant between-group differences in the integrated normalized clustering coefficient (p = 0.049) and integrated small-worldness (p = 0.035). , integrated clustering coefficient; , integrated characteristic path length; , integrated normalized shortest path length; , integrated global efficiency; , integrated local efficiency.

Nodal parameters

Table 2 lists the brain regions that showed a significant difference in any of nodal centrality metrics ( and ) of the brain functional networks between HDIs and NCs (p<0.05, uncorrected). We found that in HDIs, the nodal centrality metrics were significantly decreased in six brain regions, the bilateral middle (dorsal) cingulate gyrus (MCG.L/R), left middle frontal gyrus (MFG.L), left inferior temporal gyrus (ITG.L), right precuneus (PCUN.R), and right thalamus (THA.R), most of which belong to the cognitive control network [40], [41], [42]; but significantly increased in three brain regions, the left hippocampus (HIP.L), left inferior occipital gyrus (IOG.L), and left lingual gyrus (LING.L). These regions were rendering plotted on a cortical surface map and are shown in Fig. 3a. The result of effect sizes presented in Table 2 indicated high statistical power of the between-group comparions in nodal parameters.

Table 2

Brain regions showing abnormal nodal centrality in the heroin-dependent individuals (HDIs) compared with the normal controls (NCs).

Regions	Classification	Mean (SD)				p-value (Cohen d)
		HDIs	NCs	HDIs	NCs
HDIs<NCs
ITG.L	Association	3.50 (1.35)	5.23 (1.66)	0.15 (0.02)	0.16 (0.01)	0.0007 (1.2)	0.0005 (1.2)
MCG.L	Paralimbic	6.62 (1.99)	8.67 (1.99)	0.17 (0.02)	0.19 (0.01)	0.002 (1.0)	0.002 (1.0)
MCG.R	Paralimbic	6.63 (1.87)	8.08 (1.54)	0.18 (0.02)	0.19 (0.01)	0.004 (0.8)	0.006 (0.8)
THA.R	Subcortex	4.03 (2.34)	5.98 (2.18)	0.14 (0.03)	0.17 (0.02)	0.045 (0.9)	0.037 (0.9)
PCUN.R	Association	5.93 (1.93)	6.66 (1.99)	0.17 (0.02)	0.18 (0.02)	0.046 (0.4)	0.049 (0.4)
MFG.L	Association	5.43 (1.32)	6.58 (1.93)	0.17 (0.01)	0.18 (0.01)	—	0.048 (0.8)
HDIs>NCs
LING.L	Association	8.11 (1.50)	6.19 (1.36)	0.18 (0.01)	0.17 (0.01)	0.001 (1.3)	0.005 (1.2)
IOG.L	Association	6.51 (1.66)	5.19 (1.68)	0.17 (0.01)	0.16 (0.01)	0.012 (0.8)	0.020 (0.7)
HIP.L	Subcortex	4.50 (1.61)	3.03 (1.59)	0.15 (0.02)	0.13 (0.03)	0.014 (0.9)	0.014 (0.9)

p<0.05 (uncorrected). The symbol ‘—’ indicates no significant between-group difference. and represent the integrated nodal degree and nodal efficiency, respectively. Cohen d indicates the value of effect size. The small, medium, and large levels of the effect size are 0.2, 0.5, and 0.8, respectively, according to Cohen's definition [38]. The threshold was

Figure 3

Brain regions exhibiting abnormal integrated nodal parameters of the brain functional networks and their relationship with the duration of heroin addiction in the heroin-dependent individuals (HDIs) compared to the normal controls (NCs).

(a) Surface visualization of regions with abnormal nodal centralities using BrainNet Viewer (http://www.nitrc.org/projects/bnv/). Areas color-coded in red (blue) indicate the regions in which the values of nodal centralities corresponding to HDIs were higher (lower) than those of controls. See Table 2 for details. (b) Scatter plot of the integrated nodal parameters against the duration of heroin addiction. In the left hippocampus (HIP.L), we detected a significantly positive correlation between the integrated nodal degree and the duration of heroin addiction (p = 0.042) and a tendency toward a positive correlation between the integrated nodal efficiency and the duration of heroin addiction (p = 0.054) in HDIs. The abbreviations of regions are listed in Table S2.

Network hubs

The hub regions in the functional networks for HDIs and NCs are listed in Table S4. We found fourteen hubs in the brain functional networks of each subject group. Although the two groups had the identical number of hubs, the locations of the hubs were not completely the same. Eleven regions were shared hubs in the brain functional networks of both groups. We also found three hubs specific to HDIs, the left precuneus (PCUN.L), left postcentral gyrus (PoCG.L), and right middle frontal gyrus (MFG.R), and three hubs specific to the controls, the left middle frontal gyrus (MFG.L), right precuneus (PCUN.R), and temporal pole (TPOsup.R). We noticed that most of the shared hub regions (nine hubs) were located in the association cortices, suggesting that they had important functional roles in information transfer [14].

Functional connectivity

We utilized the NBS method to identify a single connected subnetwork with 19 regions and 19 connections, which was significantly altered in the HDIs compared to NCs (p<0.001, corrected) (Fig. 4, Table S5). We noticed that the connections in this single connected subnetwork are primarily long-distance connections linking different brain lobes. Within this subnetwork, all connections exhibited statistically significantly decreased values in HDIs (Table S5). We found that the mean connectivity value of this subnetwork correlated positively with three integrated global parameters, (r = 0.325, p = 0.069, marginally significant), (r = 0.402, p = 0.023), and (r = 0.379, p = 0.032) (Fig. 4).

Figure 4

The connected subnetwork showing decreased functional connections in the heroin-dependent individuals (HDIs) compared to the normal controls (NCs).

(a) Visualization of decreased functional connections related to heroin addiction using BrainNet Viewer (http://www.nitrc.org/projects/bnv/). The width of the line indicates the t-value of the connection comparisons between the two groups. A single, connected subnetwork containing 19 nodes and 19 connections was determined using the network-based statistic (NBS) method. (b) Scatter plot of the integrated clustering coefficient, , against the mean functional connectivity of this connected subnetwork averaged over all subjects. (c) Same as (b) but for the integrated normalized clustering coefficient, . (d) Same as (b) but for integrated small-worldness, . The abbreviations of the regions are listed in Table S2.

Correlations between network parameters and duration of heroin addiction

No significant correlations (p>0.05) were found between the integrated global parameters and the duration of heroin addiction as well as between the connections shown in Fig. 4a and the duration of heroin addiction. For the brain regions listed in Table 2, we found that the integrated degree () of the HIP.L showed a significantly positive correlation (p = 0.042), while the integrated nodal efficiency () showed a marginally significantly positive correlation (p = 0.054) with the duration of heroin addiction (Fig. 3b).

Discussion

In this study, using graph theory analysis, we constructed the functional networks, analyzed the network topological parameters, and compared the differences in these parameters of the brain functional networks between HDIs and NCs. The main findings are as follows: (1) at the global level, the heroin group showed significant decreases in the normalized clustering coefficient and in small-worldness; (2) at the nodal level, we detected significantly decreased nodal centralities primarily in regions of the cognitive control network but significant increases primarily in the HIP.L in HDIs; (3) at the connectivity level, we found a single connected subnetwork which showed significantly decreased connections in the heroin group. These findings may contribute to understanding the disrupted topological organization of whole-brain functional networks in HDIs.

The human brain is widely believed to be a complex system that requires a suitable balance between local specialization and global integration of the brain's functional activities [43]. Functional segregation and integration are two fundamental organizing principles for the human brain, a concept which is supported by the model of a small-world network characterized by a high local clustering coefficient and the shortest path length [32]. Small-world properties enable a network to maintain highly effective, specialized modular information processing as well as rapid global information transfer [44]. As has been found in previous studies of human brain functional networks [16], [17], [18], [20], [33], in this study, the whole-brain functional networks of both HDIs and NCs conserved small-worldness.

In this study, we found alterations in the global parameters of the brain functional networks of HDIs compared to NCs. Statistical analysis revealed a decreased normalized clustering coefficient in the HDIs. The normalized clustering coefficient is one of the indices that can characterize how brain networks shift to either a regular or a random network [45]. The decreased clustering coefficient in HDIs indicated that their brain functional networks may shift toward random organization. Previous studies [9], [46], [47] have suggested heroin users showed poor performance in decision making tasks compared to healthy participants. Shift toward random organization in functional network may be related with randomized decision making in HDIs. In addition, we also detected decreased small-worldness in HDIs, suggesting the topological organization in the whole-brain functional networks of HDIs was less optimal than that of the controls. Those decreases in the global network parameters in HDIs may have resulted from the decreased functional connections in a subnetwork (Fig. 4, Table S5) [19]. These connections enable spatially remote brain regions to communicate with each other and strengthen the functional integration of the brain [48]. Finding that these decreases may indicate that the brain functional integration in HDIs is disrupted.

Besides of the decreased global parameters, we also found decreased nodal centralities (integrated nodal degree and nodal efficiency) in several regions in HDIs, including the MCG.L/R, MFG.L, and PCUN.R. These regions are thought to be involved in the cognitive control network [40], [41], [42]. Previous studies have suggested drug addiction individuals exhibited deficits in neural systems associated with cognitive control [12], [49], [50], [51]. In a task-fMRI study, Kaufman et al. [49] found the dorsal cingulate cortex was less responsive during successful No-Go inhibitions in cocaine users, suggesting the drug-related dysfunction of cognitive control. Using resting-state fMRI, Ma et al. [12] found that heroin users showed reduced functional connectivity within the circuit of cognitive control, indicating the weakened strength of control in the addictive state. Recently, Liu et al. [50] studied heroin users using diffusion tensor imaging and reported that heroin users showed reduced white matter integrity in the frontal and cingulate cortex, which suggested diminished cognitive control upon craving and motivation in heroin users. Thus, our findings of decreased nodal centralities in the cognitive control regions in HDIs provided further evidence that the function of cognitive control is weakened in drug addiction [52].

Interestingly, we found the MFG and PCUN were hub regions for HDIs and NCs, but located in contralateral hemisphere, i.e., MFG.R and PCUN.L were hubs of HDIs, while MFG.L and PCUN.R were hubs of NCs. This finding may reflect the existing compensatory mechanism or neuroadaptation in the addiction brain [53], [54], [55]. For example, Jager et al. [54] found that the adolescent cannabis users showed excessive activity in the prefrontal regions during a novel task, suggesting functional compensation. Also Kanayama et al. [55] reported that the cannabis users might call upon additional brain regions not typically used for spatial working memory (such as regions in the basal ganglia) to compensate for the deficits in spatial working memory. In the present study, the brain functional networks of HDIs may need to enhance the function of MFG.R and PCUN.L to compensate for impaired function of MFG.L and PCUN.R due to their decreased nodal centralities in HDIs.

We also found that the HDIs showed increased nodal centralities in the HIP.L compared to NCs, a finding which was consistent with several previous studies [11], [52], [56]. Using graph theory analysis, Liu et al. [11] suggested that the hippocampus had a higher nodal degree in the brain of chronic heroin users. Ma et al. [56] found that heroin users showed increased functional connectivity in the hippocampus compared to controls. Baler and Volkow [52] demonstrated that the memory/learning circuit related to drug addiction is primarily located in the amygdala and hippocampus. In fact, the hippocampus is the main brain region involved in memory and learning [57], and is thought to strengthen the learning of drug-related cues which leads to drug-seeking behaviors [58]. Therefore, an increase in nodal centrality in the HIP.L may excite the expectation of the drug in HDIs.

Moreover, the integrated nodal degree and the nodal efficiency were positively correlated with the duration of heroin addiction, but only in the HIP.L. This indicated that longer the heroin use, the higher the nodal centrality of the HIP.L. Thus, the disrupted topology properties in the hippocampus may indicate that the pattern of relapse to drug-seeking behaviors that is commonly seen in HDIs is driven by abnormal memory processing.

Notably, the disrupted nodal topology can also be interpreted from the perspective of physiological aspect. Previous studies [59], [60], [61], [62], [63] have suggested that the effects of opioid drugs on the brain might depend on the opioid receptor density. The frontal cortex and cingulate cortex (anterior and middle) have the high opiate receptor-binding potentials [64] and have been reported to be commonly affected by different opoioid drugs, such as cocaine [59], nicotine [61], morphine [62], [65] and remifentanil [63], [66]. Thus, the current findings of decreased nodal centralities in the left middle frontal cortex and bilateral middle cingulate cortex in HDIs may reflect an outcome of disrupted opioidergic modulation. Actually, we cannot attribute all altered nodal centralities to opioid receptor. In current study, we also found heroin affected nodal centrality in the left hippocampus which has not been reported to include high opioid receptor. This is in line with several task and resting-state pharmacological fMRI studies [62], [65], which reported that morphine affected functional topography of hippocampus in healthy volunteers. We noticed that no change of nodal centrality has been detected in this study in at least several areas, the insula, thalamus, amygdala, and putamen, though these regions are more susceptible to the high opioid receptor [64]. These suggest that the findings of abnormal nodal centralities in HDIs partly reflect the opioid receptor distribution.

Limitations

Several limitations need to be addressed. First, due to the cross-sectional nature of this study, we can only infer that the network properties of the brain functional networks of heroin addicts are disrupted. We are not able to determine the precise relationship between heroin abuse and abnormalities of the network parameters. Second, the nodal centrality results could not survive when we adopted multiple comparisons (FDR and FEW corrections), meaning this should be considered as an exploratory analysis. To increase the statistical power, the findings need replication with a larger sample of subjects or a limited number of selected ROIs. Third, the HDIs received methadone treatment at the time of the fMRI study which might affect the brain spontaneous activity [67], [68] and the topological properties of functional network in the present study. Therefore, in the future, we should design a more rigorous experiment to exclude the effect of methadone. Fourth, we could not completely eliminate the effects of physiologic noise due to the low sampling rate (TR = 2 s), which can cause respiratory and cardiac fluctuations to impact the fMRI time series, even though a 0.01–0.08 Hz band-pass filter was used to reduce this effect. Finally, we only estimated the relations between the network properties and the duration of heroin addiction. Whether brain functional network properties are related with other clinical variables including decision-making behavior, impulsivity and consequences on daily life should be explored in the future.

In summary, we investigated the whole-brain functional network in HDIs using resting-state fMRI and a graph theory method. We found that compared to the normal controls, the whole-brain functional networks in HDIs may shift toward random organization, as indicated by a lower normalized clustering coefficient and lessened small-worldness. We also found that the nodal properties were disrupted, especially in regions of cognitive control network in HDIs. Our study indicated disruptions in the whole-brain functional networks of HDIs, findings which may be helpful for better understanding the mechanisms underlying heroin addiction.

Network-Based-Statistic Analysis (NBS).

(DOC)

Click here for additional data file.

The detailed clinical description for each heroin-dependent individual (HDIs).

(DOC)

Click here for additional data file.

The names and the corresponding abbreviations of the regions of interest (ROIs).

(DOC)

Click here for additional data file.

The mathematical definitions and descriptions of global metrics and nodal metrics.

(DOC)

Click here for additional data file.

Hub regions of the brain functional networks detected in the heroin-dependent individuals (HDIs) and normal controls (NCs).

(DOC)

Click here for additional data file.

Decreased functional connections in the heroin-dependent individuals (HDIs) as compared to the normal controls (NCs).

(DOC)

Click here for additional data file.