Topological analysis of functional connectivity in Parkinson’s disease

Parkinson’s disease (PD) is a clinically heterogeneous disorder, which mainly affects patients’ motor and non-motor function. Functional connectivity was preliminary explored and studied through resting state functional magnetic resonance imaging (rsfMRI). Through the topological analysis of 54 PD scans and 31 age-matched normal controls (NC) in the Neurocon dataset, leveraging on rsfMRI data, the brain functional connection and the Vietoris-Rips (VR) complex were constructed. The barcodes of the complex were calculated to reflect the changes of functional connectivity neural circuits (FCNC) in brain network. The 0-dimensional Betti number β0 means the number of connected branches in VR complex. The average number of connected branches in PD group was greater than that in NC group when the threshold δ ≤ 0.7. Two-sample Mann–Whitney U test and false discovery rate (FDR) correction were used for statistical analysis to investigate the FCNC changes between PD and NC groups. In PD group, under threshold of 0.7, the number of FCNC involved was significantly differences and these brain regions include the Cuneus_R, Lingual_R, Fusiform_R and Heschl_R. There are also significant differences in brain regions in the Frontal_Inf_Orb_R and Pallidum_R, when the threshold increased to 0.8 and 0.9 (p < 0.05). In addition, when the length of FCNC was medium, there was a significant statistical difference between the PD group and the NC group in the Neurocon dataset and the Parkinson’s Progression Markers Initiative (PPMI) dataset. Topological analysis based on rsfMRI data may provide comprehensive information about the changes of FCNC and may provide an alternative for clinical differential diagnosis.


Introduction
Parkinson's disease (PD) is a clinically degenerative disorder disease of the nervous system with motor and non-motor symptoms.The disease mainly affects people's motor function, such as bradykinesia, tremor, muscle stiffness or rigidity, abnormal walking gait, etc.In addition, it will also affect non-motor function, such as cognitive impairment, insomnia, depression, autonomic nerve dysfunction, and so on.The cause of PD is unclear.In the early premotor stages, the diagnosis of this disease is still difficult (Aarsland et al., 2017).Anatomical magnetic resonance imaging could not detect the loss of dopamine neurons.With the development of molecular biology, neural structure, and functional imaging technology, more and more biomarkers of PD can be discovered, providing the possibility for early diagnosis, disease monitoring, and differential diagnosis, thus achieving accurate early intervention and efficacy evaluation of the disease.Dopamine transporter single photon emission computed tomography (DAT-SPECT) and resting state functional magnetic resonance imaging (rsfMRI) are potential techniques for detecting the survival status of neurons in PD (Wang et al., 2012).Among them, DAT-SPECT can provide quantitative information about dopamine neurons, which is very useful for assessing disease severity and monitoring treatment effectiveness.However, it is expensive and carries the risk of radioactive tracer.On the other hand, rsfMRI is a non-invasive imaging technique that is relatively inexpensive and can avoid ionizing radiation.It can provide information about brain activity, which is very useful for studying the neural mechanisms of PD and evaluating treatment effectiveness.
Therefore, the study of functional connectivity (FC) based on rsfMRI is a promising method.Based on the correlation of time series, FC was preliminary explored and studied through rsfMRI (Aarsland et al., 2017).With rsfMRI, FC can be used to detect a variety of diseases such as Alzheimer's disease, PD, schizophrenia, and so on (Lee et al., 2013).Some studies have shown that motor and cognitive impairment in PD are related to abnormal functional connections (Putcha et al., 2015;Zhang et al., 2015).Graph theory-based approaches used in PD research have shown that topological properties of brain networks are disrupted, which can help identify this type of disease (Luo et al., 2015).Specifically, abnormal local and global network efficiency changes suggest clinical phenomena in PD.The above methods assume that the functional network is stable and does not change over time.A dynamic FC based method was used in recent research.For example, previous studies used rsfMRI and sliding windows to assess differences in dynamic connectivity between normal control (NC) and PD (Kim et al., 2017).Huang et al. (2020) proposed a regression method to model the dynamic correlation matrices as a linear combination of symmetric positive definite matrix to smooth the image acquisition and physiological noise.Pang et al. (2021) distinguished different motor subtypes of PD based on multilevel indices of rsfMRI and Support Vector Machine (SVM).
However, the graph theory analysis method can not describe the characteristics of higher-level complex brain networks.To study the topological characteristics of complex brain networks on a larger scale, researchers began to study Vietoris-Rips (VR) complex filtration based on persistent homology in brain networks.In topological data analysis (TDA), persistent homology is an effective tool to explore the nonlinear structure of the data.Compared with the common methods such as principal component analysis (PCA), cluster analysis, and graph theory (Li et al., 2009), TDA can effectively capture the topological information of high-dimensional data space.This kind of algorithm adopts a free threshold and solves the problem of threshold selection.It measures the topological characteristics of brain network under all possible thresholds (Choi et al., 2014;Chung et al., 2015;Lee et al., 2017).These approaches mainly associate the 0-dimensional Betti numbers β 0 with current varying thresholds.A connected component-based aggregation cost model called Integrated Persistence Features (IPF) was proposed in previous research (Kuang et al., 2019).Different from the above persistent homology feature which based on 0-dimensional Betti numbers, this paper proposes a persistent homology feature based on the 1-dimensional Betti number β 1 .To our knowledge, there is little literature investigating the 1-dimensional Betti number β 1 in PD.Our main contributions of this paper are as follows: (1) The VR complex filtering model was established based on the relationship matrix of the human brain network.The persistent homology method was used to calculate VR filtered barcodes.And then the functional connectivity neural circuits (FCNC) at different thresholds are calculated from the barcodes.(2) Two-sample Mann-Whitney U test and FDR correction are used for statistical analysis to investigate the FCNC changes between PD and NC groups.(3) Through the statistical tests on the number of FCNC in PD and NC groups, our results show that there is a significant statistical difference between the PD group and the NC group.

Basic concepts about persistent homology
The common method to reduce the dimension of data is PCA, but this method will lose some potentially valuable data more or less.Persistent homology provides us a method to find a complete data ship without dimensional reduction.
Persistent homology is an effective tool to analyze highdimensional data and explore the nonlinear structure of data.It can calculate topological features at different spatial resolutions.By identifying persistent topological features over changing scales, persistent homology provides clues for effective analysis of multi-scale networks.Its core idea is to analyze the birth and death of holes in various dimensions in a multi-scale range.To extract persistent homology features, we first need to construct VR complex.Let d(•, •) denote the distance between two points in the metric space Z.The value of δ denotes the threshold.When we change the threshold, we obtain a sequence of complexes.The VR(Z, δ) complex is defined as follows (Silva and Carlsson, 2004): 1.For vertices a and b, edge Betti intervals help describe how the homology of VR(Z, δ) changes with δ.A k-dimensional Betti interval, with endpoints (δ start , δ end ), corresponds roughly to a k-dimensional hole that appears at the threshold δ start , remains open for δ start ≤ δ < δ end , and closes at δ end .
The rank of the homology group is called Betti number (Woo et al., 2014), which is a set of important topological invariants.It uses the connectivity based on k-dimensional simplex complex to distinguish the topological space, which can well reflect the topological structure of an object.The k-dimensional Betti number is the rank of the k-th homology group and represents the number of "holes" in the k-th dimension.For example, the 0-dimensional Betti number β 0 refers to the number of connected branches (Lee et al., 2012).Similarly, the one-dimensional Betti number β 1 intuitively represents the number of one-dimensional "holes." circle centered on these points and radius of the threshold δ/2.In the process of increasing the threshold δ from 0 to the maximum, the 0-dimensional and 1-dimensional Betti numbers change constantly.When δ = 1.4, the number of connected branches is 10, that means β 0 = 10 and the number of 1-dimensional "holes" β 1 = 0 (Figure 1A).In Figure 1B, when δ = 5.2, β 0 = 1 and β 1 = 0.As the threshold δ increase, when δ = 6.3 (Figure 1C), β 0 = 1 and β 1 = 1.When δ = 8.5, some holes are "filled." At this time, there is only one connected branches (β 0 = 1) and one 1-dimensional "holes" (β 1 = 0) (Figure 1D).With the change of the threshold δ, the topological characteristics of VR complex change.This process can be represented by a barcode and a persistence diagram (Figure 2).As shown in Figure 2A, a barcode is a set of finite intervals.Each interval represents the birth and death of holes in the corresponding dimension, and these intervals are parallel to each other.K-dimensional barcodes (k = 0, 1 in Figure 2A) show us the duration of k-dimensional topological features.Generally, we regard the features with very short duration as noise, and the features with long duration as real signal features.In Figure 2B, the persistence diagram provides a multi-scale feature description.The abscissa of each point in the diagram represents the birth of the topological feature, while the ordinate represents its demise.Points away from the diagonal represent features with a long life cycle, while points close to the diagonal represent features with a short life cycle.Among them, the feature that can be maintained for a longer time is a useful persistence feature with stronger robustness.The feature with short life is more likely to be noise or detail.The red vertical line on the ordinate axis is regarded as a point at infinity, representing a topological feature that will never die.

Dataset and preprocessing
The datasets used in this article are the Neurocon dataset (Badea et al., 2017) and Parkinson's Progression Markers Initiative (PPMI) dataset (Marek et al., 2011).The Neurocon dataset includes rsfMRI data from 27 PD patients and 16 age matched NC patients, with each subject undergoing 2 repeated scans.One NC scan was subsequently excluded owing to data corruption.Finally, 54 PD scans and 31 NC scans were included in the final analyses.The Neurocon study has been approved by the Ethics Committee of the Emergency Hospital of  The Neurocon dataset is preprocessed using DPABI (DPABI Software Library v5.1) (Yan et al., 2016) as follows: 1. Removing the first 10 time points of rsfMRI data to obtain a stable signal; 2. Time correction is performed on each slice to ensure that the data on each slice corresponds to the same time; 3. Realign: eliminate the data with the maximum value of translation greater than 3 mm and the maximum value of rotation greater than 3°; 4. Registration of structural image to functional image space.Use anatomical T1 images to register to the standard Montreal Institute of Neurology MNI152 template; 5. Check the coverage of function image by Automask; 6. Bandpass filtering with a frequency of 0.01 ~ 0.1 Hz is set; 7. Normalize using EPI templates; 8. Extract region of interest (ROI) time courses using Anatomical Automatic Labeling (AAL) atlas.
The brain parcellation used in this paper is the AAL brain template (Tzourio-Mazoyer et al., 2002).There are 116 regions in the AAL template, 90 belonging to the brain, and the remaining 26 belong to the cerebellar structure.There are few studies on the cerebellum.In this paper, the classical AAL template of 90 brain regions is used to construct VR complex.

Constructing VR complex
Firstly, we preprocess the rsfMRI image as described above, and then calculate the FC between the 90 brain regions in the AAL template.Pearson correlation is probably the most commonly used scheme for calculating functional connections (Valdes-Sosa et al., 2011).We also use the Pearson correlation coefficient between vertices to construct the FC matrix in this work.FC is defined as the temporal correlation of brain region.For convenience, we defined P(j) = {P 1 (j), P 2 (j), • • •, P n (j)} (j = 1, 2, • • •, 90) as the average time signal sequence of the j-th brain regions, and n is the total number of time series.The FC matrix M = (m ij ) (i, j = 1, 2, • • •, 90) is calculated using equation (1): m corr P P P P P P To construct VR complex, we defined the distance between each two vertices using equation ( 2):  is constructed.Thus, a VR complex was further established for each subject.The VR complex take three inputs: the maximum dimension d max of any included simplex, the maximum filtration value t max , and the number of divisions.In our experiments, the maximum dimension d max = 2, which means we construct 0-dimensional simplex and 1-dimensional simplex.The maximum filtration value t max = 0.9 and the number of division is set to 450.

Discussion
The VR complex is established under multiple scales.The 0-dimensional Betti number β 0 and the 1-dimensional Betti number β 1 under each scale are calculated using Javaplex software.Figure 3 plots show the 0-dimensional and 1-dimensional Betti number curves between PD and NC.As depicted from the 0-dimensional Betti number diagram (Figure 3A), with the increase of threshold δ, the number of connected branches of PD and NC gradually decrease from the initial 90.Until the threshold δ = 0.7, they tend to be consistent and gradually decreased to 1.But β 0 in PD is always greater than that in NC when δ ≤ 0.7.In the 1-dimensional Betti number curves (Figure 3B), the topological features was significantly different between the PD brain network and the NC brain network under different thresholds.
In Table 2, two sample Kolmogorov Smirnov (K-S) tests were also conducted from Figures 3A,B to compare whether there is a significant difference in distribution between the PD group and the NC group.
From Table 2, it can be seen that the maximum absolute differences in the cumulative probability of Betti numbers in 0-dimension and 1-dimension are 0.049 and 0.082, respectively.Assuming a significance level of 0.05, as the probability p-values (both 0.000) are less than the significance level, it can be concluded that there is a significant difference in the Betti number curve between the PD group and the NC group, regardless of whether it is 0-dimensional or 1-dimensional.Therefore, by comparing the Betti numbers at all scales, the differences between these two groups can be detected.
To more intuitively show the difference between PD and NC groups, comparison of the 1-dimention FCNC in the PD group and the NC group with δ = 0.4 in the Neurocon dataset are shown in Figure 4.Among them, the coronal and sagittal views are shown in Figure 4A and Figure 4B respectively, and their corresponding relationship matrices are given in Figure 4C.Note that, the different colors in Figure 4 only represent different FCNC.
The number of FCNC in 90 brain regions was analyzed and a histogram was drawn in Figure 5.It can be seen intuitively, the average   number of FCNC in PD and NC groups is significantly different in some brain regions.To further analyze these differences, the mean and standard deviation of the number of FCNC in each brain region were calculated for the PD and NC groups under thresholds of 0.7, 0.8, and 0.9 in Neurocon dataset, respectively.Then, two sample Mann Whitney U test and false discovery rate (FDR) correction were performed to detect differences between the two groups.The brain regions with statistical differences are shown in Table 3.In PD group, under threshold of 0.7, the number of FCNC involved was significantly differences and these brain regions include the Cuneus_R, Lingual_R, Fusiform_R and Heschl_R.In addition to the aforementioned brain regions, there are also significant differences in brain regions in the Frontal_Inf_Orb_R and Pallidum_R, when the threshold increases to 0.8 and 0.9 (p < 0.05 and FDR correction).Comparison of the average number of FCNC between PD and NC groups in AAL 90 brain regions.

0-dimension 1-dimension
10.3389/fnins.2023.1236128 Frontiers in Neuroscience frontiersin.org We also analyzed the statistical differences in the number of FCNC with medium lengths in the Neurocon and PPMI datasets (Table 4).From Table 4, it can be seen that there is a significant difference between the PD group and the NC group for medium length of FCNC.Specifically, at thresholds of 0.8 and 0.9, there was a significant difference between the two groups in the Neurocon dataset (length = 8).When the thresholds are 0.7, 0.8, and 0.9, there was a significant difference between the two groups in the PPMI dataset (length = 9).All results were corrected through FDR correction.In these two datasets, at different thresholds, the mean of the PD group is always greater than that of the NC group.
At different thresholds, we can also visually see the significant differences between the two groups through the box-plot in the Neurocon dataset.From the box-plot Figure 6, it can be seen that  when the threshold is there is no significant difference in the median between the two groups.However, when the thresholds are increased to 0.8 and 0.9, the median of the PD group are significantly higher than that of the NC group.The above results show that there are significant differences in the characteristics of FCNC in some brain regions, and the medium length of FCNC in PD patients show significant changes.
By analyzing the FCNC, we found significant statistical difference in the Frontal_Inf_Orb, Cuneus, Lingual gyrus, Fusiform, Pallidum and Heschl areas.A similar result was obtained in previous studies regarding the function connectivity strength of the SMN and VAN (Caspers et al., 2021;Tsuboi et al., 2021) and our findings provide guidance for further studies on the pathogenesis of early PD.Interestingly, increases in FC within SMN (Pang et al., 2021) have been observed upon a dopaminergic challenge in PD patients.Both anterior central gyrus and transverse Nie gyrus belong to SMN.Caspers et al. (2021) pointed out that PD is accompanied by the loss of functional connection of SMN, whether within the network or in the interaction with other networks.The lesions of globus pallidus can have symptoms such as increased muscle tension, decreased movement and static tremor (such as Parkinson's syndrome).Globus pallidus plays an important role in the regulation of motor function.It is not only the relay nucleus between the caudate putamen (CPU) and subthalamic nucleus (STN), but also integrates the inhibitory afferent from CPU and the excitatory afferent from STN, neocortex and thalamus, thus affecting the efferent signal of basal ganglia (Esposito et al., 2013).Globus pallidus stimulation can be used to improve the brain connectivity in order to treat advanced PD (Tsuboi et al., 2021).

Conclusion
This study applies persistent homology to the brain functional networks of PD.This work provides some new insights into the evolution of functional network in the progression of PD and may provide evidence for the study of preclinical biomarkers of PD.We observed that there are significant differences in the characteristics of FCNC in some brain regions, and the medium length of FCNC in PD patients show significant changes.Topological analysis based on rsfMRI data may provide comprehensive information about the changes of FCNC and may provide an alternative for clinical differential diagnosis.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers.Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Figure 1
Figure1gives an example of Betti number changes with δ = 1.4,5.2, 6.3, and 8.5, respectively.Take 10 points randomly and draw a

FIGURE 3
FIGURE 3 Betti number curves of PD and NC groups in 0 and 1 dimensions.(A) The 0-dimensional Betti number diagram.(B) The 1-dimensional Betti number diagram.

FIGURE 4
FIGURE 4 Comparison of the FCNC in the PD group and the NC group with δ = 0.4 in the Neurocon dataset.(A) Coronal view.(B) Sagittal view.(C) The relation matrix.The first row: the FCNC in the PD group.The second row: the FCNC in the NC group.Different colors are just for the purpose of easy observation.

Table 1 .
Retrieve rsfMRI data from the PPMI dataset for 154 PD patients and 24 age matched NC patients.Each study in the PPMI dataset was approved by the Human Experimental Ethics Standards Committee before registration, and each subject signed a written informed consent form.This study obtained the right to use the PPMI database data.

TABLE 1
The details of the Neurocon dataset.

TABLE 2
Two-sample K-S tests.

TABLE 3
Statistically significant differences in the numbers of FCNC in the involved brain regions.

TABLE 4
Statistically significant difference between the two groups in the number of medium length of FCNCs.