Integration of temporal & spatial properties of dynamic functional connectivity based on two-directional two-dimensional principal component analysis for disease analysis

Subjects

The rs-fMRI data of 45 autism spectrum disorders (ASD) patients and 47 healthy controls (HC) groups were scanned using a 3-T Siemens Allegra scanner at the NYU Langone Medical Center. Due to the difference in medical device, collection protocol, etc to mitigate data heterogeneity, scans were exclusively considered from 45 individuals diagnosed with ASD and 47 individuals within the age range of 7 to 15 years. The imaging was conducted at NYU Langone Medical Center. All subjects under consideration exhibited minimal head movement, with displacement <1.5 mm or angular rotation <1.5° in any of the three directions. The ASD subjects were diagnosed according to the autism criteria in the Diagnostic and Statistical Manual of Mental Disorders (4th Edition, Text Revised) (DSM-IV-TR) (American Psychiatric Association, 2000). For more details of the data collection, exclusion criteria, and scan parameters, please visit the ABIDE website. The main scanning parameters used in this dataset include the flip angle = 90, 33 slices, TR/TE = 2,000/15 ms, 180 volumes, and voxel thickness = 4 mm.

The rs-fMRI data of 52 end-stage renal disease (ESRD) patients and 49 HC groups were scanned using GE Signa HDX 3.0T MRI equipment at hospital of Qingdao University. Inclusion criteria for the ESRD group: (1) Diagnosed by the Affiliated Hospital of Qingdao University as Chronic Kidney Disease Stage 5, glomerular filtration rate less than 15 ml $\cdot {min}^{-1}\cdot {\left(1.73m^2\right)}^{-1}$. (2) Age 18–70 years old; (3) No contraindications to magnetic resonance examination and claustrophobia. Exclusion criteria: (1) history of severe traumatic brain injury; (2) intracranial organic pathology, intracranial organic disease, such as tumour, infarction, haemorrhage, etc.; (3) cerebrovascular disease, such as cerebral arteriovenous malformation, smoky disease, etc. and smog disease; (4) history of mental illness; and (5) history of substance abuse (drug ethanol or cigarettes); (6) hypertension, diabetes mellitus, coronary artery disease, heart failure, liver and kidney failure, etc. failure, liver or kidney failure, and other serious systemic diseases; (7) incompatibility with the test subjects who are incompatible with the test or unable to effectively complete the MRI scan. The scan parameters are TR = 3,000 ms, TE = 40 ms, FA = 90°, 25 slices, thickness = 5 mm, gap = 0 mm, matrix size = 96 ×96 and FOV = 24 cm × 24 cm.

Data preprocessing

Data preprocessing is performed for all rs-fMRI data by using a standard pipeline, including time correction, head motion correction, spatial normalization and smoothing using MATLAB’s DPABI toolbox. Specifically, (1) conversion of data format; (2) removal of data from the first 10 sampling points for each subject in order to exclude the effects of initial magnetic field instability and the subject’s initial emotional instability; (3) correction of time layer, head motion. Subjects with excessive head movements are excluded; (4) spatial normalization in terms of MNI standard space with the resolution of 3 × 3 × 3 mm²; (5) ventricle signals, cerebrospinal fluid signal and white matter signal were regressed out as nuisance signals; (6) to filter out various physiological noises to reduce the effects of low-frequency linear drift and physiological noise, the data were processed to remove linear trends and band-pass filtered (0.01−0.08 Hz).

In this experiment, after data pre-processing, the brain was divided into 116 brain regions for each subject according to the AAL template, and the average of the BOLD signal of all voxels in each brain region was calculated to obtain the average time series of the 116 brain regions.

Construction of DFC

For each subject, let ${\mathrm{z}}_i=\left({\mathrm{z}}_{i1},{\mathrm{z}}_{i2},\dots ,{\mathrm{z}}_{i\mathrm{N}}\right)$ (i = 1, 2, …,116) denotes the average rs-fMRI time series of the ith brain region, where N denotes the number of image time points. Dynamic FC is implemented using a sliding window approach, assuming that the length of the sliding window is T, and the step length between every two adjacent windows is S. Then the time series of length N is divided into K partially overlapping subseries with each other, where $\mathrm{K}=\frac{\left[\left(\mathrm{N}-\mathrm{T}\right)\right]}{\mathrm{S}}+1$.

Letting ${\mathrm{z}}_i\left(k\right)=\left[{\mathrm{z}}_{i1}\left(k\right),{\mathrm{z}}_{i2}\left(k\right),\dots ,{\mathrm{z}}_{i\mathrm{T}}\left(k\right)\right]\left(k=1,2,\dots ,\mathrm{K}\right)$ denotes the kth time subseries of z_i, then just need to list all sliding windows of FC. FC is often modeled as a FC network (FCN), with each brain ROI as a node in the network, and the strength of FC between a pair of brain ROIs as an edge. The FCN of the kth time subseries can be generated by a symmetric submatrix $\mathrm{D}\left(\mathrm{k}\right)$, defined as: (1)${\displaystyle D\left(k\right)={\left[{\rho }_{ij}\left(k\right)\right]}_{1\le i,j\le K}}$ where ${\rho }_{ij}\left(k\right)$ denotes the Pearson’s correlation between the average time subseries $z_i\left(k\right)$ and $z_j\left(k\right)$, defined as: (2)${\displaystyle {\rho }_{ij}\left(k\right)=\frac{\sum _{t=1}^T\left(z_{it}\left(k\right)-\overline{z_i\left(k\right)}\right)\left(z_{jt}\left(k\right)-\overline{z_j\left(k\right)}\right)}{\sqrt{\sum _{t=1}^T{\left(z_{it}\left(k\right)-\overline{z_i\left(k\right)}\right)}^2}\sqrt{\sum _{t=1}^T{\left(z_{jt}\left(k\right)-\overline{z_j\left(k\right)}\right)}^2}}}$ where $\overline{z_i\left(k\right)}$ and $\overline{z_j\left(k\right)}$ denote average value of the average time subseries $z_i\left(k\right)$ and $z_j\left(k\right)$. The submatrix series ${\left\{D\left(k\right)\right\}}_{k=1}^K$ can describe the temporal change of the connectivity strength for all ROI pairs.

The Fourier transform

Discrete Fourier Transform, as a mathematical tools, is commonly used in signal processing. The Discrete Fourier Transform can convert a signal into a linear combination of a set of sinusoidal functions of different frequencies, which have the property of being easy to implement and observe.

The Discrete Fourier Transform of vector $x\left(k\right),k=\left\{1,2,\dots ,K\right\}$ is defined as: (3)${\displaystyle f\left(u\right)=\sum _{k=1}^{K}x\left(k\right)W_K^{\left(k-1\right)\left(u-1\right)}}$ where u denote the frequency-domain pixel coordinate of C_r, u = {1, 2, …, K}, W_K is defined as: (4)${\displaystyle W_K=e^{\frac{-2\pi i}{K}}.}$

(2D)² PCA

${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ is a matrix feature extraction technique based on two-dimensional principal component analysis (2DPCA). 2DPCA can extract whole-matrix features by finding optimal projection directions from the covariance matrix. First introduce how 2DPCA learns an optimal projection matrix reflecting the information between rows of matrixes. Letting A ∈ R^m×n be a matrix, X ∈ R^n×d be a projection matrix with n projection directions. Projecting A onto X to obtain the matrix Y = AX. Suppose that there are M training matrixes, denoted by $A_k\left(k=1,2,\dots ,M\right)$, and denote the covariance matrix as: (5)${\displaystyle X=\frac{1}{M}\sum _k^M{\left(A_k-\overline{A}\right)}^T\left(A_k-\overline{A}\right)}$ where $\overline{A}$ be denoted as: (6)${\displaystyle \overline{A}=\frac{1}{M}\sum _k^M{A}_k}$ The optimal value for the projection matrix X is composed by the vectors X₁, X₂, …, X_d of G corresponding to the d largest eigenvalues. The value of d can be controlled as follows: (7)${\displaystyle \frac{\sum _{i=1}^d{\lambda }_i}{\sum _{i=1}^n{\lambda }_i}\ge \theta }$ where $\lambda =\left\{{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n\right\}$ is a sequence of eigenvalues in ascending order for G and θ is a pre-set threshold.

2DPCA also can learn an optimal matrix Z reflecting the information between columns, and then projects A onto Z to obtain the matrix B = Z^TA, and denote the covariance matrix as: (8)${\displaystyle Z=\frac{1}{M}\sum _k^M\left(A_k-\overline{A}\right){\left(A_k-\overline{A}\right)}^T}$ (2D)²PCA is a way to simultaneously use the projection matrices X and Z. Projecting A onto X and Z, yielding a matrix C: (9)${\displaystyle C=Z^{T}AX}$ The matrix C in Eq. (9) is the matrix after projection.

Comparison methods

In this section we describe two comparison methods used in subsequent experiments: the central moment method and the topological index method.

Construction of the proposed framework

In this subsection, we describe the proposed framework building process in detail. Figure 1 depicts the main steps of the feature extraction framework, and it can be seen that there are four main steps: (1) Construction of the dynamic network. (2) Construction of the new network. (3) (2D)²PCA learns the best projection matrix. (4) Projecting the new network and selecting features.

Specifically, we first construct the DFC ${\left\{D\left(k\right)\right\}}_{k=1}^K$ on the continuous rs-fMRI time series by the sliding-window strategy. Then in order to retain both temporal and spatial dynamic information, we transform the 3D DFC into a 2D new network, so that the spatial and temporal information to be stored in the row and column structure of the new network.

Figure 2 illustrates the process of new network construction. For each subject, we transform the submatrix series ${\left\{D\left(k\right)\right\}}_{k=1}^K$ into a matrix A ∈ R^M∗K, where M denotes the number of upper triangular elements of matrix $D\left(k\right)$. The upper triangular elements of matrix $D\left(k\right)$ are pulled into a column vector as the kth column of matrix A. Each row of matrix A reflects the change of each functional connectivity series ${\left\{{\rho }_{ij}\left(k\right)\right\}}_{k=1}^K$ over time. To solve the problem that the correspondence of dynamic FC subnetworks from the same window of different subjects cannot be established, we convert the time series ${\left\{{\rho }_{ij}\left(k\right)\right\}}_{k=1}^K$ from the time domain to the frequency domain by Fourier transform, which is calculated by Eq. (3). After transformation, the transformation of functional connectivity with time is represented by a linear combination of ‘trigonometric functions of different frequencies, so we can analyze the problem from a new perspective, which eliminates the problem of time not corresponding.

Furthermore, we use ${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ to learn the best projection matrix, which is calculated by Eqs. (5) and (8). The new network is projected by Eq. (9) with the optimal projection matrix to select features that integrate spatial and temporal properties of the new network. ${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ is widely used for image matrix dimensionality reduction in face recognition and palmprint recognition, which can learn two optimal matrixes from a set of training matrixes reflecting information between rows and between columns of training matrixes, so ${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ can integrate properties in both directions.

Finally, the features extracted by ${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ will have high dimensionality, so an effective feature selection method is conducive to accurate diagnosis. For dichotomous classification tasks, t-test can be effective in selecting significantly different features between subjects and healthy controls. Therefore, we use a t-test to select the features with p-values less than a certain threshold value.

Methods for comparison

We first compare our proposed method with the central moment method (denoted as Central). In the central moment method, the Pearson correlation between ROI based on the time series of each window is first computed, then the brain networks between ROI are constructed, and the multistep central moment of the brain networks of all the windows solved are used as features for the classification task. Then we compare our proposed method with the method that use topological index of functional connectivity to classify (defined as T-index). In the topological index method, the topology is generated by retaining as nodes the functional connections in the dynamic functions that are greater than a certain threshold, and then the topological index of the topology is derived for binary classification.

Moreover, we compare our proposed method with two methods for extracting unidirectional information using two-dimensional principal component analysis (2DPCA) (Wang, Wang & Yan, 2018; Wee et al., 2016; Xu et al., 2019), which include (1) the method using 2DPCA in the temporal direction (denoted as Fourier-temporal), and (2) the method using 2DPCA in the spatial direction (denoted as Fourier-spatial). In this way, it is demonstrated that ${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ can integrate temporal and spatial information improves the performance of classification. It is worth noting that both methods use the Fourier transform to eliminate the effects of time mismatch.

In order to demonstrate that the Fourier transform can play a role in eliminating the time mismatch problem. We also compare our proposed method with that of ${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ which does not use Fourier transform (denoted as Temporal-spatial). All the methods use the same T-test as the proposed method for feature selection, and then the features are fed into a SVM classifier for classification.

Experimental setting

In the experiment, we constructed two classification tasks, including (1) the classification task for autism spectrum disorder (ASD) and normal control group (NC), and (2) the classification task for end-stage renal disease (ESRD) and NC. We used six evaluation indexes, namely, classification accuracy (ACC), true positive rate (TPR), true negative rate (TNR), positive predictive value (PPV), negative predictive value (NPV), and F1 score, to comprehensively evaluate the advantages and disadvantages of the parameters and the classification effect of the methods. Considering that the number of samples used in the experiments is more than 50 and the parameters such as sliding window length, sliding step size, information content thresholds for temporal and spatial directions, and p-value thresholds need to be adjusted in the feature extraction framework proposed in this article. After referring to the relevant literature (Shi et al., 2019), we selected the ten-fold cross-validation method, which appears in the literature, and did not choose the leave-one-out method, which is more computationally intensive. To avoid possible errors in the experiments, the whole ten-fold cross-validation was repeated 10 times, dividing the samples randomly each time. Finally, the average of the 10 times of recognition accuracy was used as the final accuracy.

In the construction of the network, in order to avoid erroneous experimental results with arbitrarily determined window width and step size, we constructed brain networks with different window widths and step sizes, window width parameter T ∈[20, 30, 40, 50, 60] and sliding step size S ∈[2, 4, 6, 8, 10]. Studies have shown that window width and step size in this range produce robust results in image acquisitions and topological properties of brain networks. We take the DFC derived from all the different parameters and take the average value as the final DFC. In the feature selection phase, we choose features with statistical significance p < 0.01. After that, the features are fed into a SVM classifier for training and classification, and the six metrics mentioned in the previous paragraph are calculated (Yang et al., 2022; Yao, Becker & Kendrick, 2021).

Classification performance

Table 1 summarizes the classification performance of each of the six methods in the two classification tasks. From Table 1, we can see that our proposed method performs better than comparison methods. For example, the proposed method acc achieves 88.1% and 76.1% in the two classification tasks and the highest accuracy of the comparison method is at 86.1% and 72.8%. Table 1 shows that the method integrating temporal and spatial information of the network performs better than the method using one attribute alone, which implies that the temporal and spatial information are two complementary pieces of information, so integrating the brain through ${\left(2\mathrm{D}\right)}^2\mathrm{PCA}$ network’s two kinds of information can further improve the classification performance.

In addition, we could observe that, from Table 1, that using Fourier transform to transform the DFC from the time domain to the frequency domain has a higher classification accuracy than the method that does not use Fourier transform, so we can think that Fourier transform can eliminate the effect of time mismatch to a certain extent and bring better classification performance.

In the process of constructing the dynamic brain network, we used the “sliding window” strategy, so there are two hyperparameters, window width and step size, and we subsequently evaluated the effect of hyperparameters on the classification performance in the classification experiments of nephropathic and normal human beings. Figure 3 shows the effect of different hyperparameters on the classification performance in the four methods, namely Fourier-spatial, Fourier-temporal, Central and our proposed method. It can be concluded that there is an effect of hyperparameters on the classification performance, and our proposed method has more stable and higher classification performance on different combinations of hyperparameters, and has the highest classification accuracy of 88.1% when T = 50 and S = 8.