Data-Driven Parcellation Approaches Based on Functional Connectivity of Visual Cortices in Primary Open-Angle Glaucoma

Abstract

Purpose: Functional changes have been observed between diseased and healthy subjects, and functional brain atlases derived from healthy populations may fail to reflect functional characteristic of the diseased brain. Therefore the aim of this study was to generate a visual atlas based on functional connectivity from primary open-angle glaucoma (POAG) patients and to prove the applicability of the visual atlas in functional connectivity and network analysis.

Methods: Functional magnetic resonance images were acquired from 36 POAG patients and 20 healthy controls. Two data-driven approaches, K-means and Ward clustering algorithms, were adopted for visual cortices parcellation. Dice coefficient and adjusted Rand index were used to assess reproducibility of the two approaches. Homogeneity index, silhouette coefficient, and network properties were adopted to assess functional validity for the data-driven approaches and frequently used brain atlas. Graph theoretical analysis was adopted to investigate altered network patterns in POAG patients based on data-driven visual atlas.

Results: Parcellation results demonstrated asymmetric patterns between left and right hemispheres in POAG patients compared with healthy controls. In terms of evaluating metrics, K-means performed better than Ward clustering in reproducibility. Data-driven parcellations outperformed frequently used brain atlases in terms of functional homogeneity and network properties. Graph theoretical analysis revealed that atlases generated by data-driven approaches were more conducive in detecting network alterations between POAG patients and healthy controls.

Conclusions: Our findings suggested that POAG patients experienced functional alterations in the visual cortices. Results also highlighted the necessity of data-driven atlases for functional connectivity and functional network analysis of POAG brain.

Brain atlases that provide network node for functional connectivity and network analysis are crucial, and many atlases could achieve this process.¹ There are mainly two types of brain atlas. One is based on anatomic landmarks or structural connectivity²; traditional brain atlases, such as the Brodmann atlas and automated anatomic labeling atlas, which are based on cerebral structural characteristics (e.g., cytoarchitecture), belong to this type of atlas.^3–5 The other is the functional brain atlas,² which uses information derived from functional images to parcellate the brain into several subregions.^2–5 Functional magnetic resonance imaging (fMRI) acquired during resting-state, known as resting-state fMRI, could reveal both functional homogeneity and functional connectivity.^2–5 Brain regions based on functional connectivity profiles could provide credible evidence for the establishment of network nodes.² Several state-of-the-art studies have used functional connectivity profiles derived from resting-state fMRI to generate a functional brain atlas.^6–10 Data-driven techniques that have been used to generate functional brain atlases include independent component analysis,^11–13 Gaussian mixture model,¹⁴ and clustering algorithms.¹⁵ 

Visual cortices, which consists of the primary visual cortex (also known as the striate cortex) and higher visual cortex (also known as the extrastriate cortex), are located in the occipital cortex on both sides of the talus fissure of the brain.¹⁶ Parcellation of visual cortices is crucial for neuroscience studies of varies visual diseases.¹⁵ A visual atlas based on brain function rather than cytoarchitecture may provide an efficient way to study disease-related functional changes in the brain.^15,17,18 In addition, the functional atlas is more likely to be generated from resting-state fMRI data.^17,18 As a result, more and more studies on parcellation of visual cortices have used functional connectivity profiles from resting-state fMRI to yield subregions instead of using structural information.^18,19 

Primary open-angle glaucoma (POAG) is a progressive optic neuropathy associated with retinal ganglion cell loss and optic nerve damage.^20–23 Gupta et al.²⁴ have discovered that POAG not only leads to optic nerve atrophy, but also changes the visual cortices and lateral geniculate body. Subsequently, several studies have proved that POAG is a degenerative disease of the central nervous system, analogous to Alzheimer disease and Parkinson disease.^25–27 There are a growing number of studies on POAG patients’ central nervous system via neuroimaging techniques, especially fMRI.^23,27–30 However, previous studies on functional connectivity and network analysis of POAG have mainly used atlases generated from healthy populations, which may fail to reflect functional characteristics of the diseased brain.^28,29 To date, there has been no functional atlas of the visual cortices generated from POAG patients and for the study of the POAG brain. 

Data-driven parcellation approaches based on resting-state fMRI could assess functional variation across different subjects, or between healthy and diseased subjects, which are suitable means for generating functional brain atlases for diseased subjects.² Therefore we aim to generate a functional atlas of the visual cortices for POAG patients, which reflects functional characteristics of the POAG brain. We further investigated altered network patterns of POAG patients using the proposed atlas to prove the applicability of such atlas in functional connectivity and network analysis. Due to relatively high performance reported by previous studies,^15,31 two data-driven approaches, K-means and Ward algorithm, were selected to generate the atlas. 

Materials and Methods

Subjects

The protocol of this study was approved by the institutional review board of Shandong First Medical University and Taian City Central Hospital according to the tenets of the Declaration of Helsinki. All subjects gave their written, informed consent before participating in this study. 

Participants were recruited with the following criteria: (1) right-handed; (2) more than 12 years of formal education; (3) age 40 to 60 years old; (4) and were POAG patients based on a clinical diagnosis of primary glaucoma, and at least one eye with elevated intraocular pressure (IOP) (higher than 21 mm Hg); healthy controls (HCs) were included if they had normal IOP. Exclusion criteria were: (1) advanced-stage glaucoma, secondary glaucoma, or other oculopathies; (2) hypertension, impaired glucose tolerance, diabetes, or other metabolic disease; (3) history of psychiatric or neurologic diseases; and (4) MRI contraindications. Finally, 36 POAG patients and 20 age- and gender-matched HCs were recruited. All POAG patients were in the early- and intermediate-stage. POAG patients underwent a detailed ophthalmology examination, including IOP measurement, retinal nerve fiber layer (RNFL) thickness measurement, and optic disc evaluation. Supplementary Table S1 lists the demographic and clinical information of the POAG group and HCs. 

Connectivity Matrix and Parcellation Approaches

The fMRI acquisition and preprocessing are shown in Supplementary material section 1. After data preprocessing, a connectivity matrix was constructed before data-driven parcellations. The rows of the connectivity matrix denote voxels in the visual cortices, the columns denote voxels in the gray matter. In this study, the visual cortices, which contained primary and higher visual cortices, were determined as follows: primary visual cortex was defined as Brodmann area 17, higher visual cortex was defined as the region that contained Brodmann area 18 and 19. Masks of the visual cortices and whole brain were generated. 

The connectivity matrix characterizes functional connectivity between visual cortices and the brain gray matter. To construct the matrix, we used visual cortices as the first mask. Each time course within the visual cortices was normalized to have zero mean and unit length.^8,9,32,33 This process was implemented by the following formula^17,34:  

\begin{equation}{\frac{{{v_i} - \overline {{v_i}} }}{{\left\| {{v_i} - \overline {{v_i}} } \right\|}}_2},\,i = 1,\,2,\, \ldots \,N\end{equation}

(1)

where v_i denotes the time course of the ith voxel in the visual cortices, \(\overline {{v_i}} \) denotes the mean value of this time course, ‖.‖₂ denotes L2-norm of a vector. Then we used the brain gray matter as the second mask and normalized the time courses of each voxel in the second mask using Equation (1). Assume v_i is the ith voxel in the visual mask, and w_j is the jth voxel in the gray matter mask. The Pearson correlation coefficient between the two time courses v_i and w_j could be simply calculated by multiplication of the two normalized time courses^17,34:  

\begin{equation}{\rm{co}}rr\left( {{v_i},{w_j}} \right) = {v_i} \cdot {w_j}\end{equation}

(2)

 

Assume V is the fMRI data within the visual cortices with normalized time courses, and W is the fMRI data in the gray matter with normalized time courses. The individual connectivity matrix could be calculated by VW^T as shown in Equation (2). 

Then we used individual connectivity matrices to generate a group-averaged connectivity matrix. We first converted the individual connectivity matrices into z maps by Fisher r-to-z transformation to increase the normality of the distribution of the correlations,^17,34 then all the individual z-maps were averaged to yield the averaged z-map. In the end, we converted the averaged z-map back to the averaged connectivity matrix via inverse Fisher r-to-z transformation. 

To improve the spatial contiguity of the parcellation, spatial coordinates (with spatial information in three directions) of each voxel in the visual cortices were combined with the group-averaged connectivity matrix. State-of-the-art studies have reported that K-means and Ward algorithm outperformed alternative approaches in terms of functional homogeneity and reproducibility.^15,31 Therefore the two data-driven approaches, K-means and Ward clustering, were applied to the combined matrix to get the final parcellations. Supplementary Table S2 lists the details of the two approaches, and Figure 1 demonstrates the parcellation procedures of each approach. It was worth mentioning that no threshold was applied to control negative and weak correlations as in previous studies^8,34 because a hard threshold might yield clusters with comparable shapes and sizes, which were unlikely to be functional subregions.³⁴ In addition, negative correlation also contained useful information, especially for disease-related parcellation. 

Figure 1.

 

The procedures of two data-driven parcellation approaches.

Evaluation Metrics

To assess reproducibility and homogeneity of data-driven parcellation approaches, we adopted four metrics including dice coefficient, adjusted Rand index (ARI), homogeneity index, and silhouette coefficient (SC).^32,33 Dice coefficient and ARI were applied for assessing reproducibility, homogeneity index and SC were applied for evaluating functional homogeneity. In addition, we evaluated functional network characteristics constructed by the underlying approaches using nodal and global properties. 

The detailed description of dice coefficient and ARI is given in Supplementary material section 2. To calculate these two metrics, the total 36 POAG subjects were randomly divided into two groups, with one group containing 18 subjects and the other containing the remaining subjects. The group-averaged connectivity matrices of the two groups were calculated as introduced earlier. After that, we parcellated the combined matrices and generated parcellation results (X and Y) of the two groups. The parcellation results X and Y from the two groups were used to compute dice coefficient and ARI. This random division was repeated 10 times, and dice coefficient and ARI were calculated 10 times. The final dice coefficient and ARI were then calculated as the mean value of these 10 random repeats. 

Function homogeneity of subregions was crucial for functional network analysis, so homogeneity index and SC were used to evaluate functional homogeneity. The detailed description of homogeneity index and SC is given in Supplementary material section 3. To compare functional homogeneity between data-driven approaches and frequently used brain atlases, we selected four brain atlases, namely, the Brodmann atlas, and three functional atlases from Craddock et al.,⁸ Dosenbach et al.,¹⁰ and Power et al.⁹ Homogeneity index and SC were calculated for the underlying approaches in the POAG group. 

To compare network properties between data-driven approaches and the four frequently used brain atlases, we followed the procedure by Arslan et al.¹⁵ Once group-level parcellation for POAG patients was done, functional connectivity was then estimated among the time courses by calculating Pearson correlation coefficient. An N × N correlation matrix for each subject was generated, where N is the parcellation number. Fisher r-to-z transformation was performed to improve the normality of the correlation matrices. A weighted network for POAG patients was obtained and sparsity threshold was used to preserve the top 20% of the edges and to reduce threshold effects on network properties.¹⁵ Then a graph theoretical analysis was performed to investigate topological properties of the network using GRETNA (version 2.0, http://www.nitrc.org/projects/gretna/).³⁵ Due to granularity of parcellation,¹⁵ we only applied nodal degree, efficiency to evaluate nodal characteristics, and small-world index and rich club index to evaluate global characteristics when the parcellation number was more than 3. The detailed description of network properties is shown in Supplementary material section 4. 

Graph Theoretical Analysis Between POAG Patients and HCs

Overall, the main purpose of data-driven parcellation is to define nodes for functional connectivity and network analyses.¹² Therefore we evaluated the underlying parcellation approaches and frequently used brain atlas using a standard graph theoretical analysis of brain functional networks between POAG and HCs. The visual atlas generated from POAG patients via K-means algorithm and the four frequently used brain atlases were used. Network construction was the same as introduced earlier. Four topological properties, namely nodal degree, efficiency, small-world index, and rich club index, were used to assess brain network alterations in POAG patients at both global and nodal levels. Two sample t-test was used to assess the difference of network properties between the two groups. 

Results

Visual Cortices Parcellation

To observe functional reorganization between POAG patients and HCs, we also parcellated the visual cortices for HCs. Figure 2 displays parcellation results of POAG patients and HCs when parcellation number was 3, 6, and 9. 

Figure 2.

 

Parcellation results of the two data-driven approaches when the parcellation number is set to be 3, 6, and 9. (a) Parcellation results for 36 POAG patients. (b) Parcellation results for 20 HCs.

Figure 2 illustrates that both algorithms could achieve parcellation of visual cortices. It was found that the functional subregions of POAG patients’ visual cortices were more fragmented than HCs, and POAG patients’ parcellation patterns were asymmetric compared with those of HCs. The visual cortices could be parcellated into the anterior and posterior subregions, which was different from the traditional boundaries between primary visual cortex and higher visual cortex. 

Evaluation Metrics

Reproducibility and functional homogeneity results are shown in Figure 3. In terms of reproducibility, we found that K-means algorithm yielded better results than Ward algorithm in both dice coefficient and ARI. Moreover, reproducibility results evaluated by dice coefficient and ARI of these two approaches decreased with the increase of parcellation numbers. For homogeneity index, K-means and Ward clustering obtained comparable results. When it comes to SC, K-means algorithm had a better performance. Both SC and homogeneity index increased with the increase of parcellation number. 

Figure 3.

 

Group-level reproducibility and homogeneity results in POAG group (mean ± SD). (a) Dice coefficient for 10-time repeats. (b) ARI for 10-time repeats. (c) Homogeneity index. (d) SC.

For comparisons between data-driven approaches and frequently used brain atlases, functional atlases of Dosenbach et al.¹⁰ and Power et al.⁹ had the best performances in homogeneity index, but had the worst performances in SC. In addition, both SC and homogeneity value of the Brodmann atlas and the Craddock et al.⁸ functional atlas were smaller than those of the two data-driven parcellation approaches. Homogeneity results for HCs using parcellations from the POAG group is demonstrated in Supplementary Figure S1. 

The results of network properties are given in Figure 4. In terms of nodal degree and nodal efficiency, the two data-driven parcellation approaches prevailed over frequently used brain atlases. As for global metrics, K-means and Ward approaches performed better than the three atlases in terms of small-world index and rich club index. 

Figure 4.

 

Network properties of the two data-driven parcellation approaches and frequently used brain atlas. (a) Nodal degree. (b) Nodal efficiency. (c) Small-world index. (d) Rich club index.

Results of Graph Theoretical Analysis

Figure 5 shows the results of graph theoretical analysis using data-driven atlas generated from POAG patients via K-means algorithm and the four frequently used atlases. The comparisons between data-driven atlas and the four frequently used atlases were conducted under the same resolution. 

Figure 5.

 

Graph theoretical analysis between POAG patients and HCs using data-driven atlas and frequently used brain atlases. We were not able to calculate the small-world index and rich club index for the two groups using Brodmann atlas due to granularity of parcellation.

POAG patients showed decreased nodal degree and rich club index, and increased nodal efficiency and small-world characteristic compared with HCs. Furthermore, in most cases, the differences of network properties between POAG patients and HCs were more significant using data-driven atlas than using frequently used atlases. A detailed comparison between POAG patients and HCs using data-driven atlas across all resolutions is shown in Supplementary Figure S2. 

Discussion

Brain parcellation of healthy subjects has been carried out by many research groups worldwide,³⁶ and several patterns have been already learned.¹⁵ However, knowledge on parcellation of different brains, such as diseased brain, has been limited. Because parcellation of diseased brain can reveal functional characteristics of the diseased brain,¹⁵ in this study, we parcellated visual cortices of the POAG brain to generate a visual atlas for POAG subjects for the first time to our knowledge. Parcellation results revealed functional reorganization in the visual cortices of POAG patients. Graph theoretical analysis demonstrated altered network patterns of POAG patients, and also suggested that data-driven approach was necessary for studies of functional connectivity and functional network of POAG patients. 

Functional Alterations in the Visual Cortices

In the present study, parcellation results revealed that HCs had a symmetric parcellation pattern in the bilateral hemispheres with large and equal clusters, which indicated normal brain function. POAG patients’ parcellation patterns were asymmetrical, some subregions were small and fragmented clusters. Two mechanisms could explain this, one is lateralization of brain functional changes, in other words, changes in the central visual cortices does not simply mirror the peripheral damage in a one-to-one manner.^30,37 The other is the difference between the two eyes in terms of IOP, RNFL thickness, in other words, functional changes of the visual cortices may not be controlled by exact same factors. 

Visually assessed results showed that V1 could be roughly identified by data-driven approaches, especially when parcellation number was greater than or equal to 6. The boundaries of subregions in the higher visual cortex were not consistent with cytoarchitecture. The results also showed that visual cortices, including V1 and higher visual cortex, could be functionally divided into anterior and posterior subregions, which was consistent with previous studies.^7,18,38 State-of-the-art studies have demonstrated that the visual cortices could be divided into peripheral and foveal areas belonging to different resting-state networks with eccentricity-dependent functional connectivity.^11,18 In line with previous findings, our results might reflect distinctively biased processes in the peripheral and foveal networks.¹¹ From the sagittal view, we found that the visual cortices of HCs could be roughly divided into two parts: the superior part near the parietal lobe and the inferior part near the temporal lobe, which might reflect the dorsal and ventral pathway in visual information transmission.¹⁶ However, the parcellation patterns of POAG patients could not be clearly divided into superior and inferior parts, which might be due to some defects in the process of visual information transmission. 

Graph theoretical analysis demonstrated differences in network properties between the two groups. Studies have revealed that in the POAG brain, several visual areas had decreased spontaneous brain activities and decreased connectivities with the rest of the brain compared with HCs.^25–27 This might account for the decreased nodal degree and rich club index in POAG patients compared with HCs. Previous studies have reported functional reorganization in the visual cortex in the POAG brain,^23,28,29 in other words, the original functional connections might be blocked and some long-range connections might be generated in the visual networks, leading to small-world topology for the visual network and enhancement of network efficiency. In addition, several studies have reported functional compensation in the POAG brain with increased functional connectivities between normally functioning regions.^28,29 The increased small-world index and nodal efficiency may also be attributed to functional compensation. 

Evaluating Metrics

In terms of evaluating metrics, K-means clustering had better performance than Ward clustering in terms of reproducibility, which was consistent with the previous study.¹⁵ The functional atlas of Dosenbach et al.¹⁰ and Power et al.⁹ achieved better scores in homogeneity index, which was because these two atlases defined each subregion as a sphere region-of-interest (ROI) with a radius of 5 mm; smaller subregions tended to have a higher homogeneity index. In other words, these two atlases were discontinuous atlases. This could also be seen from the results of SC. Because SC not only assesses within-parcel similarity but also assesses interparcel dissimilarity, making it more sensitive to spatial contiguity.¹⁵ The worst performance of the Dosenbach et al.¹⁰ and Power et al.⁹ atlases were attributed to their discontinuity. Besides, it was easy to find that data-driven approach performed better than Brodmann atlas and the Craddock et al.⁸ atlas in SC and homogeneity index, which implied that data-driven approach could better characterize functional organization of the POAG brain than traditional brain atlases that were generated from healthy populations.³ 

Among different network properties, nodal degree is the number of edges attached to the node, which reflects nodal information communication ability in the network.³⁵ Nodal efficiency characterizes the efficiency of parallel information transfer in the network.³⁵ Data-driven parcellation approaches had higher nodal degree and efficiency, which reflected higher information communication and transmission in the network. Based on the earlier described results, it could be concluded that functional features captured by the data-driven atlas may be more representative of the POAG brain.¹⁵ 

Applicability of Data-Driven Atlas in fMRI Studies of POAG

In graph theoretical analysis, all atlases revealed the same type of differences (POAG < HC or POAG > HC) in network properties between the two groups. Also, the application of data-driven visual atlas facilitated the differences in graph theoretical analysis, which is important for future studies. However, under certain circumstances the brain ROIs by Dosenbach et al.¹⁰ and Power et al.⁹ showed more significant results between the two groups. As discussed earlier, these two types of brain ROIs possessed high functional homogeneity inside each ROI, which might help to facilitate the differences in graph theoretical analysis. 

Data-driven approaches parcellated voxels with similar functional connectivity profiles into the same cluster.^15,31 Because POAG patients showed functional alterations in the visual cortices,^25–27 data-driven algorithm would detect the heterogeneity and assigned relevant voxels into a different cluster. As for HCs, brain functions during resting-state were more homogeneous. Therefore data-driven atlas could characterize functional reorganization of the POAG brain, and by using data-driven atlas generated from diseased patients, it was easier to locate differences in functional connectivity and network properties between diseased and healthy subjects. In addition, the intrasubject variability for POAG patients and HCs using the data-driven atlas were comparable to those utilizing frequently used brain atlas, indicating its generalization ability. Therefore it is suggested that the data-driven atlas has potential use in the population-based studies of functional connectivity and functional network of POAG patients. 

There are several limitations needed to be addressed. First, we only parcellated the visual cortices of POAG patients because visual cortices were the main ROI for visual diseases. Second, small sample size, age variability, and disease severity may affect the generalization of the data-driven approaches. Future studies should focus on a larger sample size. Third, the current study was not able to provide a recommended parcellation number for the visual cortices. 

Conclusions

In this study, we parcellated the visual cortices of POAG patients via data-driven parcellation algorithms and found visual cortices could be functionally divided into subregions different from anatomic boundaries. We also demonstrated functional reorganizations and altered network patterns in the visual cortices in the POAG brain compared with HCs, which may provide insight into the pathology of the disease. Data-driven parcellation prevailed over frequently used brain atlas in functional homogeneity and network properties. In addition, our findings suggest that data-driven parcellation approaches with varying resolution, higher network properties, and higher homogeneity may be more suitable for studies of functional connectivity and functional network of the POAG brain than brain atlases generated from healthy populations. 

Acknowledgments

Supported by the Key Research and Development Program of Shandong Province (2017GGX201010; WL), the Natural Science Foundation of Shandong Province (ZR2016HM73; JQ), and the Academic Promotion Programme of Shandong First Medical University (2019QL009). Jianfeng Qiu was supported by Taishan Scholars Program of Shandong Province (TS201712065). 

Disclosure: H. Qu, None; Y. Wang, None; T. Yan, None; J. Zhou, None; W. Lu, None; J. Qiu, None