In total,45 subjects divided into two groups (23 PSD, 22 PSND) participated in this study. This study was approved by the Ethics Committee of the Zhujiang Hospital of Southern Medical University, China. The subjects were all patients admitted to the Department of Rehabilitation Medicine of the Zhujiang Hospital of Southern Medical University from August 2017 to February 2019. The inclusion criteria included: All subjects met the criteria of the World Health Organization diagnosis of cerebral infarction; single infarction (3–5 cm). The exclusion criteria included: 1) patients with history of neurological disease (e.g. epilepsy); 2) patients with malingering, or any psychiatric disorders (e.g. mood and anxiety disorders, schizophrenia and psychosis). 3) antidepressant use at stroke onset or a family history of mental disorders; The following information was collected from each subject: demographics (i.e., age, sex, level of education, and living alone) and stroke severity as measured by NIHSS, MMSE, BI, and HAMD. Each subject received a detailed description and full instructions of the experimental procedure and signed an informed consent form. Clinical interviews were conducted by experienced neuropsychologists according to the DSM-V criteria. In the PSD group, participants had to meet DSM-V criteria for depressive disorder, with HAMD score of 17 or higher. Other participants were assigned to the PSND group. HAMD-24 was used to assess the severity of depression.
Experimental data were obtained with a Philips 3.0 Tesla magnetic resonance scanner (Royal Philips Electronics, Eindhoven, the Netherlands). During data acquisition, subjects were instructed to relax with their eyes closed, and to keep their heads still. Foam padding and earplugs were used to reduce head motion and scanner noise. The helmet of the magnetic resonance apparatus was used to limit movement of the patient's head. All subjects underwent T1-weighted magnetic resonance imaging and rs-fMRI scanning. High-resolution T1-weighted anatomical images were also acquired using a magnetization-prepared rapid gradient-echo sequence. TR/TE = 8.1/3.8 ms; flip angle = 90; field of view = 256×256 mm2; in-plane matrix = 256×256; voxel size =1×1×1 mm3, no slice gap, 176 sagittal slices) for each subject. The structural sequence took 6 mins. The DTI scan uses a single-shot spin echo-plane echo sequence (DW-EPI). (TR/TE = 12500/112 ms; flip angle = 90°; field of view =230×230 mm2; in-plane matrix = 128×128; voxel size = 1.88×1.88×2 mm3, slice thickness 3.0 mm). Blood oxygenation level–dependent functional imaging was acquired using a gradient echo-planar imaging sequence as follows: (TR/TE = =2000/35ms; flip angle = 90; field of view = 230×230 mm2; in-plane matrix = 64×64;, no slice gap; 33axial slices, slice thickness=2mm).The fMRI sequence took 8mins.
We studied the network topological changes in PSD and PSND. Subject-level structural connectome was derived from diffusion MRI data using probabilistic fiber tracking based on 90 regions of interests (see Methods for details). Subject-level functional connectome was derived from task-free fMRI data based on pairwise Pearson's correlations between the 90 regions of interests (24). Graph theoretical global-wise and nodal-wise metrics were computed for both connectomes. Furthermore, SC-FC coupling was evaluated as the Pearson's correlation between the structural matrices and functional matrices both at the individual and group levels. Alterations in structural and functional network topology metrics as well as SC-FC coupling measures were compared across groups and subsequently associated with HAMD factor score. (Fig.1)
TBSS Image processing
Diffusion-weighted images were preprocessed with functional MR imaging of the brain (FMRIB) Software Library (FSL) (http://www.fmrib.ox.au.uk/fsl). Following the TBSS guideline, we normalized individual FA volumes of the two groups to the MNI space via affine registration. A mean FA image was created by averaging all the registered FA images, and from the mean FA image, the FA skeleton which represents the center of all WM tracts was created. Then, the FA threshold (>0.2) was set on the skeleton to exclude the peripheral brain areas, including WM voxels only. All the registered FA images were further projected onto the skeleton, resulting in a 4D file of all skeletonized images from the individual subject. Then, we generated the mean diffusivity (MD) data by applying the FA nonlinear registration to the individual parametric maps and projecting them onto the skeleton.
Construction of the structural network
To construct the structural network, diffusion tensors were reconstructed based on DTI theory (25). Following DTI reconstruction, a streamline-based fiber tracking algorithm (26) was applied on voxel-wise diffusion tensors with the following tracking parameters: random whole-brain seeding, 200,000 reconstructed streamlines, anisotropy threshold of 0.15, angular threshold of 45°, and streamline length between 30 and 300 mm. All DTI reconstructions and streamline-based fiber tracking procedures were performed in PANDA (https://www.nitrc.org/projects/panda/). To construct the structural brain network, the 90 cerebral regions from the AAL template were used as nodes of a network (24). We quantified the network edge between two distinct AAL regions by calculating multiple DTI metrics along the interconnected streamlines. Because not all DTI metrics are positively associated with the strength or integrity of WM connections, to distinguish between natures of network edges, we referred to DTI metrics as SC metrics only if they are designed to reflect the strengths of WM connections.
Construction of the functional network
To construct the functional network, rs-fMRI data of individual patients were preprocessed using the DPARSF toolbox (27). Preprocessing comprised the removal of the first 10 volumes, slice timing correction, coregistration to 3D T1WI, nuisance signal regression (head motion, WM signals and cerebrospinal fluid signals), nonlinear spatial normalization using T1WI, and band-pass filtering (0.01–0.1 Hz). Images with head motion exceeding half of the voxel width (1.5 mm) were excluded from the study. On completing preprocessing, the mean time series of each AAL region was obtained by averaging the voxel-wise BOLD signals within the selected cerebral region (24).
Multiparametric graph theoretical analysis
Following the construction of each individual network, graph theoretical analysis was employed to characterize its structural and functional network topologies. A threshold connectivity metric is usually applied to filter out unwanted spurious edges of a network prior to graph theoretical analysis. We employed a multiple-sparsity thresholding method, rather than a single threshold, to reduce variations caused by different thresholding values (28). Different ranges of thresholding values were used for constructing the structural and functional networks.
For each individual sparsity threshold, the sparsity thresholding process was applied to the streamline count matrix such that connections with few streamlines were removed from the graph to match the designated graph sparsity. This process was applied for every sparsity threshold value, yielding a set of thresholded SC matrices corresponding to each sparsity value. These thresholded SC matrices were then used as masks to filter out the unwanted entries in other structural network matrices. Finally, structural network metrics were calculated from each thresholded matrix, and averaged metric values across all sparsity thresholds were used for subsequent analyses. Similarly, this method was used for functional network analysis;
Small-world properties were measured by the Cp and the Lp, which computed the respective extent of interconnectivity of a network at local and global levels (29). The small-worldness was used to measure the balance between integration and segregation of network (30). Besides, efficiency metrics were also used to measure the small-world properties. The Eglob quantifies the extent of information transmission at the global networks, and the Eloc quantifies the extent at individual node levels (31).
We investigated the topological properties of brain networks at both global and nodal levels using the GRETNA toolbox (32). The degree for each of the nodes was quantified and used here as a reliable index of regional connectivity. The small-world coefficient σ was used to evaluate small-world behavior. For the whole brain network, we also calculated local and global efficiency.
Individual correlation matrices were transformed to the binary format at a wide range of network sparsity levels for extensive evaluation. Network sparsity measures the percentage of the number of existing edges in all possible connections and is used as threshold. Small-world properties were compared between groups at each sparsity level.
A module in the complex network is defined as a subset of nodes tightly connected within the modules but sparsely connected between the modules. The modularity of a network was computed using a modified greedy optimization algorithm (33, 34), details of which are provided in Supplementary materials. The modules are non-overlapping, with each node assigned to only one module. The process of modularity optimization does not need to specify either the number or the size of the modules.
To measure the potential regional role played by each brain region, the DC for each region (node) were defined as the indices of their inter-module connection density (35). For example, for a node in modules, the DC will be close to 0 if all weights are largely intramodular (36).
We first computed global topological properties (SW, γ, λ, σ) of weighted functional and structural connectivity networks for each subject. We also calculated the AUC for global topological properties, which provides a summarized scalar for brain topological properties independent of single threshold selection. Then, a nonparametric permutation test method was performed to detect the significant group differences in global topological properties. A significance threshold of p<0.05 (FDR-corrected) was used for testing global topological properties, except AUC of each global topological property for which an uncorrected threshold of p<0.01 was used since no multiple comparisons were performed. before the permutation tests, age and gender for each network properties were regressed out as covariates by multiple linear regression analyses. Furthermore, the coupling of functional-structural connectivity was compared by using permutation tests.