## 2 Study Design
### Clustering WM Networks  
Firstly, the WM mask at the group level was made according to the segmentation results of T1-weighted images. The steps were as follows: For each subject, we binarized each voxel as WM, GM or CSF based on the probability value of the T1 segmentation images. Then, these masks were averaged among subjects to obtain the percentage of each voxel that was classified as WM or GM. For WM, a voxel with a mean value greater than 60% was defined as a WM voxel, indicating that at least 60% of subjects recognized that voxel as WM. Then, to further remove GM clumps deep in the brain, subcortical areas of the Harvard-Oxford template were removed from the WM mask [25]. Finally, in order to ensure consistency with the functional image, the WM mask was registered in the functional space and resampled to a voxel size of 3 \(\times\) 3 \(\times\) 3 mm\({}^{3}\). The resulting WM mask contained 17,564 voxels. In addition, gray matter mask generation is similar to WM, and we used a relatively loose threshold (20%) to define GM voxel, ending up with a GM mask that contains almost all GM voxels. Finally, we made sure there is no overlap between GM and WM masks.  
Secondly, the correlation matrix at the group level was calculated for clustering WM functional networks. The 17,564 voxels of WM mask were resampled to 4353 nodes using the inter-changing grid strategy [26], and the Pearson's linear correlation coefficient between the time series of each WM voxel. The time series of the sampled node voxels was calculated to obtain the correlation matrix (17,564 \(\times\) 4353). To obtain an unbiased group level correlation matrix, two averages were performed. Firstly, Matrix-I was obtained by averaging the correlation matrices of all subjects in PS group. Then, the correlation matrices of all subjects in PIS group were averaged to obtain Matrix-II. Finally, we took the average of Matrix-I and Matrix-II to get Matrix-III. Because the number of people in PS and PIS groups is not equal, if we directly average all subjects, the results may be weighted more heavily for groups with more subjects, resulting in bias.  
Finally, the cluster analysis is carried out on WM. In order to obtain different WM functional networks, we use an iterative clustering algorithm called K-means clusteringalgorithm. Firstly, K-means clustering is performed on matrix Matrix-III, and the number of clusters K is traversed from 2 to 22. Previously, excellent researches indicated that clustering according to different features should provide approximately similar results for numbers of clusters [27; 28; 29]. In order to ensure the stability of clustering results, Matrix-III (17,564 \(\times\) 4353) is randomly divided into four submatrices, namely Matrix-III.1, Matrix-III.2, Matrix-III.3 and Matrix-III.4, with the size of each submatrix being 17,564 \(\times\) 1088. The same cluster analysis was performed for each submatrix separately, and the similarity between cluster results was measured with Dice coefficient. The averaged Dice's coefficient was used to assess the stability of the number of clusters. Finally, the maximum K value whose Dice coefficient is greater than 0.9 is selected as the final number of clusters.  
In view of the reliable and highly consistent results of clustering GM fMRI signals in many studies, this study did not conduct separate clustering analysis of GM voxels. Instead, a GM map in previous literature was selected, which was obtained by using the K-means clustering algorithm mentioned in this study. The results were consistent with the existing brain function map.