Relapse risk revealed by degree centrality and cluster analysis in heroin addicts undergoing methadone maintenance treatment

Abstract

Background

Based on hubs of neural circuits associated with addiction and their degree centrality (DC), this study aimed to construct the addiction-related brain networks for patients diagnosed with heroin dependence undertaking stable methadone maintenance treatment (MMT) and further prospectively identify the ones at high risk for relapse with cluster analysis.

Methods

Sixty-two male MMT patients and 30 matched healthy controls (HC) underwent brain resting-state functional MRI data acquisition. The patients received 26-month follow-up for the monthly illegal-drug-use information. Ten addiction-related hubs were chosen to construct a user-defined network for the patients. Then the networks were discriminated with K-means-clustering-algorithm into different groups and followed by comparative analysis to the groups and HC. Regression analysis was used to investigate the brain regions significantly contributed to relapse.

Results

Sixty MMT patients were classified into two groups according to their brain-network patterns calculated by the best clustering-number-K. The two groups had no difference in the demographic, psychological indicators and clinical information except relapse rate and total heroin consumption. The group with high-relapse had a wider range of DC changes in the cortical−striatal−thalamic circuit relative to HC and a reduced DC in the mesocorticolimbic circuit relative to the low-relapse group. DC activity in NAc, vACC, hippocampus and amygdala were closely related with relapse.

Conclusion

MMT patients can be identified and classified into two subgroups with significantly different relapse rates by defining distinct brain-network patterns even if we are blind to their relapse outcomes in advance. This may provide a new strategy to optimize MMT.

Introduction

Heroin addiction is a chronic disease characterized by compulsive drug seeking and use. In clinical practice, the maintenance of pharmacotherapy is a strong determinant of the outcome, but adherence to treatment remains problematic due to that relapse is common (Volkow, Jones, Einstein, & Wargo, 2019). A deep understanding of the mechanism of addiction is important and may provide an opportunity of seeking new pharmacotherapeutic targets.

By means of neurobiological imaging, the procedure from hedonic drug use to compulsive drug seeking has been well described and elaborated with a four-circuit addiction model composed of reward circuit, motivation/drive circuit, memory/learning circuit, and control circuit (Volkow, Fowler, & Wang, 2003). The reward circuit mainly comprises the ventral tegmental area (VTA) and nucleus accumbens (NAc); the motivation/drive circuit comprises the orbitofrontal cortex (OFC), thalamus, caudate and putamen; the memory/learning circuit comprises hippocampus, amygdala; the control circuit comprises anterior cingulate cortex (ACC), inferior prefrontal cortex (IPC), dorsolateral and prefrontal cortex (dlPFC). It was proved that between the four circuits existed functional unbalance interplay, such as overactivating of reward circuit and motivation/drive circuit and lack of interaction of control circuit with reward circuit and motivation/drive circuit (Volkow et al., 2003; Volkow, Wang, Tomasi, & Baler, 2013). Recently, studies supplemented that insular cortex was an essential structure for addiction behavioral maintenance representing the interoceptive effects of drug taking and making this information available to conscious awareness, memory and executive functions (Naqvi, Gaznick, Tranel, & Bechara, 2014). Overall, the brain addiction-related topological organization model has been proposed and it, as a network comprising a number of modules made of nests sub-modules (circuits), is mediated by inter-modular connections of hubs.

Although the four-circuit model provides a relatively comprehensive neural framework for understanding the main neurobiological processes through which biological and sociocultural factors contribute to resilience against or vulnerability to drug use, it cannot completely explain the remarkable differences in the outcomes among addicts with pharmacotherapy. Methadone maintenance treatment (MMT), the most effective and accepted therapy for heroin dependence, for example, has situations reported like illicit drug use or retention rates altering widely among clinics (Gauthier, Eibl, & Marsh, 2018; Zhou & Zhuang, 2014). Clinical neuropsychology studies demonstrated that multiple factors, e.g. age, gender, genetic variation, methadone dose, social experience, stress, trauma and environmental factors, may individually or combinedly contribute to the gap of MMT outcomes and further suggested that those factors could be used as predictors for MMT outcomes to a certain extent (Andersen, 2019; Coller, Barratt, Dahlen, Loennechen, & Somogyi, 2006; Darker, Ho, Kelly, Whiston, & Barry, 2016; Jaremko, Sterling, & Van Bockstaele, 2015; Lin, Hung, Peng, Chao, & Lee, 2015; Lister, Brown, Greenwald, & Ledgerwood, 2019; Luo et al., 2017; Proctor et al., 2015; Zhou, Li, Wei, Li, & Zhuang, 2017). However, few of the predictors or strategies was applied into clinical practice because of diversity, inconsistency and subjectivity of the studies' results.

Considering that the fundamental neural activity pattern underlies the relapse risk gap among MMT patients, scientists applied neuroimaging techniques to search the neuronal substrates associated with an eventual outcome to MMT. Li et al., studied the relationship between the percentage of positive urine drug screens and the functional connectivity (FC) strength between salience network, default mode network (DMN) and executive control network in MMT patients, and suggested the disrupted coupling between the salience network and DMN and between the left executive control network and DMN was associated with relapse behavior (Li et al., 2018). With independent component analysis, Li and her colleagues compared the FC of DMN between the relapsers and abstainers during MMT, revealing the potential predictive value of DMN concerning heroin relapse under MMT (Li et al., 2015b). Some other neuroimaging indexes were also addressed for the association with relapse to heroin use in MMT, such as fractional anisotropy and axial diffusivity in the posterior limb of the internal capsule (Li et al., 2016c), regional homogeneity value in caudate (Chang et al., 2016), and cue-induced activation in NAc/subcallosal cortex (Li et al., 2015a).

Although several relapse-related brain regions and neuroimaging indices were identified by studying the neuronal activity characteristics in the neuroimaging researches mentioned above, the practicality of these results was restricted from their retrospective evaluation design which means that participants' relapse outcomes were known in advance. Thus, a well-designed prospective cohort study using neuroimaging techniques at the condition of being completely blind to relapse outcomes was necessarily needed. Furthermore, an addiction-related brain network should be a suitable index of analysis for predicting patients' response to MMT, because the pathological alteration led by heroin use is wide at whole-brain level.

According to prior researches, the topological patterns in addiction-related networks were made of a set of submodules-neural circuits, consisting of some key nodes densely connected to each other and inter-modular connector hubs mentioned above. Degree centrality (DC), as topological features of graphs meaning the number of functional relationships of a given voxel with the rest of the network, can describe the addiction-related graphical models of the brain better (Buckner et al., 2009; Guo et al., 2016; Luo et al., 2017; Wang, Jiao, Zhang, & Lin, 2017; Zhou et al., 2014; Zuo et al., 2012). And it was conceived as a stable and sensitive parameter in test−retest in the face of BOLD signal noise that was difficult to overcome by other methods (Zuo et al., 2012). So it has been widely applied to examine node characteristics of intrinsic connectivity networks and the topology of brain hub changes with the disease at the global voxel level, for example, to investigate the association of alcohol dependence with changed intrinsic functional hubs and to reveal the attribution of several neural systems to the impaired motor behavior (Luo et al., 2017). Even DC was suggested as an index for attentional impulsivity in a recent study on codeine-containing cough syrups dependence (Hua et al., 2018).

In this current study, based on the key nodes and inter-circuits brain hubs disclosed by prior studies on addiction and those using DC value, we constituted the addiction-related brain network of MMT patients and identified the relapse-related brain network organization pattern with cluster analysis. Our aim was to prospectively identify the poor MMT responders with neuroimaging data-driven method and reveal the brain regions contributing to relapse behavior. It was to be expected that the prospective identification of poor responders to MMT with cluster analysis based on neuroimaging data would be useful in modifying intervention strategies for improving the outcomes of MMT patients in future.

Methods

Participants

This study was approved by the Institution Board of the Fourth Military Medical University and conducted in accordance with the Declaration of Helsinki. Sixty-two patients diagnosed with heroin dependence receiving MMT (MMT patients) (male; mean age, 35.8 years; age range, 22–53 years) and 30 matched male healthy controls (HC) (mean age, 34.8 years; age range, 19–48 years) were recruited in this study. All participants were fully informed about the details of the experiment and signed the written consents for their involvement. Recruited from the outpatient of Xi'an Methadone Substitution Treatment Center, all the MMT patients met the DSM-IV criteria for heroin dependence, received MMT for no less than 6 months, and had at least 1 month of stable daily methadone dose before entering the study. Other inclusion criteria included: (1) right handed; (2) aged 18–55 years old; (3) heroin use history of more than 12 months. They were excluded if they had current or past psychiatric medical illness other than heroin dependence. HC were recruited from local community by advertising and they were free of DSM-IV-TR Axis I disorders and had no history of drug use. The exclusion criteria for all participants included: (1) current or past major medical illnesses or current use of prescription medications; (2) dependence other than psychoactive substances (except nicotine); (3) head trauma or neurological illness history; (4) MRI examination contraindication.

Evaluation measures

The demographics information (age, educational level), heroin use history (lifetime use and use dosage per day) and methadone treatment information (duration and methadone daily dosage) of each participant were collected during the enrollment interview. Before the MRI scan, all participants were asked to complete the questionnaire survey under psychiatrist guidance, including Beck depression inventory (BDI), Hamilton Anxiety Scale (HAMA). Furthermore, subjective heroin craving was assessed using a 0–10 visual analog scale, asking, “To what extent do you feel the urge to use heroin?” (0 indicating the least craving and 10 the strongest).

According to the reports from MMT fix-clinics in different cities of China, the MMT retention rate varied greatly (Cao et al., 2014; Jiang et al., 2014; Zhang et al., 2013). It was notable that studies on retention in MMT among opioid-dependent patients in Xi'an revealed the dropout and relapse mostly occurred between 12 and 60 months after MMT initiation (Wei et al., 2013; Zhou & Zhuang, 2014). Thus, the follow-up duration in this study was set as 26 months.

After MRI scan, the MMT participants were followed for 26 months. During each appointment, the participant underwent a monthly structured interview assessing illicit drug use and a urine drug test. As opiates (mainly heroin, 38.1%) and synthetic drugs (mainly methamphetamine, 60.5%) were the main drugs of abuse in China (China Anti-Drug Network, annual report in 2017), relapse was defined as any use of heroin or methamphetamine identified by positive urine drug test (Morphine/Methamphetamine Diagnostic Kit, Guangzhou Jianlun Biological Technology Co., Ltd). In addition, the participants who missed the appointments without any contact with the study coordinator or who were loss of follow-up were considered to be relapsed. Finally, the accumulative total amount of positive urine tests and a missed interview was collected to calculate MMT relapse rate. Our coordinator contacted MMT participants 3 days before each appointment and encouraged them to continue the MMT program. There were 3 MMT patients lost in the follow-up period, two patients lost to follow up after MRI scan and one lost after the first clinical structural interview. There were 451 relapse counts totally, among them 216 were due to missing appointments and 139 due to loss to follow-up.

Brain imaging methodology and data analysis

Magnetic resonance imaging acquisition

MRI data were acquired with a 3.0T GE-Signa HDxt MRI scanner using an eight-channel head coil (GE Healthcare, Milwaukee, USA) in the Department of Radiology, Tangdu Hospital, the Fourth Military Medical University. Alcohol, tea, caffeine, any drug or medicine were prohibited in the 12 h before MRI scanning. All the MRI scanning were conducted 4 h after the MMT patients took their daily dosage (methadone peak plasma level) to avoid drug withdrawal symptoms (Wolff, Hay, Raistrick, & Calvert, 1993).

Lying supine with the head fixed by a belt and foam pads, participants were instructed to keep their heads still, close their eyes and not to think anything specific. A routine structure MRI scanning was conducted to exclude gross cerebral pathology. After that, the resting-state fMRI (rs-fMRI) data were obtained using a gradient-echo planar imaging (GRE-EPI) pulse sequence, settings as: 30 axial slices, TR = 2000 ms, TE = 30 ms, flip angle = 90°, FOV = 256 mm × 256 mm, slice thickness = 5.0 mm, skip = 0 mm, matrix = 64 × 64. This session lasted for 5 min and 10 s, including a 10-s dummy scan period at the beginning of the scan. For each participant, 150 echo-planar volumes were collected during the rs-fMRI scan. And then a high-resolution structure session was performed with contiguous slices to cover the whole brain in a steady state using an axial fast spoiled gradient recalled echo (3D-FSGPR) for spatial normalization of the data sets to a standard atlas, scanning parameters: TR = 7.8 ms, TE = 3.0 ms, TI = 450 ms, FOV = 256 mm × 192 mm, slice thickness = 1.0 mm, skip = 0 mm, matrix = 256 × 256. After scanning, all participants reported they were awake during all scans.

Imaging data processing

Data preprocessing was conducted with SPM 8 software (http://www.fil.ion.ucl.ac.uk/spm) and DPABI (http://rfmri.org/dpabi). Prior to data preprocessing, the first five volumes were discarded to improve signal stabilization. Then, the remaining rs-fMRI data were corrected for slice timing and realigned for motion correction. The standard Montreal Neurological Institute (MNI) template was used for spatial normalization with a resampling voxel size of 3 mm × 3 mm × 3 mm. Participants with head motion more than 2 mm of translation and 2° of rotation in any of the x, y and z axes or mean frame-displacement more than 0.2 were excluded. Next, nuisance signals (including Friston 24-head motion parameters, the white matter and cerebrospinal fluid signals) were extracted and regressed out from the time series of every voxel to reduce the effects of nonneuronal signals. After a linear trend of the time course was removed, a band-pass filter (0.01–0.1 Hz) was applied. Then, scrubbing signal spikes based on the method of Power and colleagues was conducted (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012). Specifically, a threshold of Frame-Wise Displacement = 0.5 mm was used and the bad volume and the volumes including 2 before and 1 after were removed. Two MMT participants and one HC were excluded from the next analysis because of excessive head motion or mean frame displacement.

Degree centrality (DC) calculation

Weighted DC measures were calculated using the “REST-DC” toolkit in the REST V1.8 packages (Zuo et al., 2012) as previously described (Di Martino et al., 2013; Liu et al., 2015). Briefly, to obtain each participant's graph, Pearson correlation coefficients were computed between the time series of all pairs of brain voxels within gray matter mask. Each voxel represented a node in the graph, and each significant functional connection (i.e. Pearson correlation) between any pair of voxels were the edges. As a result, a N × N undirected adjacency matrix (N was the count of all voxels within the gray matter mask) was obtained to construct the whole-brain FC matrix for each participant, each element of the matrix represents the Pearson correlation coefficient between pairs of voxels. To eliminate the possible spurious connectivity (Li et al., 2016b), the adjacency matrix was further thresholded with each Pearson correlation coefficient at r > 0.25. The weighted DC of a voxel was calculated as the sum of the connections' strength between a given brain voxel and all other voxels.

The formula is: $DC_i = \mathop {\sum {a_{ij}} }\limits_{j = 1}^N$.

Where the DC _i represents the weighted DC value for given voxel i, and a _ij means the Pearson correlation coefficient between voxel i and voxel j.

Furthermore, standardized weighted DC (DC Z-score, zDC) maps were acquired by subtracting the mean value, and then divided by the standard deviation within the whole gray matter mask (Takeuchi et al., 2015; Zuo et al., 2012) according to the Z-score standardization formula below:

$$Z_i = \displaystyle{{{\rm D}{\rm C}_i-\mu } \over \sigma }, \;\,1 \le i \le N.$$

Where Z _i represents the Z-scored DC value of voxel i, μ and σ are the mean and standard deviation of the DC measures across all N voxels. Finally, the resulting DC maps were spatially smoothed with a 4-mm FWHM Gaussian kernel.

Addiction-related network constitution

Based on the current addiction neurobiological theoretical framework, 10 pairs of brain nodes and hubs were chosen from “four circuits model” and addiction theory of insula and set as regions of interest (ROIs) to constitute brain addiction-related network for MMT patients. These regions included the bilateral NAc, amygdala, ACC, caudate, OFC, hippocampus, insular, putamen, thalamus, and DLPFC. The ROIs were determined with FSL https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL (Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012; Smith et al., 2004; Woolrich et al., 2009) according to Harvard–Oxford cortical and subcortical structural atlases based on MNI 152 1 mm template. The threshold probability was set as >50%. The DLPFC ROI was obtained by merging BA9 and BA46 of WFU PickAtlas http://fmri.wfubmc.edu/software/PickAtlas (Lancaster et al., 2000; Maldjian, Laurienti, Kraft, & Burdette, 2003). All ROI masks were resampled into 3 mm functional space. For characterizing the topological properties of brain networks, the DC value was extracted from the 10 pairs of ROIs of the 60 MMT patients after regressed out total methadone consumption, and used to make a 60 × 20 matrix M, representing the addiction-related network for subjects.

Cluster analysis and cluster validation and stability test

The clustering analysis was used to classify MMT patients' matrix M into subgroups with R (https://www.R-project.org/) and several R packages (factoextra, NbClust). The procedure included six primary steps: (1) assessing clustering tendency of the matrix M with the Hopkins' statistic and a visual approach; (2) identifying the optimal number K of clustering with Euclidean distance of ward D2 aggregation algorithm; (3) grouping the MMT participants according to K-means clustering algorithm; (4) clustering efficiency test with Silhouette analysis (Sun et al., 2015); (5) clustering stability test with Fowlkes–Mallows (FM) index; (6) grouping stability test with the method provided by Hening C (online Supplementary-Materials).

Statistical analysis

To investigate the DC features coupling with subgroups, the differences in DC between MMT subgroups and HC were conducted with two-sample t test on the whole brain level (TFCE corrected p < 0.05 and cluster size K > 10). The difference in relapse rate between subgroups was calculated via Pearson's χ² test with Yates' continuity correction. The differences in age, education level, smoking habit, heroin use information, methadone use history and behavioral score between MMT subgroups were calculated using two-sample t test or analysis of variance (>2 subgroups). The results were considered as significant if p value were less than 0.05.

Regression analysis

To explore which brain region contributed to relapse significantly, a sphere ROIs with 3 mm radius was made according to the coordination of peak voxel of brain region with significantly different inter-subgroup DC value. Robust Poisson regression (Tsou, 2006) was run to predict the relapse rate based on the z-scored DC signal strength for each ROI in 60 MMT participants, in which relapse count was set as dependent variable and only one brain DC value as independent variable in the model at each time. The statistical threshold was adjusted to 0.05/n after Bonferroni correction (n = the number of brain regions with significantly different inter-subgroup DC values).

Results

Clustering analysis

The results of clustering are shown as an ordered dissimilarity image (ODI) of 60 × 20 matrix M (Fig. 1). According to the index of clustering tendency, Hopkins statistic value for M was 0.38, far less than 0.5, indicating the clusterability of the matrix M. Besides, it was also implied by the ODI that the matrix M was of significant clusterability. The optimal clustering number K was determined as 2 according to Wards' algorithm and the calculations of 7 clustering indices out of 26. Clustering efficiency and stability tests also suggested the 2-cluster solution was the most stable and efficient one (online Supplementary-Materials). Finally, the 60 MMT participants were classified into two subgroups, 29 MMT patients (48.3%) in subgroup-1 and 31 (51.7%) in subgroup-2.

Fig. 1. Ordered dissimilarity image of matrix M. The color level is proportional to the value of the dissimilarity between observations. Objects belonging to the same cluster are displayed in consecutive order. The dissimilarity matrix image confirmed that there is a cluster structure in the HD participants' data set. Two main subgroups (subgroup1 and subgroup2) were identified. Red: High similarity, Blue: Low similarity.

Differences in demographic data, clinical information between MMT subgroups

There was a significant difference in the relapse rate between subgroup-1 and subgroup-2 [χ² = 22.578, df = 1, p = 2.018 × 10⁻⁰⁶; ϕ = 0.13, 95% confidence interval (CI)  = (0.07–0.017)], subgroup-1 higher than subgroup-2 (Fig. 2). Besides, subgroup-1 was significantly higher than subgroup-2 in the total heroin consumption [p = 0.047; Cohen's d = 0.58, 95% CI = (0.06–1.11)]. No significant differences existed in age, education level, methadone use information, BDI, HAMA, craving score, smoking duration, and daily cigarette number between the two subgroups (p > 0.05) (Table 1).

Fig. 2. Violin plot of relapse rate of HD subgroups. HD subgroup1 had a significant higher relapse rate than subgroup2 (p < 0.05).

Table 1. Demographic and clinical information of study subjects and subgroups

Note: HC, healthy control; HD, patients diagnosed with heroin dependence; S1, subgroup-1; S2, subgroup-2; NA, not applicable; s.d., standard deviation; *, statistically significant; HD-HC, HD v. HC, and so on.

Differences in DC value between MMT subgroups and HC group

As a result of circular analysis, the areas with significant differences in DC value between MMT subgroups and HC group reported were almost the ROIs used for cluster analysis. Compared with HC, MMT subgroup-1 had more and wider distribution of brain areas with changed DC, including the bilateral caudate, putamen, thalamus, left amygdala, hippocampus, vACC, and OFC with decreased DC, the bilateral paracentral gyrus and right postcentral gyrus with increased DC. On the contrary, subgroup-2 showed a less brain regions with altered DC, the bilateral calcarine and lingual gyrus with decreased DC value. When compared to subgroup-2, subgroup-1 had 17 brain regions with decreased DC, including the bilateral amygdala, caudate, OFC, thalamus, hippocampus, NAc, putamen, vACC, and left insular (TFCE corrected p < 0.05 and cluster size K > 10) (Fig. 3).

Fig. 3. DC differences between HD and control and between subgroup 1 and subgroup 2 (TFCE corrected p < 0.05 and cluster size K > 10). Subgroup1 had a widespread altered density of connectivity involving in cortical–striatal–thalamic circuit at the whole-brain level, whereas a limited alteration in the subgroup2. As a result of circular analysis, the areas with significant differences in DC value between MMT subgroups and HC group reported were almost the ROIs used for cluster analysis (the upper, middle, and lower rows, respectively).

Regression analysis

The DC value of the bilateral vACC and NAc, hippocampus and right amygdala of MMT patients had a significant negative relationship with relapse rate. For every unit increased in the DC value of left vACC, right vACC, left NAc, right NAc, left hippocampus, right hippocampus and right amygdala, relapse rate decreased by 53.57, 58.11, 55.47, 43.41, 42.69, 50.89 and 41.4%, respectively (p = 6 × 10⁻⁰⁶, Exp(B) = 0.4643; p = 0, Exp(B) = 0.4189; p = 1 × 10⁻⁰⁶, Exp(B) = 0.4453; p = 3 × 10⁻⁰⁶, Exp(B) = 0.4659; p = 0.001 Exp(B) = 0.5731; p = 5.3 × 10⁻⁰⁵ Exp(B) = 0.4911; p = 2 × 10⁻⁰⁴, Exp(B) = 0.526) (Fig. 4).

Fig. 4. Estimates for signal of six brain ROIs in the models showed by Robust Poisson regression coefficient plot with 95% confidential intervals.

Discussion

This present study constructed the addiction-related brain network for MMT patients and at the condition of being completely blind to relapse rate in advance, we, by clustering the data set of the network, identified two MMT subgroups with different relapse-level defined by distinct patterns of addiction-related brain network abnormalities. The subsequent statistical analysis of clinical information confirmed that no differences were found in the demographic characteristics, methadone use or BDI and HAMA score between the two subgroups, and that subgroup-1 was the poor responder to MMT because it had a higher relapse rate. Here subgroup-1 was called as high relapse group (HRG) and subgroup-2 low relapse group (LRG). We think the prospective identification of HRG of MMT patients is the highlight of this study, because it may provide an opportunity for improved treatment tailoring and/or based on it, healthcare providers can adjust more aggressive adjunct therapies for the high-risk populations. In the following discussion part, we do not discuss the areas with significant DC value between the MMT subgroups and HC group to avoid emphasizing circular analysis.

Our study demonstrated that HRG had a widespread altered density of connectivity involving a cortical-striatal-thalamic circuit at the whole-brain level, whereas a limited one in the LRG. This finding implied that the brain network pattern of HRG identified is characterized by alteration of the importance of hubs extensively. Two factors may account for such global deficits.

The first one might be related to chronic heroin use. Heroin exposure cause brain structural injuries, including cell apoptosis (Tramullas, Martinez-Cue, & Hurle, 2008), mitochondrial dysfunction (Feng et al., 2013), synaptic defects (Garcia-Fuster et al., 2008), spongiform leukoencephalopathy (Pirompanich & Chankrachang, 2015), disturbance of neurogenesis development (Bayer et al., 2015), neuronal loss (Feng et al., 2013). In line with these findings, gray matter nucleus volume loss (Muller et al., 2015), abnormal cortical thickness (Li et al., 2014) and white matter structural connectivity damage (Liu et al., 2008; Wollman et al., 2015) were discovered in morphological studies. Correspondingly, disrupted activities of neuronal clusters, brain regions, neuronal circuits and brain networks in heroin hijacked brain were revealed in functional imaging studies (Chang et al., 2016; Jiang et al., 2011; Li et al., 2016a; Wang et al., 2013; Zhang et al., 2017). Actually, the pathological alterations resulted from heroin use is not equally distributed over the brain, but preferentially affect the hub regions (Volkow & Morales, 2015; Volkow, Wang, Fowler, & Tomasi, 2012; Volkow, Wang, Fowler, Tomasi, & Telang, 2011). Thus, the brain hubs with prominent reduced DC in HRG may reflect the neuropathological features associated with chronic heroin use.

From the view of neurodevelopment, the large-scale changed DC in hubs of high relapse risk network pattern could be interpreted as the consequence of complex interactions between genetic, epigenetic and environmental factors. Genetic factor accounts for approximately half of the risk for addiction (Volkow & Li, 2005), and a series of gene variants are associated with addiction vulnerability by impacting on synaptic plasticity (Oliver et al., 2018; Randesi et al., 2018), receptor binding affinity (Shi et al., 2002), dopamine receptor density (Noble, 2000), neuron cell sensitivity to the substance of abuse (Alia-Klein et al., 2011), cortical neurotransmission and neurogenesis (Crews & Vetreno, 2011), and gray matter volume shrinkage in chronic substance use (Gitik et al., 2016). In addition, the epigenetic factors also contribute to the risk for addiction. People who experienced early life stress are at significantly high risk for the development of addiction caused by stress-induced neuro-plastic changes (Nestler, 2014). Other factors associated with addiction vulnerability include neurodevelopment imbalance between circuits supporting reward-seeking behaviors and self-control during adolescence (Giedd, 2008), early exposure to drugs of abuse (Nestler, 2014) and psychological traits of impulsivity, excitement seeking and stress reactivity (Bickel, Jarmolowicz, Mueller, Gatchalian, & McClure, 2012; Holmes, Hollinshead, Roffman, Smoller, & Buckner, 2016; Jasinska, Stein, Kaiser, Naumer, & Yalachkov, 2014). Consistent with the results, the longer heroin use duration in HRG indicated a younger heroin use age relative to LRG, which might imply the influence of gene and environment (Su et al., 2015).

The areas with significantly reduced DC in the primary and secondary visual cortex in the HRG might reflect the target of action of methadone. The previous study has demonstrated that methadone can depress visual function by acting on the visual parts, including midbrain and thalamic visual nuclei (Rothenberg, Peck, Schottenfeld, Betley, & Altman, 1979). Then, the relationship between the disrupted neural activity of occipital and long-term methadone use is confirmed by electroencephalograph techniques (McGlone et al., 2008; Wang, Kydd, & Russell, 2016). Although we did not observe the similar DC change of visual cortex in the LRG, we still think it is worth exploring the possible methadone-associated alteration of visual cortex function with imaging technology in the future. The increased DC value of paracentral and postcentral lobe in the HRG is inconsistent with those of prior opiate study, such as decreased glucose metabolism, reduced response to heroin-related cue and GO/NOGO task in paracentral and/or postcentral lobe (London et al., 1990; Mei, Zhang, & Xiao, 2010; Ye et al., 2018). It might be a reflection of methadone playing its therapeutic role in addicted brain by increasing regional blood cerebral flow (Danos et al., 1998; Jiang et al., 2011). In an early study performed in heroin addicts, postcentral gyrus showed a reduced response to drug-related cue along with a decreased craving after buprenorphine dose (Mei et al., 2010). Another study on positive therapeutic effect of acupuncture on alcohol use disorder revealed that several brain regions including postcentral gyrus were activated (Yang et al., 2017). Thereupon, we speculate the increased DC could be a reflection that treatment works. However, this phenomenon of increased DC of paracentral and postcentral region was not observed in LRG, which still calls for further investigations to understand the neuropathology underlying it.

Our regression analysis revealed that several hubs of addiction-related circuits, such as NAc, vACC, hippocampus and amygdala, were closely related to relapse. It raises the possibility of the potential therapeutic targets for addiction treatment.

NAc as a crucial hub in the network of reward, motivation and craving has rich and extensive interconnections with the VTA, amygdala, thalamus, prefrontal cortex, etc. Because of the key role in the process from voluntary drug use to compulsive drug-seeking behavior, NAc has been identified as a predictor for abstinence period following cession and a common target of pharmacological and cognitive-based interventions (Konova, Moeller, & Goldstein, 2013; Luigjes et al., 2012; Owens et al., 2017). In an earlier clinical study, a prominent improvement of 5-year retention rate was observed in former heroin addicts with bilateral NAc ablation stereotactic surgery (Li et al., 2013). Furthermore, animal studies revealed that the possible way to attenuate drug reinstatement by manipulating NAc is local activation/or activation of GABAergic interneurons in the medial prefrontal cortex via antidromic stimulation of cortico-accumbal afferents (Vassoler et al., 2013).

Hippocampus is involved in the formation of associations between drug use and special events and the regulation of reinstatement of drug-reinforced response (Francis, Chaudhury, & Lobo, 2014; Koob & Volkow, 2010; Pascoli et al., 2014). Exposure to drug-associated contextual cue contributes significantly to relapse, which is mediated by drug-induced synaptic plasticity change in the ventral hippocampus (Alvandi, Bourmpoula, Homberg, & Fathollahi, 2017; Borjkhani, Bahrami, & Janahmadi, 2018). Therefore, there has been diverse and efficacy interventions targeting the hippocampus to prevent recurrent relapse in drug-dependent animal models. Wright, V. L. successfully used alpha7 nicotinic acetylcholine receptor antagonist to block the rats' reinstatement of morphine-conditioned place preference which was mediated by the retrieval of associative drug memories in ventral hippocampus (Wright et al., 2018). Another experiment proved the reinstatement induced by stress or combination of stress and primer dose of morphine was attenuated by antagonist administration targeting D1- and D2-like receptor in CA1 region of the hippocampus, addressing the role of hippocampal structures in treatment (Nazari-Serenjeh, Rezaee, Zarrabian, & Haghparast, 2018). In addition, several studies implicated that the drug-induced increase in synaptic plasticity in hippocampus projections to NAc shell played role in the contextual-mediated associations with drug-taking, and manipulation on hippocampus influenced drug reinstatement and drug-seeking behavior (Britt et al., 2012; Pascoli et al., 2014).

vACC is associated with relapse vulnerability and therapeutic outcomes (Forster, Dickey, & Forman, 2018). A study demonstrated that the abstinent time-related FC pattern of vACC and frontopolar cortex in a group of recently abstinent stimulant use disorder patients had a shift from strengthened FC first to a greatly reduced FC afterwards, implicating the patients' relapse risk (Camchong et al., 2014). And it was found in cocaine dependents that their increasing impulsivity and higher relapse rate was associated with the aberrant FC between ventral caudate and subgenual ACC (Contreras-Rodriguez et al., 2015). It was confirmed that the drug-induced neuroadaptation causing an imbalance between emotional processing and cognitive processing did also exist in the vACC (Shapira-Lichter et al., 2018; Wang et al., 2010). In view of this, it is reasonable to presume that vACC could be a potential medical interference target for addiction treatment.

Compelling evidence demonstrated that the amygdala, specifically the basolateral nucleus appears to play a critical role in relapse related to exposure to reminders of drug consumption, drug withdraw, fear or stress (Goode & Maren, 2019; Koob, 2009). The neuronal substrates underlying it is the afferents from the VTA to the basolateral nucleus that increase its activity and the reciprocal connections with the NAc involved with the information of seeking behavior (Ambroggi, Ishikawa, Fields, & Nicola, 2008; Ford, Mark, & Williams, 2006; Grace & Rosenkranz, 2002). As the potential mediator of the conditioned cue reinstatement, the amygdala has been chosen as the target of medication in the context of relapse. Disruption of neural activity in the basolateral nucleus or central nucleus of the amygdala in a rat model of relapse after chronic cocaine or ethanol self-administration blocks the initial acquisition of cue conditioning and the subsequent relapse, as well as withdrawal-induced drinking (Bale & Vale, 2004; See, Fuchs, Ledford, & McLaughlin, 2003). These findings are also strong evidence supporting the opinion about the negative correlation between amygdala and relapse rate proposed by our study.

Although we successfully divided the MMT patients into two subgroups with statistically different relapse rate according to the distinct addiction-related network patterns, the classification could not be as good enough as we had expected. Several MMT patients might be arranged into the wrong subgroup. We think the limitations of this study may contribute to this.

The first limitation should be originated from the selection strategy of ROIs for constructing Matrix M representing addiction-related network. We were mainly based on classical addiction neurobiological theories to determine the ROIs. Despite it is relatively simpler, the ROIs we selected, perhaps, could not be the best for constructing an addiction-related brain network. It is believed that the more specific the brain network reflecting neuronal activity associated with addiction, the better the clustering analysis results. There could be other methods to select the ROIs more sensitive to relapse, which deserves our more attemptation in future. The second one that might influence the clustering results is the sensitivity of MRI index describing ROIs neuronal activity. The reason we chose DC as the MRI index was that it had stability and sensitivity when as a topological feature of graphics and few MRI indicators were suitable for clustering analysis reportedly. As for whether DC is the most suitable for clustering analysis, it is still unclear and needs more investigations. To explore the impact of other MRI parameters on cluster results is valuable. The third one is the relatively longer interval of a structured interview and a urine drug test for assessing illicit drug use, which may lead to a bias of relapse rate. Moreover, we confess that the final results could be impacted by the limitation of only having male subjects. It is necessary to use a more appropriate interview plan and gender-matched group to verify the reliability and repeatability of our findings in future.

In summary, we prospectively classified the former heroin addicts with stable MMT into a high or low relapse risk group by grouping their brain neural activity patterns with a neuroimaging data-driven method. The two groups shared the same demographic information, methadone use history and smoking use except total heroin consumption and relapse rate. In addition, it was found that HRG was more extensive in the neuropathological changes in the brain, and that NAc, vACC, hippocampus and amygdala had a role in the contribution to illegal drug use. These brain hubs may be selected as potential psychological and pharmacotherapeutic targets in future. Importantly this new strategy of identifying HRG with a data-driven approach measuring addiction-related network pattern rather than clinical phenomenology would suggest a great possibility of medication pathway and benefit the MMT patients in the optimizing treatment plan and improving therapeutic effect.