Participants
The present study included 173 university students with IGD (males/females: 98/75) and 232 HCs (males/females: 139/93). The numbers of males and females were matched between the two groups (χ2 =0.436, p =0.509). All the participants were free of any psychiatric/neurological disorders confirmed by structured psychiatric interviews (Mini International Neuropsychiatric Interview). No participant reported any illegal drug use and gambling, and excessive nicotine and alcohol uses. The diagnosis of IGD was based on the Young's Internet Addiction Test (IAT) (36) and DSM-5 criteria for IGD (37). Consistent with previous studies of IGD (4, 38-40), the diagnostic criteria for the IGD group included: 1) an IAT score >50; 2) a DSM-5 criteria score >5; and 3) spending most of their internet time in playing games. Participants with DSM-5 criteria scores <5 and IAT scores <50 were classified as the HC group. As shown in Table 1, the IGD group reported significant greater IAT scores and DSM-5 scores than the HC group.
Table 1. Demographic information and group differences
|
Items
|
IGD
(M=98, F=75)
|
HC
(M=139, F=93)
|
t
|
p
|
|
Age (years)
|
21.13±2.33
|
21.53±2.43
|
-1.663
|
0.097
|
|
IAT score
|
66.15±9.13
|
40.32±10.92
|
25.223
|
<0.001
|
|
DSM-5 score
|
6.08±1.13
|
2.51±1.37
|
26.447
|
<0.001
|
Table values: mean ± standard deviation
Abbreviations: IGD, internet gaming disorder; HC, healthy control; M, male; F, female; IAT, internet addiction test; DSM-5, the fifth edition of diagnostic and statistical manual of mental disorders.
RS-fMRI data acquisition
All participants underwent a RS-fMRI scan for 7 min using a 3 T scanner (Siemens Trio, Malvern, PA, USA) equipped for echo-planar imaging (EPI). During the scan, they were asked to open eyes and stay still and relaxed. The scan parameters were as follows: repetition time (TR) = 2000 ms, echo time (TE) = 30 ms, slice number = 33, interleaved sequence, slice thickness = 3 mm, voxel size = 3 × 3 × 3 mm3, field of view (FOV) = 220 × 220 mm2, flip angle = 90◦ and matrix of 64 × 64.
Preprocessing of RS-fMRI data
Data preprocessing was performed using DPABI V8.2 (http://rfmri.org/dpabi). The preprocessing steps included: 1) discarding the first 10 time points; 2) slice-timing; 3) realign for head-motion correction; 4) spatial normalization to Montreal Neurological Institute (MNI) standard space; 5) spatial smoothing (FWHM = 6 mm); 6) linear trend removal; 7) nuisance covariates regression, including head-motion covariates using the Friston 24-parameter model as well as signals from white matter, cerebrospinal fluid and global signals; 8) band-pass filtering with a range of 0.01-0.1 Hz.
Constructing normative models of FC
An overview of the normative modeling approach is provided in Figure 1; this approach has been described previously (22). First, ROI (bilateral NAcc)-based FC values were generated for each participant using DPABI v8.2. Second, warped Bayesian linear regression (41) was used to generate the normative models of FC at each specific voxel of the whole brain on the HC group (N=232), using age and gender as independent variables to predict FC, i.e., dependent variable, which was performed using the PCNtoolkit (Predictive Clinical Neuroscience toolkit) (42). The normative model of each voxel contains both the predicted FC and related predictive confidence at this voxel. The predictive confidence could be interpreted as centiles of variation within the norm population. Third, the normative model helps us place each participant within the normative range and then quantify their deviations (i.e., Z scores) of FC from the healthy range at each specific brain voxel. One participant-specific (i) Z score at each voxel (j) is calculated using the equation below:

This equation combines three sources of information: 1) the difference between true response (
) and predicted response (
); 2) the predicted variance of the specific voxel (
); 3) the variance of the normative data of the specific voxel (
). Using this equation, a normative probability map (NPM) with all the Z scores of all the brain voxel for each participant is generated. Thus, the NPM of one participant provides a statistical estimate of how much the FC of the participant differs from the FC of healthy normative pattern at each voxel, i.e., Z scores. Lastly, to identify brain voxels with significant deviations from healthy normative pattern across the whole brain of each participant, the NPM of each participant was thresholded using FDR (false discovery rate) at p<0.01 (22). To examine the spatial spread of the significant deviations among the IGD group, an overlap map by counting all the participant-level FDR-corrected NPMs was generated. This can be used to identify which brain regions had positive deviation (increased FC) or negative deviation (decreased FC) among the IGD group compared to the HC group.
Estimating regional deviations for each participant
To better depict the important deviations of each participant’s FC, we calculated one summary index for each participant, i.e., regional deviation to capture the maximum deviation across specific brain regions or functional networks. To calculate regional deviations, we first parcellated the whole brain into 10 functional networks (the FPN, BGN, limbic network, default mode network, medial frontal network, motor network, visual association network, visual network I and II, and cerebellum network) from Shen’s brain functional atlas (32) (https://bioimagesuiteweb.github.io/ webapp/connviewer.html?species5human). Then, the trimmed mean of 1% of top positive/negative deviations across each functional network was taken as the regional deviation for each region of each participant. In the end, the greater trimmed mean deviation between positive and negative deviations was taken as the representative deviation for each participant in one functional network.
Cluster analysis
Cluster analysis was applied to identify subtypes of IGD participants. In the current study, the data-driven K-means cluster analysis (43) was conducted using the “kmeans” function in MATLAB to classify the IGD participants based on their regional deviations of the cognitive control and reward networks, i.e., FPN and BGN. The optimal number of clusters was determined by two methods: 1) the average Silhouette score for each value of cluster numbers, and the solution with the highest Silhouette score was considered as the optimal clustering configuration (44); 2) the within-cluster sum of squares (WCSS) and the elbow method, where the point of inflection indicated a suitable number of clusters (45).