Introduction
Methamphetamine use disorder (MUD) has emerged as a critical public health crisis and a leading contributor to drug-related mortality (Xing et al., Reference Xing, Mohiuddin, Bhuiyan, Hossain, Al-Yafeai, Mosa Faisal, Goeders, Conrad, Vanchiere, Patterson, Kevil and Bhuiyan2024). Chronic methamphetamine exposure induces severe mental and physical health consequences, marked by widespread cortical thinning in critical brain regions, including the prefrontal, cingulate, parietal, and temporal cortices (Koester et al., Reference Koester, Tittgemeyer, Wagner, Becker, Gouzoulis-Mayfrank and Daumann2012; Nie et al., Reference Nie, Zhao, Wen, Luo, Ju, Ren, Wu and Li2020; Petzold et al., Reference Petzold, Dean, Pochon, Ghahremani, De La Garza and London2022). Neuroimaging studies demonstrate that the severity of cortical volume loss escalates with both duration and cumulative drug exposure, establishing this structural decline as a progressive neuropathological process (Joo et al., Reference Joo, Lee, Hwang, Kim, Cheon, Lee, Kim, Yurgelun-Todd, Renshaw, Yoon and Lyoo2023; Lyoo et al., Reference Lyoo, Yoon, Kim, Lim, Choi, Kim, Hwang, Jeong, Cho, Chung and Renshaw2015; Nie et al., Reference Nie, Ghahremani, Mandelkern, Dean, Luo, Ren, Li and London2021; Ruan et al., Reference Ruan, Zhong, Yang, Fan, Zhuang, Du, Jiang and Zhao2018). However, the mechanisms driving the regionally selective vulnerability of cortical architecture in MUD remain unresolved. Identifying the neurobiological determinants of this patterned morphological degeneration could reveal critical therapeutic entry points for halting disease progression.
Emerging evidence suggests that the formation and maintenance of addiction are driven by disruptions in brain connectome networks (Joutsa et al., Reference Joutsa, Moussawi, Siddiqi, Abdolahi, Drew, Cohen, Ross, Deshpande, Wang, Bruss, Stein, Volkow, Grafman, van Wijngaarden, Boes and Fox2022; Zhang & Volkow, Reference Zhang and Volkow2019). Patterns of cortical thinning in MUD are highly organized and circumscribed by specific networks, supporting the notion that connectome architecture plays a significant role in the pathological process. For instance, brain volume loss in MUD appears to be highly organized and confined to specific functional networks such as default mode network (DMN) (Joo et al., Reference Joo, Lee, Hwang, Kim, Cheon, Lee, Kim, Yurgelun-Todd, Renshaw, Yoon and Lyoo2023). Another study reported that cortical thinning and the breakdown of white matter integrity occur concurrently in patients with MUD (Lyoo et al., Reference Lyoo, Yoon, Kim, Lim, Choi, Kim, Hwang, Jeong, Cho, Chung and Renshaw2015). The functional connections of region showing structure alterations and the medial prefrontal cortex could predict relapse status (Geng et al., Reference Geng, Hu, Gu, Salmeron, Adinoff, Stein and Yang2017). These findings support that brain network organization may be linked to structural alterations in MUD.
Network-dependent hypothesis supported that the morphological abnormalities observed in the brain may lie in the influence of the brain network or connectome (Jucker & Walker, Reference Jucker and Walker2018; Seguin, Sporns, & Zalesky, Reference Seguin, Sporns and Zalesky2023). Specifically, neural connections might help spread pathological processes. This occurs by transporting toxic markers to distant cortical areas. Alternatively, if the network architecture is impaired, it can disrupt normal communication between brain regions (Wannan et al., Reference Wannan, Cropley, Chakravarty, Bousman, Ganella, Bruggemann, Weickert, Weickert, Everall, McGorry, Velakoulis, Wood, Bartholomeusz, Pantelis and Zalesky2019). Under this hypothesis, various models have been developed to explore the impact of brain networks on regional morphological changes under diseased conditions (Zhou, Gennatas, Kramer, Miller, & Seeley, Reference Zhou, Gennatas, Kramer, Miller and Seeley2012). One such model is the nodal stress model, which hypothesizes that highly interconnected regions appear to be more vulnerable to structural damage due to activity-related ‘wear and tear’ metabolic demands (Brown et al., Reference Brown, Deng, Neuhaus, Sible, Sias, Lee, Kornak, Marx, Karydas, Spina, Grinberg, Coppola, Geschwind, Kramer, Gorno-Tempini, Miller, Rosen and Seeley2019). For example, studies show that progressive cortical atrophy is spatially correlated with node centrality in several diseases. These include neurodegenerative diseases such as Parkinson’s disease (Zeigha et al., Reference Zeighami, Ulla, Iturria-Medina, Dadar, Zhang, Larcher, Fonov, Evans, Collins and Dagher2015) and Alzheimer’s disease (Zhou et al., Reference Zhou, Gennatas, Kramer, Miller and Seeley2012). Psychiatric conditions such as major depressive disorder (Ha et al., Reference Han, Fang, Zheng, Li, Zhou, Sheng, Wen, Liu, Wei, Chen, Chen, Cui, Cheng and Zhang2023) and schizophrenia (Georgiadis et al., Reference Georgiadis, Lariviere, Glahn, Hong, Kochunov, Mowry, Loughland, Pantelis, Henskens, Green, Cairns, Michie, Rasser, Catts, Tooney, Scott, Schall, Carr, Quide and Kirschner2024) also show this correlation. In addition, studies believed that white matter connectome architecture can also contribute to the progression of atrophy, where connected regions exhibit similar atrophy patterns (here, termed nodal-neighbor common atrophy model). This model has revealed regional atrophy was linked to atrophy of connected neighbors under conditions such as schizophrenia (Shafiei et al., Reference Shafiei, Markello, Makowski, Talpalaru, Kirschner, Devenyi, Guma, Hagmann, Cashman, Lepage, Chakravarty, Dagher and Misic2020) and frontotemporal dementia (Shafiei et al., Reference Shafiei, Bazinet, Dadar, Manera, Collins, Dagher, Borroni, Sanchez-Valle, Moreno, Laforce, Graff, Synofzik, Galimberti, Rowe, Masellis, Tartaglia, Finger, Vandenberghe, de Mendonca and Initiative2023). Network-dependent hypotheses provide a valuable framework for exploring the development of cortical morphology in healthy individuals, the progression of diseases, and their heterogeneity (Jiang et al., Reference Jiang, Palaniyappan, Luo, Chang, Zhang, Tang, Zhang, Li, Zhou, Yu, Li, An, Zhou, Huang, Tsai, Lin, Cheng, Wang, Yao and Feng2024; Li et al., Reference Li, Long, Sheng, Du, Qiu, Chen and Liao2024; Liang et al., Reference Liang, Sun, Liao, Lei, Xia, Duan, Zeng, Li, Xu, Men, Wang, Tan, Gao, Qin, Tao, Dong, Zhao and He2024). However, limited research has focused on the role of these models in cortical morphology alterations associated with MUD.
In this study, we aimed to explore the associations between cortical thickness reduction in individuals with MUD and the healthy brain connectome, delineating a network-dependent spatial pattern of cortical thickness thinning. First, we spatially correlated nodal centrality with the cortical thinning patterns in patients with MUD. Second, we examined the relationship between regional cortical thinning and the collective cortical thickness reduction of neighboring regions, weighted by the healthy brain connectome. The relationship between network-weighted cortical abnormality and clinical phenotype scores was estimated by partial least squares (PLS) regression. To assess the effects of drug use duration and dosage on network-based cortical thickness alterations, we performed clustering analysis on MUD subgroups using individual network-weighted cortical abnormality scores. Building on prior studies of neurodegenerative and psychiatric disorders, we hypothesized that cortical thickness reductions in MUD are shaped by connectome architecture. Specifically, brain regions with high network centrality may be more susceptible to morphometric alterations, and such abnormalities may propagate to connected areas via structural pathways.
Materials and methods
Study samples
The study recruited 139 individuals with MUD from the Kangda Voluntary Drug Rehabilitation Centre and Hunan Brain Hospital in Hunan Province, China. We also recruited 119 healthy controls (HCs) from local communities via online advertisements. All participants were right-handed and aged between 18 and 45 years. The patients with MUD met the Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition, Text Revision (DSM-IV-TR) criteria for substance dependence (methamphetamine) as determined by the Structured Clinical Interview for DSM Disorders. Exclusion criteria included learning disabilities, neurological or psychiatric conditions (excluding substance-induced symptoms), history of traumatic brain injury (>10 minutes unconsciousness), familial predisposition to heritable mental disorders, recent neuromodulatory interventions within 3 months, pregnancy, MRI contraindications, and dependence on non-nicotine substances.
The entire study protocol was approved by the institutional review board of the Second Xiangya Hospital (No. S163, 2011) and Sir Run Run Shaw Hospital, Zhejiang University School of Medicine (No. 2023-826-01), in accordance with the Declaration of Helsinki (1975, revised 2008). Written informed consent was obtained from all participants prior to their involvement in the study.
Clinical measurements
The sociodemographic characteristics and drug use patterns of each participant were evaluated, including sex, age, years of education, duration of drug use (in months), and lifetime drug consumption. Among the patients with MUD, we used a 10-point visual analogue scale to assess methamphetamine craving, where 0 represents the lowest craving level and 10 indicates the highest intensity of craving (Mezinskis, Honos-Webb, Kropp, & Somoza, Reference Mezinskis, Honos-Webb, Kropp and Somoza2001). Additionally, the severity of psychosis, depression, and anxiety symptoms was measured using the positive and negative syndrome scale (PANSS), the 24-item Hamilton depression rating scale (HDRS), and the 14-item Hamilton anxiety rating scale (HARS) (Kay, Fiszbein, & Opler, Reference Kay, Fiszbein and Opler1987; Maier, Buller, Philipp, & Heuser, Reference Maier, Buller, Philipp and Heuser1988; Williams, Reference Williams1988), respectively. To ensure the reliability of these assessments, all psychiatrists involved in the study underwent extensive training in administering PANSS, HDRS, and HARS before the study began.
MRI data acquisition
The T1-weighted anatomical images were acquired using a Siemens Magnetom Trio 3 T scanner (Allegra; Siemens, Erlangen, Germany) with a 64-channel head coil. The scanning parameters were as follows: a 3D MPRAGE sequence with a 256 × 256 mm2 field of view, 1 mm slice thickness with no gap, repetition time (TR) of 2000 ms, repetition time (TE) of 3.7 ms, flip angle of 8, and a total of 176 slices.
Image processing and cortical thickness calculation
All neuroimage data were preprocessed via FreeSurfer software (version 7.1, https://surfer.nmr.mgh.harvard.edu). The preprocessing of 3D-T1 images employed a standard auto-reconstruction algorithm that included several key steps: motion correction, intensity normalization, non-brain tissue removal, automated transformation into Talairach space, segmentation of white and gray matter, volumetric structure delineation, and tessellation. In addition, automated topology correction and surface deformation techniques were applied. Cortical thickness (CT) was estimated at each surface vertex by calculating the shortest distance between the gray and white matter surfaces. For each participant, CT measurements were obtained for 68 cortex regions according to the Desikan-Killiany parcellation atlas (Desikan et al., Reference Desikan, Ségonne, Fischl, Quinn, Dickerson, Blacker, Buckner, Dale, Maguire, Hyman, Albert and Killiany2006). The quality of cortical parcellations was visually assessed, with histograms generated for all regions to enable a comprehensive visual examination. No subjects were excluded from the analysis.
We identified cortical thinning patterns associated with MUD by employing a general linear model to assess regional CT differences between patients and healthy controls. Covariates such as sex, age, years of education, and intracranial volume (ICV) were included in the analysis. Bonferroni correction was used for multiple comparisons across 68 brain regions, resulting in an adjusted significance threshold of p < 0.05/68. Statistical t-maps were extracted and converted to z statistics for following analysis.
Structural, morphological, and functional network reconstruction
We utilized the ENIGMA toolbox (https://enigma-toolbox.readthedocs.io/en/latest/) to derive resting-state functional connectivity and structural connectivity matrices. These matrices, sourced from healthy adults (n = 207, 40% males, mean age = 28.73) in the Human Connectome Project (HCP; http://www.humanconnectome.org/), have been employed in prior studies exploring the impact of brain networks on cortical alterations (Georgiadis et al., Reference Georgiadis, Lariviere, Glahn, Hong, Kochunov, Mowry, Loughland, Pantelis, Henskens, Green, Cairns, Michie, Rasser, Catts, Tooney, Scott, Schall, Carr, Quide and Kirschner2024; Hettwer et al., Reference Hettwer, Lariviere, Park, van den Heuvel, Schmaal, Andreassen, Ching, Hoogman, Buitelaar, van Rooij, Veltman, Stein, Franke, van Erp, Group, Group, Group, Group, Group and Valk2022). A group-average functional connectivity matrix was generated by calculating pairwise correlations between 68 cortical regions of Desikan-Killiany atlas. Additionally, the group-average structural connectivity matrix was computed by applying a distance-dependent threshold to the number of streamlines connecting these regions. For more details, see Supplementary Materials.
For the morphological network analysis, CT measurements were obtained from 68 cortical regions, as defined by the Desikan-Killiany parcellation scheme, using data from our healthy control group. A group-level morphological network was constructed by calculating Pearson correlations between the CT values of each pair of regions. Fisher r-to-z was applied to the correlation coefficients to improve normality.
Network analysis of cortical alterations
The modal stress model and nodal-neighbor common atrophy model were modeled the relationship between normal connectome across three modality and cortical thickness changes in MUD (Supplementary Figure S1). In the nodal stress model (Zhou et al., Reference Zhou, Gennatas, Kramer, Miller and Seeley2012), we predicted that higher degree centrality of a node would exhibit larger cortical thinning in MUD. Weighted degree centrality was calculated using structural, morphological, and functional connectivity data by summing the values of all cortico-cortical connections for each region, respectively. Spatial correlation method was used to examine relationship between MUD-related cortical thinning patterns (t-value map) and normative weighted degree centrality profiles using Pearson correlation coefficients. The significance of associations of nodal centrality and cortical thinning were tested by spatial autocorrelation permutation tests (spin tests) with 5000 times (Alexander-Bloch et al., Reference Alexander-Bloch, Shou, Liu, Satterthwaite, Glahn, Shinohara, Vandekar and Raznahan2018). Specifically, the centroids of the 68 cortical regions defined by the Desikan-Killiany atlas were projected onto a spherical surface. The surface was then randomly rotated across the sphere, preserving the spatial topology of the brain. After each rotation, regional values were reassigned according to the closest parcel on the rotated sphere. This process was repeated 5000 times to generate a null distribution that maintains spatial dependencies. The empirical p-value was computed as the proportion of permuted correlations that exceeded the observed correlation.
In the nodal-neighbor common atrophy model (Shafiei et al., Reference Shafiei, Markello, Makowski, Talpalaru, Kirschner, Devenyi, Guma, Hagmann, Cashman, Lepage, Chakravarty, Dagher and Misic2020), we predicted that nodes exhibiting more severe cortical changes are connected to neighboring regions with greater cortical thinning. At the group level, we used healthy connectivity networks from three modalities to estimate the morphological abnormalities of the neighbors for each region. The structurally weighted neighbor cortical thinning level of a region was then calculated as the average cortical thinning values of all the brain regions structurally connected to that node:
$$ {T}_i=\frac{i}{N_i}{\sum}_{j\ne i,j=1}^{N_i}{t}_j $$
where Ti represents the average neighbor cortical thinning value of structural neighbors of a given node i, tj is the cortical thinning of the jth neighbor structurally connected to node i, and Ni is the total number of regions structurally connected to node i. Normalization by the term Ni ensures that the estimated neighbor cortical thinning value is independent of the nodal degree. The calculation of neighbor cortical thinning excludes self-connections (i.e., i ≠ j).
The morphological-weighted connected neighbor cortical thinning level was assessed using Equation (1), with the exception that cortical thinning values were weighted by the strength of morphological connectivity between nodes i and j (SMNij):
$$ SMN-{T}_i=\frac{i}{N_i}{\sum}_{j\ne i,j=1}^{N_i}{t}_j\times {SMN}_{ij} $$
where SMN-Ti represents the average neighbor cortical thinning value of morphologically defined neighbors of a given node i, and SMNij is the mean morphological connectivity between nodes i and j.
Similarly, the functional-weighted connected neighbor cortical thinning level was assessed using Equation (2), with the exception that cortical thinning values were weighted by the strength of functional connectivity between nodes i and j (FCij):
$$ FC-{T}_i=\frac{i}{N_i}{\sum}_{j\ne i,j=1}^{N_i}{t}_j\times {FC}_{ij} $$
where FC-Ti represents the average neighbor cortical thinning value of morphologically defined neighbors of a given node i, and FCij is the mean functional connectivity between nodes i and j.
For structural-, morphological- and functional-weighted neighbor cortical thinning estimations, neighbors of a region were determined as a node structurally connected to that region. Pearson correlation coefficients were employed to assess the relationship between nodal cortical thinning and its structurally, morphological- and functionally weighted connected neighbors cortical changes level, respectively. The spatial autocorrelation-preserving permutation tests was determined the significance of associations of nodal centrality and cortical thinning (5000 repetitions).
Epicenter identification
To identify potential epicenters of cortical thinning in MUD, we adopted a nodal-neighbor cortical thinning ranking method (Shafiei et al., Reference Shafiei, Markello, Makowski, Talpalaru, Kirschner, Devenyi, Guma, Hagmann, Cashman, Lepage, Chakravarty, Dagher and Misic2020). Specifically, for each brain region, we first ranked the degree of cortical thinning (based on the group-level MUD versus control thickness difference) in ascending order. Separately, for each region, we computed the morphology-weighted average cortical thinning level of its structurally connected neighbors using the morphological covariance network and ranked these values in the same way. The final epicenter likelihood for each region was determined by averaging its own atrophy rank and its neighbor-based rank. Regions with the lowest mean ranks were considered the most probable epicenters, as they showed high levels of cortical thinning themselves and were surrounded by similarly affected neighbors. To assess the statistical significance of these rankings, we employed a spin permutation test with 5000 repetitions, which preserves the spatial autocorrelation structure of the cortical surface by rotating the cortical thinning map. Given that the morphological network best explained cortical thinning patterns in both the nodal stress and nodal-neighbor models (see Results), the epicenter analysis was conducted based on morphology-weighted neighbor thinning values.
Individual network-weighted cortical abnormality model
To examine the associations between network-weighted cortical alterations (participants × brain regions) and clinical factors (participants × measures) in patients, a multivariate statistical algorithm, PLS regression, was employed (McIntosh & Mišić, Reference McIntosh and Mišić2013). Personalized network-weighted cortical abnormality was defined using connectome-weighted cortical thinning. Specifically, we estimated individual cortical thinning maps for patients with MUD by measuring deviations from the normal distribution. A multiple regression model was developed to predict regional cortical thickness, incorporating sex, age, education level, and ICV as covariates. This model was then applied to estimate the cortical thickness for each patient with MUD. Regional W-scores were generated by calculating the difference between the observed and predicted cortical thickness, divided by the standard error of the model fit in controls, providing individualized cortical thickness maps. A positive W-score means cortical thinning in patients with MUD in relative to HC and vice versa. Personalized network-weighted cortical abnormality of a region was determined by the normalized W score of all regions that were connected to region i by a morphological connection.
$$ SMN-{W}_i=\frac{i}{N_i}{\sum}_{j\ne i,j=1}^{N_i}{w}_j\times {FC}_{ij} $$
where SMN-Wi represents the individual neighbor cortical thinning value of morphologically defined neighbors of a given node i, and W is individualized cortical thickness thinning.
PLS analysis was performed to investigate the association between individual network-weighted cortical abnormalities and clinical factors (including craving, PANSS, HARS, and HDRS). Network-weighted cortical abnormality maps served as predictor variables, while clinical scale scores were used as response variables. The significance of the variance explained by the PLS components was evaluated using permutation test with 5000 iterations. To model the effect of brain atrophy on clinical symptoms, Pearson correlation analyses were performed between network-weighted cortical abnormalities (brain scores) and the behavioral phenotypes (behavior scores) identified in the significant latent variable. The behavior loadings are denoted by the Pearson correlation between behavior phenotypic data and the brain scores.
Associations of network-based cortical abnormality and drug use
Previous studies have indicated that network-level gray matter abnormalities progress from the early stages of illness to more severe stages in neurodegenerative and psychiatric diseases (Chopra et al., Reference Chopra, Segal, Oldham, Holmes, Sabaroedin, Orchard, Francey, O’Donoghue, Cropley, Nelson, Graham, Baldwin, Tiego, Yuen, Allott, Alvarez-Jimenez, Harrigan, Fulcher, Aquino and Fornito2023; Wannan et al., Reference Wannan, Cropley, Chakravarty, Bousman, Ganella, Bruggemann, Weickert, Weickert, Everall, McGorry, Velakoulis, Wood, Bartholomeusz, Pantelis and Zalesky2019; Zhou et al., Reference Zhou, Gennatas, Kramer, Miller and Seeley2012). Building on these findings, we hypothesized that drug use, such as duration and dose of drug, might influence cortical thickness alterations at the network level, resulting in heterogeneous subgroups. To investigate the impact of drug use on network-based cortical thickness changes, we performed k-means clustering to identify subgroups related to drug use. K-means clustering divides the dataset into k non-overlapping groups by minimizing variance within clusters and maximizing differences between clusters, based on a predefined distance metric. This method has been widely used in psychiatric research for data-driven biotype discovery (Ha et al., Reference Han, Fang, Zheng, Li, Zhou, Sheng, Wen, Liu, Wei, Chen, Chen, Cui, Cheng and Zhang2023; Li et al., Reference Li, Long, Sheng, Du, Qiu, Chen and Liao2024). We first constructed a 139 × 139 similarity matrix by calculating the Pearson correlation coefficients between SMN-weighted cortical thinning profiles for all MUD participants. This matrix captures the spatial similarity in cortical thinning patterns between every pair of individuals, where higher values indicate greater similarity. Based on this similarity matrix, k-means clustering using Euclidean distance were applied to identify potential subgroups. To enhance robustness and reproducibility, the clustering procedure was repeated 100 times with random initializations. The optimal number of clusters (k) was determined by calculating silhouette scores for values of k ranging from 2 to 10. The differences in methamphetamine duration and dose between the subgroups were assessed using a two-sample test.
Results
Sample characteristics
Table 1 presents the sociodemographic characteristics of patients with MUD and HC. There was no significant group difference in sex and age between individuals with MUD and HC. However, patients with MUD showed lower education level than HC. Patients with MUD showed higher anxiety and depression level in relative to HC.
Table 1. The demographic and clinical characteristics of patients with methamphetamine use disorder and HC
Note: Data are showed as the mean ± SD. HC, healthy control; HARS, Hamilton anxiety rating scale; HDRS, Hamilton depression rating scale; MUD, methamphetamine use disorder; PANSS, positive and negative syndrome scale; VAS, visual analogue scale.
Cortical thickness reductions in patients with MUD
Compared to HC, a total of 28 cortex regions were found to exhibit significant cortical thickness thinning in patients with MUD (p < 0.05, Bonferroni corrected). We did not find any increased cortical thickness in methamphetamine dependent groups. As illustrated in Figure 1, pronounced cortical thinning was observed in the frontal and temporal lobes, particularly in the superior frontal gyrus and middle temporal gyrus. To investigate whether this cortical thinning was predominantly associated with specific brain networks, we employed two brain atlases: (1) the Yeo atlas (Yeo et al., Reference Yeo, Krienen, Sepulcre, Sabuncu, Lashkari, Hollinshead, Roffman, Smoller, Zöllei, Polimeni, Fischl, Liu and Buckner2011), defined by intrinsic functional connectivity, and (2) the Von Economo-Koskinas atlas (Scholtens, de Reus, de Lange, Schmidt, & van den Heuvel, Reference Scholtens, de Reus, de Lange, Schmidt and van den Heuvel2018), based on cytoarchitectonic classes. Mean cortical alteration values were calculated within each brain system according to both atlases. A permutation test with 5000 repetitions was used to assess the significance of cortical alterations at the network level. The results revealed that patients with MUD exhibited greater cortical thinning in the frontoparietal (FPN) and DMNs. For the Von Economo-Koskinas atlas, more significant cortical changes were observed within the frontal and parietal cytoarchitectonic classes.
Figure 1. Cortical thickness reductions in patients with methamphetamine use disorder (MUD). (A) Maps displaying cortical morphological changes in both MUD patients and the control group. (B) Mean cortical alteration values assessed within resting-state networks based on the Yeo atlas. (C) Mean cortical alteration values calculated using the Von Economo-Koskinas atlas. The brain networks according to the Yeo atlas include: VIS (visual network), SMN (somatomotor network), DAN (dorsal attention network), FPN (frontoparietal network), VAN (ventral attention network), and DMN (default mode network).
Brain connectome shape cortical thinning
We employed two models to investigate the relationship between brain networks and cortical morphological alterations in patients with MUD. As described in the Figure 2, in the stress nodal model, we found that greater nodal degree centrality was significantly associated with more pronounced cortical thickness changes. This significant relationship was observed across different brain network modalities, including structural (r = 0.330, pspin = 0.007), morphological (r = 0.583, pspin < 0.001), and functional (r = 0.361, pspin = 0.028) degree centrality. In the nodal-neighbor common atrophy model, cortical thinning in a specific region was significantly linked to the cortical thinning of its connected neighbors across structural (r = 0.368, pspin = 0.020), morphological (r = 0.618, pspin < 0.001), and functional (r = 0.480, pspin = 0.003).
Figure 2. Network-based cortical thickness alterations. If regional morphometric changes depend on network connectivity, nodes with higher morphometric alterations exhibit high levels of normative network centrality and are connected to neighbors with greater nodal atrophy. (A) Schematic of epicenter identification: A region was defined as an epicenter if it exhibited high cortical thickness changes, and its connected neighbors also showed significant changes. (B) The T-value map of cortical thickness alterations is associated with structural, morphological, and functional connectivity of degree centrality. (C) Nodal thinning is related to alterations in its neighbors, defined by structural connectivity, and weighted by both morphological and functional connectivity. (D) The spatial distribution of the putative epicenters in MUD is shown. The upper panel illustrates the mean rankings mapped onto the brain surface, while the lower panel depicts the statistical significance of these rankings based on spin tests.
Disease of epicenters
Using a network-based node-neighbor ranking method, we identified an epicenter likelihood pattern (Figure 2A). This method determines the epicenter probability of a region by calculating the mean ranking based on both the nodal cortical thinning value and the thinning values of its neighboring regions. We found that regions in the bilateral superior frontal gyrus, the right rostral middle frontal gyrus, the right pars orbitalis, and the left middle temporal gyrus exhibited significantly higher mean rankings (Figure 2D, pspin < 0.001). This finding was consistent across both structural- and functional-weighted neighbor thinning analyses.
Mapping individual network-weighted cortical thinning on clinical symptoms
The PLS analysis identified the first latent variable, which accounted for 11.66% of the covariance between network-based cortical abnormality maps and clinical scores (p = 0.043, permutation test). As shown in Figure 3, this latent variable highlighted a cortical thinning pattern predominantly in prefrontal areas, closely resembling the epicenter likelihood distribution. Individuals exhibiting this thinning pattern experienced more severe psychosis, anxiety, and depression symptoms. This pattern was not related to craving.
Figure 3. The associations of individual network-weighted cortical alterations (W-scores) and clinical symptoms. (A) The first latent variable derived from a partial least squares (PLS) regression analysis linking individual W-scores. (B) The relationship between brain scores and behavior scores, where brain scores are defined by expression of individual W-scores on the first latent variable, and behavior scores represent the behavioral phenotype on the first latent variable. (C) The bar plot depicts the behavioral loadings, which were calculated based on the relationship between behavioral phenotypic data and brain scores.
Network-based cortical thinning identified two MUD biotypes
We employed an unsupervised approach to investigate the impact of drug use on cortical thinning at the network level. As shown in Figure 4, K-means clustering analysis identified two distinct MUD biotypes. Subgroup 1 (n = 94) displayed patterns of cortical thinning, whereas subgroup 2 (n = 45) exhibited a tendency toward cortical thickening. Subgroup 1 had significantly longer durations of methamphetamine use (t = 2.064, p = 0.040) and higher doses of methamphetamine use (t = 2.896, p = 0.004) compared to subgroup 2.
Figure 4. Clustering analysis. (A) The distribution of w-scores in subgroup 1 and subgroup. (B) The group difference for each subgroup in the duration of methamphetamine use. (C) The group difference between two subgroups in the total dose of methamphetamine use. Error bars represent the standard deviation of the data.
Discussion
The current study employed both the nodal stress model and the nodal-neighbor common atrophy model to investigate cortical thinning patterns in patients with MUD at the network level. The findings revealed that cortical thinning was constrained by healthy brain networks across distinct modalities. We observed widespread cortical thickness reduction in patients with MUD, accumulating in specific systems. Moreover, nodes exhibiting greater cortical thinning were associated with higher nodal degree centrality and connected neighbors’ cortical alterations. The frontal and temporal regions emerged as key epicenters. The individual network-based patterns of cortical alteration correlated with clinical symptom variability, distinguishing MUD subtypes associated with drug use.
Our findings align with previous studies (Joo et al., Reference Joo, Lee, Hwang, Kim, Cheon, Lee, Kim, Yurgelun-Todd, Renshaw, Yoon and Lyoo2023; Lyoo et al., Reference Lyoo, Yoon, Kim, Lim, Choi, Kim, Hwang, Jeong, Cho, Chung and Renshaw2015), showing widespread cortical thickness reduction in MUD patients, especially in the prefrontal and temporal areas. Cortical thinning was organized by specific functional brain networks, including the FPN and DMN. It was also structured by cytoarchitectonically defined frontal and parietal regions. These loci anchored within these functional networks, mainly in the prefrontal areas, were linked to inhibitory control and motivation (Friedman & Robbins, Reference Friedman and Robbins2022; Munakata et al., Reference Munakata, Herd, Chatham, Depue, Banich and O’Reilly2011). Dysfunctions in these behaviors contribute to the initiation and maintenance of addiction (Ceceli, Bradberry, & Goldstein, Reference Ceceli, Bradberry and Goldstein2022; Goldstein & Volkow, Reference Goldstein and Volkow2011). Combining these insights with epicenter analysis, we identified the dorsolateral prefrontal cortex as the likely spatial origin of cortical thinning in MUD, serving as a central hub for the network-driven spread of cortical changes. This suggested that damage to the prefrontal cortex is a critical factor in both vulnerability to and resilience against drug addiction (Ersche et al., Reference Ersche, Meng, Ziauddeen, Stochl, Williams, Bullmore and Robbins2020).
Our results further elucidated the relationship between brain networks and cortical thickness alterations, providing evidence that cortical morphological changes in MUD may be constrained by network architecture through both the nodal stress and nodal-neighbor models. According to the nodal stress hypothesis, highly interconnected nodes were positively correlated with greater degrees of deformation, suggesting that atrophy in patients with MUD is shaped by the topology of the brain connectome. In line with our findings, previous studies have reported alterations in small-world properties and modularity across various substance use disorders, including nicotine (Fan et al., Reference Fan, Zha, Liu, Wei, Wang, Song, Lv, Ren, Hong, Gou, Zhang, Chen, Zhou, Pan and Zhang2023), alcohol (Ottino-González & Garavan, Reference Ottino-González and Garavan2022), methamphetamine (Liu et al., Reference Liu, Li, Zhang, Wang, Wang, Chen, Zhu, Shi, Wang and Li2022), and heroin (Yuan et al., Reference Yuan, Qin, Liu, Guo, Dong, Sun, Zhang, Liu, Wang, Wang, Li, Yang, von Deneen, Gold, Liu and Tian2010). Additionally, cortical thickness alterations in MUD were influenced by white matter architecture, as significant relationships were found between nodal cortical changes and those of their connected neighbors. This may be due to impaired transneuronal transport of trophic factors, leading to brain tissue loss (Zalesky et al., Reference Zalesky, Pantelis, Cropley, Fornito, Cocchi, McAdams, Clasen, Greenstein, Rapoport and Gogtay2015). Brain-behavior analysis revealed that individuals with greater deviations in cortical thickness in the prefrontal areas from the normal model were associated with more severe clinical symptoms, such as psychosis and negative phenotypes. This suggests that network-based morphological alterations not only help trace disease epicenters but also map cortical alterations underlying clinical symptoms (Li et al., Reference Li, Wang, Mao, Huang, Qiu and Wang2023). Altogether, these findings indirectly support the idea that cortical morphological alterations in MUD patients align with the transneuronal network spread model, providing an explanation for the widespread, irregular, and discontinuous pattern of cortical thickness reduction in MUD. It is important to note that our findings remain correlational, making it difficult to disentangle the causal relationship between cortical thickness and network architecture in MUD. Beyond network-level factors, accumulating evidence suggests that neurobiological processes such as neuroinflammation and excitotoxicity may also play a role (Li et al., Reference Li, Ramos-Rolón, Kass, Pereira-Rufino, Shifman, Shi, Volkow and Wiers2024). For example, an early PET study using a Translocator Protein (TSPO) tracer reported elevated microglial activation (an indicator of neuroinflammation) in the midbrain and insular cortex of individuals with MUD (Sekine et al., Reference Sekine, Ouchi, Sugihara, Takei, Yoshikawa, Nakamura, Iwata, Tsuchiya, Suda, Suzuki, Kawai, Takebayashi, Yamamoto, Matsuzaki, Ueki, Mori, Gold and Cadet2008). Decreased glutamine levels have been observed in the medial prefrontal cortex of MUD individuals relative to HC, which may reflect altered glutamatergic neurotransmission (Wu et al., Reference Wu, Qi, Long, Liao, Wang, Xie, Liu, Hao, Tang, Yang, Liu and Tang2018). These findings underscore the importance of future multimodal studies to elucidate how inflammatory and excitotoxic processes interact with structural networks in driving cortical degeneration in MUD.
We further identified two distinct biotypes of MUD based on individual cortical thickness thinning scores at the network level. Consistent with our hypothesis, these two biotypes exhibited significant differences in drug use variables, such as the duration and dose of methamphetamine use. Individuals with longer histories and higher doses of methamphetamine use showed more extensive cortical alterations, suggesting a progressive network spread of gray matter changes in MUD. The dose-effect analysis further supports this notion, as it revealed that global mean cortical thickness in chronic methamphetamine users was significantly negatively correlated with the lifetime dose of drug use (Joo et al., Reference Joo, Lee, Hwang, Kim, Cheon, Lee, Kim, Yurgelun-Todd, Renshaw, Yoon and Lyoo2023; Lyoo et al., Reference Lyoo, Yoon, Kim, Lim, Choi, Kim, Hwang, Jeong, Cho, Chung and Renshaw2015). These findings are reminiscent of observations in neurodegenerative disorders and schizophrenia (Chopra et al., Reference Chopra, Segal, Oldham, Holmes, Sabaroedin, Orchard, Francey, O’Donoghue, Cropley, Nelson, Graham, Baldwin, Tiego, Yuen, Allott, Alvarez-Jimenez, Harrigan, Fulcher, Aquino and Fornito2023; Leuzy et al., Reference Leuzy, Binette, Vogel, Klein, Borroni, Tonietto, Strandberg, Mattsson-Carlgren, Palmqvist, Pontecorvo, Iaccarino, Stomrud, Ossenkoppele, Smith, Hansson and Alzheimer’s Disease Neuroimaging2023; Wannan et al., Reference Wannan, Cropley, Chakravarty, Bousman, Ganella, Bruggemann, Weickert, Weickert, Everall, McGorry, Velakoulis, Wood, Bartholomeusz, Pantelis and Zalesky2019), where patients in the later stages of illness tend to experience greater brain structural alterations. This suggests that the severity and duration of drug exposure play a critical role in the extent of cortical thinning observed in MUD.
The current study has several limitations. First, our imaging findings were derived from cross-sectional data, which limit our ability to make causal inferences about the progression and potential propagation of cortical thinning in MUD. As such, we cannot determine whether the observed morphometric alterations reflect consequences of prolonged drug exposure, predisposing factors, or dynamic disease processes. To validate the network spread hypothesis and clarify the temporal sequence of structural changes in relation to brain connectome architecture, future longitudinal studies are essential. Second, we derived group-level normative structural and functional connectivity matrices from HCP samples. Future research should investigate individual connectivity-based models of atrophy to evaluate how disease-specific connectomes influence brain atrophy. Third, all findings are based on a single type of addictive drug and a single dataset. Although we used a relatively large sample, the results should ideally be validated independently using a distinct MUD dataset. Additionally, the generalizability of network spreading of gray matter alterations in addiction should be further examined using different types of drugs. Fourth, individuals with MUD exhibited lower education levels than HC, which could potentially influence cortical thickness measurements. Previous studies have shown that variations in socioeconomic status are a key risk factor for substance use disorders and are associated with alterations in brain structure and function (Evans-Lacko et al., Reference Evans-Lacko, Aguilar-Gaxiola, Al-Hamzawi, Alonso, Benjet, Bruffaerts, Chiu, Florescu, de Girolamo, Gureje, Haro, He, Hu, Karam, Kawakami, Lee, Lund, Kovess-Masfety, Levinson and Thornicroft2018; Tian et al., Reference Tian, Young, Zhu, Xu, He, Chen, Hao, Jiang, Qiu, Chen and Qin2021). Future research is encouraged to further investigate the specific impact of educational attainment and other sociodemographic variables on brain alterations in MUD.
In conclusion, our results suggest that cortical connectivity networks can explain the irregular distribution of cortical thickness reductions in MUD. We identified the prefrontal areas as the likeliest ‘disease epicenters’. Network-based cortical thickness reductions were also correlated with the severity of individual clinical symptoms and could distinguish distinct diagnostic MUD subtypes associated with drug use. These findings provide network mechanistic insights into cortical morphological changes in MUD and addiction.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/S0033291725102067.
Acknowledgments
We sincerely thank all the participants who volunteered for this study.
Funding statement
This study was funded by the National Natural Science Foundation of China (grant no. U22A20302) and the STI 2030 – Major Projects (grant no. 2022ZD0211200). The funding sources had no involvement in the study’s design, data collection and analysis, decision to publish the findings, or manuscript preparation.
Competing interests
All authors unanimously declare the absence of any biomedical financial interests or potential conflicts of interest.
Open access
