An imbalance between functional segregation and integration in patients with pontine stroke: A dynamic functional network connectivity study.

BACKGROUND: Previous studies on brain functional connectivity have revealed the neural physiopathology in patients with pontine stroke (PS). However, those studies focused only on the static features of intrinsic fluctuations, rather than on the time-varying effects throughout the entire scan. In the present study, we sought to explore the underlying mechanism of PS using the dynamic functional network connectivity (dFNC) method.
METHODS: Resting-state functional magnetic resonance imaging (fMRI) data were collected from 58 patients with PS and 52 healthy controls (HC). Independent component analysis (ICA), the sliding window method, and k-means clustering analysis were performed to extract different functional networks, to calculate dFNC matrices, and to estimate distinct dynamic connectivity states. Additionally, temporal features were compared between the two groups in each state to explore the brain's preference for different dynamic connectivity states in PS, and global and local efficiency were compared among states to explore the differences of topologic organization across different dFNC states. The correlations between clinical scales and the temporal features that differed between the two groups also were calculated.
RESULTS: The dFNC analyses suggested four recurring states; in two of these states, the PS group showed a different duration from that of the HC group. Patients with PS spent significantly more time in a sparsely connected state (State 1), which was characterized by relatively low levels of connectivity within and between all brain networks. In contrast, patients with PS spent significantly less time in a highly segregated state (State 2), which was characterized by high levels of positive connectivities within primary perceptional domains and within higher cognitive control domains, and by high levels of negative inter-functional connectivities (inter-FCs) among primary perceptional and higher cognitive control domains. Additionally, the dwell time in State 2 was positively correlated with HC group's long-term memory scores in the Rey Auditory Verbal Learning Test (RAVLT-L), whereas there was no correlation between the State-2 dwell time and RAVLT-L scores in the PS group. Furthermore, the sparsely connected state and the highly segregated state mentioned above had the highest global efficiency and the highest local efficiency among the four states, respectively.
CONCLUSIONS: In summary, we observed a preference in the aberrant brain for dynamic connectivity states with different network topologic organization in patients with PS, indicating abnormal functional segregation and integration of the whole brain and confirming the imperfection of functional network connectivity in patients with PS. These findings provide new evidence for the dynamic neural mechanisms underlying clinical symptoms in patients with PS.

Introduction

In brainstem infarcts, the pons is more susceptible than other structures, with an incidence of 73.3–81.8%; damage to the pontine structures can lead to hemiparesis, dysarthria and dysphagia, ataxia, and other clinical symptoms, seriously affecting the patient's quality of life(Pontine Infarcts and Hemorrhages, 2012, Maeshima et al., 2012, van Mierlo et al., 2015). However, the mechanism of functional deficits in patients with pontine stroke (PS) remains unclear. Further exploration of the neural physiopathology process in patients with PS is needed to improve our understanding of clinical manifestations, promote diagnostic accuracy, and provide appropriate treatment options in patients with PS.

Functional connectivity (FC) can be defined as ‘the temporal correlation of a neurophysiological index measured in different brain areas’ (Friston et al., 1993). In resting-state functional magnetic resonance imaging (fMRI), the blood oxygen level-dependent (BOLD) signal is considered to be a neurophysiological index, exhibiting low-frequency spontaneous fluctuations in the resting brain(Greicius, 2008). Some studies have explored FC alterations in patients with PS compared to healthy control (HC) patients. For example, a seed-based analysis (SBA) of early pontine infarction showed increased voxel-mirrored homotopic connectivity (VMHC) in the hippocampus/amygdala and frontal pole, and the VMHC was negatively correlated with the later motor assessment(Shan et al., 2018). An independent component analysis (ICA) of PS reported decreased FC in the left medial prefrontal gyrus in the anterior default mode network, right precuneus and right posterior cingulate cortex in the posterior default mode network, and left middle cingulate gyrus in the sensorimotor network(Chen et al., 2019). In addition to SBA and ICA, another promising technique for studying functional connectivity is graph theoretic analysis, which provides a new perspective for evaluating the property of integration and separation of brain networks(Smitha et al., 2017). However, few researchers (to our knowledge) have used this method to explore brain topological properties in patients with PS, thus ignoring the information transfer efficiency, which emphasizes the information flow across the entire network or across the portion of a network related to a particular function(Kawagoe et al., 2017). Additionally, the studies mentioned above assumed that FC remained constant, without considering the time-varying effects over the entire resting-state scanning.

Recent studies have indicated that resting-state FC (rsFC) is not static but constantly dynamic during scan acquisition(Gonzalez-Castillo et al., 2015, Hutchison et al., 2013, Calhoun et al., 2014, Chang and Glover, 2010). The most common approach to exploring dynamic FC is the sliding window technique, in which a time window is selected and then shifted by a fixed number of data points in the whole-scan series(Filippi et al., 2019). ICA is a data-driven method based on a blind source-separation algorithm(McKeown, 1998, Bell and Sejnowski, 1995, Boly et al., 2008); which does not require prior knowledge for the identification of regions of interest and can identify spatially segregated and functionally homogeneous intrinsic connectivity networks (ICNs)(Chen et al., 2018). The combination of the sliding window approach and ICA can identify spatial maps of each windowed BOLD signal(Jones et al., 2012, Premi et al., 2019), permitting measurement of the dynamic connectivity between ICNs, i.e., dynamic functional network connectivity (dFNC). The dFNC approach has been shown to capture valuable information that would be missed in static FNC, providing higher overall accuracy than the static FNC approach(Rashid et al., 2016, Nomi et al., 2016, Cetin et al., 2016). Recently, researchers have found that dFNC may be associated with some neurophysiological processes in the human brain, including the aging process(Chen et al., 2018, Xia et al., 2019), imagination and creativity(Beaty et al., 2018); spontaneous cognition and behavioral transitions(Allen et al., 2014). Furthermore, a growing number of studies has suggested that the properties of dFNC may serve as distinctive characteristic markers for neuropsychiatric diseases. For example, Bhinge et al. revealed that patients with schizophrenia tended to remain in a state corresponding to a hyperconnected brain network(Bhinge et al., 2019); Zhu et al. found decreased cerebellar-cerebral dynamic FC of the cerebellar subregions connecting with the executive, default-mode and affective-limbic networks in patients with major depressive disorder(Zhu et al., 2020); and Klugah-Brown et al. reported that patients with epilepsy demonstrated decreased dFNC between the frontoparietal system and other systems in all states(Klugah-Brown et al., 2019). Recent years also have witnessed the use of graph theoretic analysis as a powerful tool to characterize network states of dFNCs(Zoltowski et al., 2014). For example, migraine patients were found to have significantly lower information transfer efficiency in both global and local dFNCs(Tu et al., 2019). Additionally, patients with major depression disorder also represented reduced global and local efficiencies compared to healthy controls when the topological properties of dFNC states were analyzed(Zhi et al., 2018). However, to date, few studies (to our knowledge) have investigated the abnormalities of time-varying brain activity and topological properties in patients with PS.

In this study, we sought to investigate the alteration of time-varying brain activity and topological properties in patients with PS. First, ICA was used to extract independent components (ICs) of each resting-state network (RSN). Second, a sliding window method and k-means clustering were applied to create covariance matrices of all RSNs and to identify recurring whole-brain network connectivity states. Subsequently, graph theoretic analysis was used to evaluate the topological properties of each connectivity state. Third, the temporal features of dFNC (such as the reoccurrence fraction, the dwell time in each state, and the total transition number between all states) were calculated and compared between the PS and HC groups, and the changes of network topological properties (including global efficiency and local efficiency) were evaluated among all four states. Fourth, we explored the correlations between the temporal features that had significant intergroup differences and behavioral assessment scores. We suspected that patients with PS have abnormal temporal features and topological properties compared with HC patients, and that the abnormalities observed by the dFNC method may correlate with functional impairment, improving our understanding of the neurophysiological processes in patients with PS.

Materials and methods

Participants

The sample evaluated in this study consisted of 110 participants from two groups, including 52 patients with chronic pontine stroke (28 left- and 24 right-sided PS) and 58 age-, sex-, and education-matched HC patients. All subjects were recruited at two medical centers, the First Affiliated Hospital of Zhengzhou University and Tianjin Medical University General Hospital. The local medical research ethics committees of the respective hospitals approved the study, and all participants signed informed consent. The inclusion criteria were as follows: (1) first-onset ischemic stroke, (2) single lesion of ischemic infarct involving the pons, (3) right-handedness before stroke, (4) greater than six months since PS onset, and (5) age ranging from 40 to 80 years old. Exclusion criteria were as follows: (1) recurrent stroke, (2) other brain abnormalities in the past, (3) other brain structure damage identified by MRI examination, (4) history of alcohol or drug dependency, (5) MRI contraindications or poor physical conditions that might influence image acquisition, and (6) a score of greater than one on the Fazekas scale(Fazekas et al., 1987) for white matter hyperintensity when assessed via T2 fluid-attenuated inversion recovery (T2 FLAIR) images.

Behavioral measures

Motor function was assessed using the Fugl-Meyer Test (FMT), including upper and lower limb scores. Memory function was estimated in each subject using the Rey Auditory Verbal Learning Test (RAVLT). For RAVLT, the participants listened to 15 words with a presentation rate of one word per second and were requested to recall as many words as possible. The procedure was repeated five times, and the five recall trials were summed into one score (short-term recall score of RAVLT, RAVLT-S). After 20 min, the participants were asked to recall as many of the 15 above-mentioned words as possible (long-term recall score of RAVLT, RAVLT-L).

MR imaging data acquisition

All images of subjects were collected on two same-type scanners (Discovery MR750 3.0 Tesla, General Electric, Milwaukee, WI, USA) using the same parameters at both medical centers. All participants completed a 6-minute resting-state scan, during which they were asked to lie quietly in the scanner with their eyes closed. For each subject, resting-state fMRI data were obtained using a gradient echo single-shot echo-planar imaging sequence with the following parameters: repetition time (TR) / echo time (TE) = 2000 ms / 30 ms, thickness = 3.0 mm, gap = 1.0 mm, field of view (FOV) = 240 mm × 240 mm, matrix = 64 × 64, slices = 38, and volumes = 180. Sagittal three-dimensional T1-weighted images (3D-T1WI) were acquired by a brain volume (BRAVO) sequence with the following parameters: TR/TE = 8.1 ms / 3.1 ms, thickness = 1.0 mm, no gap, FOV = 256 mm × 256 mm, matrix = 256 × 256, slices = 176, and voxel size = 1 mm × 1 mm × 1 mm. T2 FLAIR images were acquired with the following parameters: TR/TE = 8400 ms / 155 ms, thickness = 5.0 mm, FOV = 240 mm × 240 mm, and slices = 21.

Preprocessing of functional magnetic resonance imaging data

Before data preprocessing, we flipped the imaging data with right hemisphere lesions from right to left along the median sagittal line. Thus, the left side corresponded to the ipsilesional (affected) hemisphere, and the right side corresponded to the contralesional (unaffected) hemisphere for all patients. The probability map of the lesion location is shown in Fig. 1. Firstly, the 3D-T1WIs were spatially normalized to a standard Montreal Neurological Institute (MNI) space (3 mm × 3 mm × 3 mm). Then, the lesions were artificially delineated by an experienced neuroradiologist on the normalized 3D-T1WI layer by layer with the MRIcron software (http://www.mccauslandcenter.sc.edu/mricro/mricron/), and this process generated a lesion mask for each patient with PS. Subsequently, the lesion masks of all patients with PS were overlapped to obtain an average individual lesion mask. Finally, the average lesion mask was overlaid on the standard MNI template. All MRI data were preprocessed using the Data Processing & Analysis for Brain Imaging (DPABI; http://rfmri.org/DPABI) and Statistical Parametric Mapping 8 program (SPM8; http://www.fifil.ion.ucl.ac.uk/spm). Firstly, the raw DICOM images of each participant were converted to Neuroimaging Informatics Technology Initiative (NIFTI) format files, and the first 10 time points were removed to ensure steady-state magnetization. Secondly, the remaining 170 volumes were processed with slice timing correction and motion correction to reduce displacement between volumes. Finally, the functional images were spatially normalized to an MNI space (3 mm × 3 mm × 3 mm) using Diffeomorphic Anatomical Registration Through Exponentiated Lie algebra (DARTEL) and smoothed with a Gaussian kernel of 8 mm × 8 mm × 8 mm full-width at half-maximum to increase the signal-to-noise ratio. To eliminate the influence of motion on dFNC, participants with a maximum displacement above 1.5 mm and maximum rotation above 1.5 degrees were excluded before further analysis. In practice, no subjects were excluded.

Fig. 1

The overlapped lesion map of 52 S patients. The color bar indicates the number of subjects with lesions in each voxel. Z axis from Z = 26 to Z = 50 in MNI coordinates, with an incremental interval of 4. The imaging data with right hemisphere lesions were flipped from right to left along the median sagittal line. Abbreviations: L = left; R = right.

Identification of resting-state networks

Pre-processed images were decomposed into different functional networks using the Group ICA of fMRI Toolbox software (GIFT version 4.0b; ), which is based on blind source-separation technology. The ICA analysis was conducted in three stages: (a) data reduction, (b) application of the ICA algorithm, and (c) back-reconstruction for each subject. Firstly, based on all available fMRI data (for both the PS and HC groups), the number of ICs were auto-estimated by the GICA software, yielding 23 components; all subjects’ data then were reduced in dimensionality using two stages of principal component analysis (PCA). Secondly, Group ICA was performed using the infomax algorithm(Bell and Sejnowski, 1995) to extract independent spatial maps and time courses for every component. To improve the stability and consistency of the separation, this analysis was iterated 100 times using the ICASSO(Himberg et al., 2004) algorithm. Lastly, the individual-level components were obtained from back-reconstruction(Erhardt et al., 2011). The intensity values of connectivity within each independent component were represented by the z-scores, which reflect the degree of correlation between the time series of a given voxel and the mean time series of its corresponding component.

Of the 23 ICs, 13 were considered meaningful based on the following criteria: (a) peak coordinates of spatial maps located primarily in gray matter; (b) no spatial overlap with vascular, ventricular, or susceptibility artifacts; (c) time courses dominated by low-frequency signals (ratio of powers below 0.1 Hz to 0.15–0.25 Hz in the spectrum); and (d) time courses characterized by high dynamic range (a range difference between the minimum and maximum power frequencies)(Allen et al., 2014). The locations of these RSNs were consistent with those obtained in previous studies(Smith et al., 2009, Seeley et al., 2007). The spatial maps of these ICs are presented in Fig. 2.

Fig. 2

Spatial maps of the 13 independent network components. They were grouped into 8 RSNs based on their anatomical and functional properties including the auditory network (AUN), the visual network (VIS), the sensorimotor network (SMN), the attention network (ATN), the default mode network (DMN), the executive control network (ECN), the salience network (SAN), and the cerebellar network (CB).

Post-processing of time courses

Following the previous studies(Allen et al., 2014, Tu et al., 2019, Mennigen et al., 2019), we performed additional post-processing steps on time courses of the 13 ICs considered meaningful, including: (a) removing linear, quadratic, and cubic detrends; (b) regressing out six realignment parameters and their temporal derivatives (in x-, y-, and z-directions as well as pitch, roll, and yaw); (c) low-pass filtering of 0.15 Hz; and (d) removing spikes to ensure that artifactual spikes did not negatively impact the signal analysis. Finally, we used the residual time courses to perform dFNC analyses. To create the static functional connectivity matrix, pair-wise Pearson’s correlations were calculated using the post-processed time courses between ICs throughout the entire scanning, and then converted to z values with Fisher’s z-transformation, as shown in Fig. 3.

Fig. 3

Static functional network connectivity matrix of all pairwise correlations between the selected 13 ICs time-courses averaged over subjects. Hot and winter colors indicate positive and negative correlations. Color bar correlation values range from −0.5 to 0.5. Abbreviations: AUN = the auditory network, VIS = the visual network, SMN = the sensorimotor network, ATN = the attention network, DMN = the default mode network, ECN = the executive control network, SAN = the salience network, CB = the cerebellar network.

Controlling for head motion

We adopted the following steps to minimize the potential effects of head motion on dFNC. Firstly, participants with a maximum displacement above 1.5 mm and maximum rotation above 1.5 degrees were excluded. Secondly, the “artifactual” components were identified and removed using ICA, which has been shown to successfully remove motion-related noise from fMRI data(Uddin, 2017). Thirdly, the selected IC time-courses were despiked, the data points with a root mean square of the frame-wise displacement > 0.5 mm were defined as spikes based on the 3Ddespike algorithm as implemented in Analysis of Functional NeuroImages (AFNI)(Cox, 1996), and these points were interpolated with a third-order spline fit to the clean portions of the data(Mennigen et al., 2019, Mennigen et al., 2019). This approach is similar to the “scrubbing” method(Power et al., 2012) with the advantage that it does not eliminate volumes that would disrupt the temporal continuity that is vital for a dFNC analysis(Nomi et al., 2017). Previous research has shown that despiking decreases outlier impact on FC analyses(Allen et al., 2014). Lastly, motion was regressed out of the data during dFNC processing using six realignment parameters. Together, these strategies provide an effective combination of linear and non-linear motion reduction. In addition, we repeated our main analyses after removing subjects with mean framewise displacement exceeding 0.5 mm, thereby ensuring that results were not driven by head motion (please see detail information in Supplementary material 3).

dFNC computation

DFNC analysis was examined using a sliding-window approach in the dFNC toolbox of GIFT. Based on the technique employed in previous studies(Allen et al., 2014, Damaraju et al., 2014), we selected a length (22 TRs) with a Gaussian value of σ = 3 TRs, in a step of 1 TR, resulting in 148 windows. For each window, we estimated the functional connectivities between ICs from a regularized inverse covariance matrix; in this analysis, a penalty on the L1 norm (i.e., the sum of the absolute values of the elements of the precision matrix) was imposed in the graphic LASSO framework with 10 repetitions to promote sparsity in estimation. Thus, for each subject and each window, we obtained a functional connectivity matrix with 78 (13 × 12/2) edges reflecting the time-varying FCs between the 13 ICs. Finally, a Fisher's r-to-z transformation (where r is the Pearson correlation coefficient and z is approximately normally distribution) was applied to all FC matrices to improve the normality of the correlation distribution.

k-means clustering

K-means clustering was used to detect specific FC patterns of the windowed covariance matrices. In this analysis, we use the squared Euclidean distance (500 iterates and 150 replicates) to measure the similarity between different time windows(Malhi et al., 2019). The optimal number of clusters was estimated to be four (k = 4) using the elbow criterion computed as the ratio between within-cluster distance to between-cluster distances(Allen et al., 2014), and the centroids of 4 clusters are shown in Fig. 4.

Fig. 4

The centroid of each functional network connectivity state with the total number of occurrences and percentage of total occurrences.

Dynamic properties analysis

We statistically evaluated the following temporal properties of the dFNC states: (1) the reoccurrence fraction in each state (the percentage of total time the subject spent in each state), (2) dwell time in each state (the time the subject spent in a given state without changing to another one), and (3) the total transition number between states for each subject (the summary of all transition numbers between the 4 connectivity states). A Shapiro-Wilk test was used to verify Gaussianity of three indices for each state and each group, whereafter we tested the significance of differences in the temporal properties between the PS and HC groups using a Mann-Whitney U test (the reoccurrence fraction and dwell time in each state, the total transition number). A p value of < 0.05 [false discovery rate (FDR) corrected] was considered significant.

Graph-theory parameter analysis

We applied a graph theoretic analysis to detect the topologic organization of each dFNC state with GRETNA software (http://www.nitrc.org/projects/gretna), wherein the 13 ICs and inter-IC correlations were defined as nodes and edges, respectively, resulting in a 13 × 13 connectivity matrix for each subject and each state. We applied a sparsity threshold S (the ratio of the number of actual edges to the maximum possible number of edges in a network) to sparsify all connectivity matrices. As the measures of topological properties are affected by the choice of the sparsity value, we set threshold values for sparsity with a range of 0.22 (the minimum sparsity threshold estimated using the code named gretna_get_rmax.m in GRETNA software) to 0.5 in a step of 0.01, in order to balance the prominent small-world attribute and the appropriate sparseness in brain functional networks across subjects(Watts and Strogatz, 1998). Only positive connections were considered; an edge that was larger than the threshold was specified as 1, while an edge that was smaller than the threshold was specified as 0(Zhu et al., 2018). As a result, an undirected and unweighted graph was generated. For the adjacency matrix at each sparsity threshold, we calculated global and local network efficiency of each window and each subject to explore efficiency of information transfer. To avoid the selection of a specific threshold, an area under the curve (AUC) method was applied for the topological measurements; such an approach has been widely used in previous researches(Kim et al., 2017, Li et al., 2018). Finally, we applied a Kruskal-Wallis H-test with state as a factor to assess the differences of network topology in different states.

Validation analyses

To validate the robustness of our results, we used different sliding window lengths and number of clusters to verify the robustness of our results (see detailed information in Supplementary Material 1). The additional sliding window lengths were set at 20 TRs and 30 TRs, and the numbers of clusters were set at 5 and 6.

Statistical analysis

Comparison of clinical scales between groups

The two-sample t-test was used to detect differences in age and education level between the PS and HC groups. A Pearson chi-square test was used to identify differences of gender between the PS and HC groups. For clinical variables, the two-sample t-test was used to detect differences in the scores of RAVLT between the PS and HC groups.

Brain-behavioral relationship

For each temporal feature showing significant intergroup differences (Mann-Whitney U test, p < 0.05, FDR corrected), we performed Spearman’s correlation analysis between the temporal features and the clinical variables in the PS and HC groups, respectively (p < 0.05, FDR corrected). The clinical variables included FMT and RAVLT.

Results

Demographic and clinical parameters

Detailed demographic and clinical data for the PS and HC groups are provided in Table 1. For the demographic data, no significant differences were observed between the PS and HC groups in terms of gender (x = 0.010, p = 0.921, Pearson chi-square test), age (t = 1.615, p = 0.109, two-sample t-test), years of education (t = -1.957, p = 0.053, two-sample t-test). For the clinical data, no significant differences were observed in RAVLT-S (t = -1.032, p = 0.305, two-sample t-test) and RAVLT-L (t = -1.380, p = 0.170, two-sample t-test) between the PS and HC groups. The probability map of the lesion location is shown in Fig. 1.

Table 1

Demographic and clinical data of stroke patients and healthy controls.

Variables	Stroke patients (n = 52)	Healthy controls (n = 58)	x2/t	p-Value
Gender (male/female)	30/22	34/24	0.010	0.921
Mean age (years)	58.19 ± 7.30	55.83 ± 7.98	1.615	0.109
Education (years)	9.62 ± 3.40	10.76 ± 2.72	−1.957	0.053
Lesion side (left/right)	33/29	---		---
Lesion volume (voxels)	112.20 ± 95.50	---		---
FMT-Upper	62.46 ± 8.63	---		---
FMT-Lower	31.25 ± 7.11	---		---
RAVLT-S	46.38 ± 12.65	48.53 ± 8.20	−1.032	0.305
RAVLT-L	10.32 ± 3.29	11.14 ± 2.87	−1.380	0.170

Data are presented as mean ± SD. Abbreviations: FMT = Fugl-Meyer Test; RAVLT = the Rey Auditory Verbal Learning Test; S = the short-term memory score; L = the long-term memory score; SD = standard deviation.

We identified eight resting-state networks (RSNs) from the 23 ICs (Fig. 2), including the auditory network (AUN: IC8), the visual network (VIS: IC4, IC15), the sensorimotor network (SMN: IC2, IC18), the attention network (ATN: IC16, IC22), the default mode networks (DMN: IC7, IC19), the executive control network (ECN: IC5, IC10), the salience network (SAN: IC9), and the cerebellar network (CB: IC11). The static functional network connectivity matrix is shown in Fig. 3, representing all pairwise correlations between the 13 selected IC time-courses averaged over subjects.

Dynamic functional network connectivity states

The dynamic interactions of the eight functional networks were evaluated using the sliding window and k-means clustering analyses. The 148 dFNC matrices were clustered into 4 connectivity states that were recurrent throughout the resting-state fMRI acquisition and in all subjects. The centroids of the four states are shown in Fig. 4. State 1, which accounted for 44% of all windows, was generally characterized by extensively sparse FNCs within and between all RSNs. State 2, which accounted for 15% of all windows, was characterized by a high level of positively local connectivities within primary perceptional domains (AUN, VIS, SMN, ATN) and higher cognitive control domains (DMN, ECN, SAN, CB), and a high level of negative inter-FCs between primary perceptional domains and higher cognitive control domains. State 3, which accounted for 14% of all windows, was characterized by tightly positive FNCs within and between all RSNs, except the existence of negative FNCs in CB related to the other networks. State 4, which accounted for 27% of all windows, was characterized by tightly positive connections between ATN, DMN, and ECN.

Temporal properties of functional connectivity states

Table 2 and Fig. 5 show that the PS group had an increased reoccurrence fraction in State 1 (Mann-Whitney U = 1100.5, p = 0.015; p < 0.05) as well as a decreased reoccurrence fraction and mean dwell time in State 2 (for reoccurrence fraction, Mann-Whitney U = 1115, p = 0.013; for dwell time, Mann-Whitney U = 1108, p = 0.012; p < 0.05) compared to the HC group.

Table 2

The temporal properties for each state and each group.

Variables	Stroke patients (n = 52)	Healthy controls (n = 58)	Mann-Whitney U	Z	p-Value
RF1	0.56(0.23 ~ 0.74)	0.31(0.13 ~ 0.59)	1100.5	−2.44	0.015
RF2	0(0 ~ 0.14)	0.15(0 ~ 0.33)	1115	−2.472	0.013
RF3	0.03(0 ~ 0.23)	0(0 ~ 0.16)	1362.5	−0.924	0.355
RF4	0.1(0 ~ 0.41)	0.21(0.01 ~ 0.52)	1311	−1.192	0.233
DT1	26.17(13.2 ~ 35.46)	17.71(8.75 ~ 30.33)	1213.5	−1.764	0.078
DT2	0(0 ~ 12.44)	11.2(0 ~ 19.96)	1108	−2.516	0.012
DT3	4.5(0 ~ 15.88)	0(0 ~ 14.7)	1419.5	−0.562	0.574
DT4	12.04(0 ~ 22.13)	13.25(2 ~ 26.92)	1364	−0.871	0.384
TTN	5.5(4 ~ 7)	5(4 ~ 7.25)	1452.5	−0.335	0.738

Data are presented as Median(Q1 ~ Q3) and were compared with Mann-Whitney U test (pFDR < 0.05). Abbreviations: RF = reoccurrence fraction; DT = mean dwell time. TTN = total transition number; Q1 = the 25th percentile; Q3 = the 75th percentile

Fig. 5

Temporal properties of dFNC states for the PS and HC groups. (a) Reoccurrence fraction, and (b) dwell time were plotted using violin plots. Horizontal lines indicate group medians and interquartile ranges (solid and dashed lines, respectively). The width of the violin plot bars indicates the distribution density of subjects in each group at the corresponding ordinate level. The horizontal lines above some violins indicate pairwise comparisons that demonstrated statistical significance: * P < 0.05, Mann-Whitney U test. Abbreviations: PS = pontine stroke; HC = healthy controls.

Relationship with clinical properties

In a further analysis of correlations between temporal features and clinical characteristics, we found that the dwell time in State 2 was positively correlated with the RAVLT-L score in the HC group (Fig. 6a, r = 0.323, p = 0.013, FDR corrected), whereas the dwell time in State 2 had no correlation with RAVLT-L (Fig. 6b, r = -0.005, p = 0.974) in patients with PS, indicating a decoupling between duration in State 2 and the memory function in the PS group. No significant link between temporal features and FMT scores was found.

Fig. 6

Correlations between the temporal features that existed intergroup differences and the behavior scores. (a) A significantly positive relationship was observed between the mean dwell time in State 2 and the RAVLT-L score (r = 0.323, p = 0.013, FDR corrected) in the HC group; (b) no significant relationship was observed between the mean dwell time in State 2 and the RAVLT-L score (r = -0.005, p = 0.974) in the PS group. The points where the ordinates are zero correspond to participants not spending time in this state. Abbreviations: RAVLT-L = long-term memory scores of the Rey Auditory Verbal Learning Test; HC = healthy controls; PS = pontine stroke.

Dynamic graph theory properties: Network efficiency

We calculated network topological properties of each window and each subject, and compared global and local efficiencies among states. Fig. 7 shows violin plots of global and local efficiencies in each dFNC state. In general, temporal dynamics (across 4 states) had a significant impact on global efficiency (H = 1306.3, p < 0.001, Kruskal-Wallis H-test) and local efficiency (H = 1826.1, p < 0.001, Kruskal-Wallis H-test; see detailed values in Supplementary Material 3), indicating that functional networks have different information transfer efficiencies in different states of brain connectivity. Analysis by a post hoc test indicated that State 1 had the highest global efficiency, while State 2 had the highest local efficiency in each of the four states.

Fig. 7

The (a) global efficiency and (b) local efficiency in different dynamic states are shown in violin plots. Horizontal lines indicate group medians and interquartile ranges (solid and dashed lines, respectively). The width of the violin plot bars indicates the distribution density of subjects in each state at the corresponding ordinate level. The horizontal lines above the violins indicate pairwise comparisons that demonstrated statistical significance: *** P < 0.001.

Discussion

Our study primarily investigated abnormal time-varying brain activity and network topology properties in patients with PS using the dFNC method. The dFNCs within the whole brain networks could be clustered into four reoccurring states. The patients with PS had an increased reoccurrence fraction in State 1, which showed the highest global efficiency among the four states, as manifested by extensively sparse connections within and between all RSNs. In contrast, patients with PS had a decreased reoccurrence fraction and dwell time in State 2, which showed the highest local efficiency across each of the four states, as manifested by highly segregated connections that showed tightly positive connectivities within primary perceptional domains and within higher cognitive control domains, while tightly negative inter-FCs among primary perceptional and higher cognitive control domains. In addition, a positive connection between the dwell time in State 2 and the RAVLT-L scores (evaluating participants’ memory function) was found in the HC group, while a decoupling between duration in State 2 and the memory function was observed in the PS group. These findings suggested that the patients with PS are characterized by a higher frequency of a state with sparse connections and a lower frequency of a state with segregated connections, indicating an imbalance between functional segregation and integration in the whole brain after pontine infarction.

Abnormal temporal features and topology properties in PS

Increased reoccurrence fraction in state 1

State 1 was characterized by extensively sparse FNCs within and between all RSNs; hence, we refer to State 1 as a ‘sparsely connected state’. Previous studies revealed that such weak connectivity may represent the average of a wide variety of additional states that are insufficiently distinct or frequent to be separated(Allen et al., 2014, Marusak et al., 2017). Furthermore, recent studies suggested that this state may be associated with self-referential processing and drowsiness(Allen et al., 2018). Though there is a relatively high correlation within the DMN in this state, researchers have revealed that increased connectivity within the DMN is related to self-referential processing(Monsa et al., 2018), and that a high internal focus of attention may result in an increase in DMN connectivity(Gusnard et al., 2001). Considering these results together, we speculate that this dynamic state may be related to more interoceptive awareness and reduced vigilance. In addition, we found that State 1 showed the highest global efficiency in all four states. The whole brain can be divided into different brain areas, which form a complex network and can efficiently transfer information from one region to the other, a process that can be called integration. Global efficiency is the average efficiency across all node pairs, indicating how efficiently information is integrated across the entire network(Kawagoe et al., 2017). Greater global efficiency means a higher degree of network integration and may reflect an advantage in cognitive function in healthy individuals(Bassett and Bullmore, 2009, Bassett et al., 2009, van den Heuvel et al., 2009). However, one study of older participants revealed that increased global efficiency was related to a slight decline in memory performance(Stanley et al., 2015). Older individuals constituted the majority of the subjects we selected; and the sparsely connected State 1 may represent widespread damage of the entire brain network; hence, we hypothesize that the increased occurrence of State 1 may imply impairment of the memory function in patients with PS. Previous studies have revealed that patients with PS existed impairment of both the ipsilesional and contralesional rubral branches(Guo et al., 2019), as well as reduction of gray matter volume(Jiang et al., 2017) and cerebral blood flow(Wang et al., 2019) in the cerebellum. In addition to these observations; the results of our study may provide a new clue to the functional impairment in patients with PS; in terms of whole-brain functional network connectivity. The same finding that patients had a longer duration in the sparsely connected state also has been reported in multiple brain diseases, such as depressive disorders(Zhi et al., 2018, Wang et al., 2020), schizophrenia(Yu et al., 2015, Rabany et al., 2019, Espinoza et al., 2019), Huntington's disease(Espinoza et al., 2019, Espinoza et al., 2019); and neurofibromatosis(Mennigen et al., 2019, Mennigen et al., 2019). Given the above research, such abnormalities in dFNC may be the underlying mechanism of similar dysfunctions in different brain disorders.

Decreased reoccurrence fraction and dwell time in state 2

State 2 was characterized by highly segregated connections, which showed tightly positive connectivities within primary perceptional domains (AUN, VIS, SMN, ATN) as well as higher cognitive control domains (DMN, ECN, SAN, CB), while tightly negative inter-FCs between primary perceptional and higher cognitive control domains; therefore, we called State 2 a ‘highly segregated state’. We found that State 2 had the highest local efficiency among the four states in the graph theory analysis. Local efficiency refers to the flow of information in a local network, which is a part of the whole brain network and provides a perspective on functional segregation(Smitha et al., 2017). Segregation is associated with brain modularity; higher segregation of domains means greater modularity(Bonkhoff et al., 2020). Previous research has revealed that modularity is associated with cognitive function in healthy people; including learning ability(Bassett et al., 2011); amenability to cognitive training(Gallen et al., 2016), and normal brain aging(Meunier et al., 2009). Abnormal modular organization also may occur in patients with neuropsychiatric disorders, such as Alzheimer's disease(Brier et al., 2014, Buckner et al., 2009), schizophrenia(Lynall et al., 2010, Yang et al., 2016), and traumatic brain injury and stroke(Siegel et al., 2018, Gratton et al., 2012, Caeyenberghs et al., 2017). Siegel et al.(Siegel et al., 2018) found that increased modularity was related to better recovery in patients after stroke. Our study observed decreased reoccurrence fraction and dwell time in State 2 in the PS group. Additionally, the dwell time in State 2 was positively correlated with the RAVLT-L in the HC group, whereas there was a decoupling between the duration in State 2 and the memory function in the PS group. Therefore, we speculate that the decrease in such a highly segregated state would lead to disruption of functional segregation and modular organization, further resulting in poor memory recovery in patients. An intriguing finding is that the difference in RAVLT-L between the PS and HC groups was not obvious. The reason may be that for most of the enrolled patients, the scan time interval after stroke was just over 6 months. The patients may exist a potential risk of memory impairment according to our findings, though their memory function did not decline at the time of behavioral assessment. We intend to monitor these patients’ changes of memory function in the future.

In fact, the balance between integration and separation represents a small-world property that is an optimal topology for brain networks, supporting specialized processing within tightly connected regions and facilitating efficient information transfer between disparate areas(Johnson et al., 2020). In reviewing previous studies of cerebrovascular diseases, researchers have observed abnormalities of integration and separation using dFNC analysis in patients with subcortical ischemic vascular disease (SIVD)(Fu et al., 2019) and acute ischemic stroke(Bonkhoff et al., 2020). Fu et al.(Fu et al., 2019) found that patients with SIVD spent more time in a weak and diffused dynamic state, but spent less time in a highly segregated state; the highly segregated state was manifested by tightly positive correlations within primary perceptional domains as well as higher cognitive control domains, while negative connections between the primary perceptional domains and higher cognitive control domains. Our results were similar to those of Fu et al., possibly indicating that the alteration of separation and integration are relatively stable in patients with cerebral infarction. In other work, depending on the degree of motor impairment after acute ischaemic stroke, Bonkhoff et al.(Bonkhoff et al., 2020) observed that moderately affected patients showed more time in a state with weak connections, while severely affected patients transitioned more often to a state with regionally dense connections, characterized by strongly positive intra-domain connectivity and negative inter-domain connectivity. The changes in moderately affected patients are much alike part of our findings, but the change in severely affected patients is the opposite of what we found. There are two possible reasons for this phenomenon. One is that the infarct lesions differed between the patients enrolled in the two studies. Previous studies have confirmed different alteration patterns of cerebral blood flow (CBF) and gray matter volume (GMV) in stroke patients based on the infarct location(Jiang et al., 2017, Wang et al., 2019), so we suspect that there may be also some differences in dFNC between PS and strokes in other locations. The other possible reason is that the enrolled patients were at different infarction stages. Differences in regional temporal variability among the acute, subacute, and chronic stages have been reported in stroke patients using dynamic network analysis(Hu et al., 2018). Future studies will need to investigate whether there are different organizational structures of dFNC among PS patients with infarction stages. In summary, the abnormalities of integration and separation were found both in previous studies and ours. Therefore, we suspect that abnormal duration in a sparsely connected state or a highly segregated state may be a common feature of cerebral ischemia dysfunction, which might serve as a specific diagnostic marker for the preliminary diagnosis of cerebrovascular disease.

In the validation analyses, we found that our results were well replicated when using various analytical parameters, including the sliding window lengths and number of clusters. These observations indicated that the abnormal temporal features and network topological properties in patients with PS are stable, meaning that these parameters can be used as neuroimaging markers to distinguish patients with PS from healthy people. Previous researchers also have verified results obtained using dFNC, confirming the reliability and repeatability of dFNC analysis(Shi et al., 2018, Wu et al., 2019). One interesting finding is that patients with PS had decreased total transition numbers compared with members of the HC group when the number of clusters was set at 5 or 6, while there was no significant difference in the total transition number between the two groups in the main study (k = 4). This observation suggested that, as the number of clusters increases, the brain's variability may become more readily apparent. Greater functional neural variability may lead to greater cognitive and behavioral flexibility(Marusak et al., 2017, Hutchison and Morton, 2015, McIntosh et al., 2008), and the decline in this flexibility may lead to a reduced ability for information processing.

Limitations

There are several limitations to our study. Firstly, we did not investigate the specific differences in dFNC between PS and strokes in other locations, making it unclear whether these abnormalities are common to all strokes or specific to PS. Secondly, only right-handed patients were recruited in this study; further work will be need to validate our results in left-handed patients. Lastly, the behavioral scales used in this study were limited; more scales that are sensitive to neuroimaging measures in patients with PS should be included to confirm our findings.

Conclusion

In this study, we explored alterations of the temporal features and network topological properties in PS based on dFNC analyses, providing a new perspective for the pathological process of PS. Patients with PS were found to have a longer duration in a sparsely connected state with the highest global efficiency among all reoccurring states, which might be associated with a slight decline in memory performance. In contrast, patients with PS had a shorter duration in a highly segregated state with the highest local efficiency among all reoccurring states, which might lead to poor recovery of memory function. In summary, the use of the dFNC method provided new evidence for aberrant time-varying brain activity and network properties of PS, indicating abnormal functional segregation and integration of the whole brain and confirming the imperfection of functional network connectivity in patients with PS. These findings might facilitate our understanding of the dynamic neural mechanism underlying functional impairment in PS.

CRediT authorship contribution statement

Yingying Wang: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Writing - original draft, Visualization. Caihong Wang: Conceptualization, Methodology, Formal analysis, Investigation, Funding acquisition. Peifang Miao: Methodology, Investigation. Jingchun Liu: Investigation, Resources. Ying Wei: Investigation, Resources. Luobing Wu: Investigation, Resources. Kaiyu Wang: Methodology, Software. Jingliang Cheng: Conceptualization, Resources, Data curation, Supervision, Project administration, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.