Neural Plasticity

Volume 2016, Article ID 1938292, 8 pages

http://dx.doi.org/10.1155/2016/1938292

## Node Detection Using High-Dimensional Fuzzy Parcellation Applied to the Insular Cortex

^{1}GCS fMRI, Koelliker Hospital, Turin, Italy^{2}Functional Neuroimaging and Neural Complex System Group, Department of Psychology, University of Turin, Turin, Italy^{3}Department of Psychology, University of Turin, Turin, Italy^{4}Neuroscience Institute of the Cavalieri Ottolenghi Foundation and Department of Neuroscience, University of Turin, Turin, Italy

Received 29 April 2015; Revised 29 June 2015; Accepted 6 July 2015

Academic Editor: Yong Liu

Copyright © 2016 Ugo Vercelli et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Several functional connectivity approaches require the definition of a set of regions of interest (ROIs) that act as network nodes. Different methods have been developed to define these nodes and to derive their functional and effective connections, most of which are rather complex. Here we aim to propose a relatively simple “one-step” border detection and ROI estimation procedure employing the fuzzy -mean clustering algorithm. To test this procedure and to explore insular connectivity beyond the two/three-region model currently proposed in the literature, we parcellated the insular cortex of 20 healthy right-handed volunteers scanned in a resting state. By employing a high-dimensional functional connectivity-based clustering process, we confirmed the two patterns of connectivity previously described. This method revealed a complex pattern of functional connectivity where the two previously detected insular clusters are subdivided into several other networks, some of which are not commonly associated with the insular cortex, such as the default mode network and parts of the dorsal attentional network. Furthermore, the detection of nodes was reliable, as demonstrated by the confirmative analysis performed on a replication group of subjects.

#### 1. Introduction

One powerful method for studying brain organization is the graph-theoretical approach [1]. Similar to other connectivity methods, such as seed-based functional connectivity and diffusion-tensor-imaging approaches, this method requires the definition of a set of regions of interest (ROIs) that act as network nodes [2]. Several techniques have been employed to functionally derive these nodes and their connections [3, 4]. The definition of such nodes often involves complicated functional connectivity estimation and border detection procedures [5]. Here, we suggest a relatively simple “one-step” border detection and ROI estimation procedure. In particular, we propose to take advantage of one of the characteristics of the fuzzy -mean clustering algorithm [6]. This procedure allows a fixed percentage of voxels with a borderline pattern of connectivity to be nonunequivocally attributed. The further we move from the centre towards the border of a cluster, the more the characteristics of the pattern of connectivity are intermixed with those of the neighbouring clusters (e.g., nonunequivocally determined). As was recently shown by Smith et al. [7], the maximization of spatial independence could lead to suboptimal detection of networks that share significant spatial overlaps. Our method maximizes the temporal independence because the criterion used in the clustering algorithm is the correlation between the time series of each voxel and the time series of the centre of each cluster.

Previous clustering studies [6, 8–10] performed the clustering procedure at subject level. Given the relative deficiency of time points (about 120–200 points in a six-minute run), this procedure has good reliability only for low-dimensional parcellation (e.g., with a limited number of clusters). In line with other investigations [2, 11, 12], we concatenated the time courses across all subjects to constitute a very big dataset that allowed us to obtain a higher clustering dimensionality. Far from representing a step backwards, this “fixed-effect” approach permits good estimation of a common set of clusters (and thereby of nodes) for a given group of subjects. Subsequently, the between-subject variance was taken into consideration, so that each node’s functional connectivity pattern was evaluated at the subject level and summarized using a random-effect analysis. This method is very similar to the dual regression approach [13], according to which the independent component analysis (ICA) was first applied to the concatenated dataset.

To investigate how this node detection method performs with real data, we applied our procedure to the insular surface of 20 healthy subjects scanned in a resting state. A second dataset of 18 healthy volunteers was used for replication testing. We chose the insular surface because the insula is a complex and pivotal [14] brain area in which different inputs from the body and the external world are highly integrated [15]. This brain region has been parcellated by using different measures, such as resting-state functional connectivity [8, 10, 14], task-related functional connectivity [16, 17], and diffusion tensor imaging [18, 19], into two [8, 16, 20–23], or three [9, 10], or more clusters [4, 5, 17, 24], each of which has a unique pattern of connectivity. A recent paper by Kelly et al. [25] demonstrated a convergence between resting state, task-based functional connectivity, and anatomical coactivations at several different parcellation levels (from two to 12 clusters per side; with more than 12 clusters, reliability dropped by about 50%), thus supporting a common hierarchical structure within the insular cortex.

Given these premises, it would be of great interest to test this new node detection method in search of a resting-state functional parcellation of the insular surface. As an additional consideration, we suppose that high-dimensional clustering [11] can make it possible to demonstrate the existence of a more complex pattern with “echoes” [12] of several brain networks nested within the two main insular patterns previously reported.

#### 2. Methods

##### 2.1. Main Group

Main group consists of twenty healthy right-handed volunteers (10 females, with a mean age of 32.6 ± 11.2). Replication group comprises eighteen healthy right-handed volunteers (nine females, with a mean age of 25.3 ± 4.2). All subjects were free of neurological or psychiatric conditions, were not taking any medication known to alter brain activity, and had no history of drug or alcohol abuse. Handedness was ascertained with the Edinburgh Inventory [26]. We obtained written informed consent from every subject, in accordance with the Declaration of Helsinki. The study was approved by the institutional committee on ethical use of human subjects at the University of Turin.

##### 2.2. Task and Image Acquisition

Images were acquired during a resting-state scan on a 1.5 Tesla INTERA scanner (Philips Medical Systems). Functional weighted images were acquired using echoplanar (EPI) sequences, with a repetition time (TR) of 2000 ms, an echo time (TE) of 50 ms, and a 90° flip angle. The acquisition matrix was 64 × 64, with a 200 mm field of view (FoV). A total of 200 volumes were acquired, with each volume consisting of 19 axial slices; slice thickness was 4.5 mm with a 0.5 mm gap, while in-plane resolution was 3.1 mm. Two scans were added at the beginning of functional scanning to achieve steady-state magnetization before acquiring the experimental data. A set of three-dimensional high-resolution T_{1}-weighted structural images was acquired, using a fast field echo (FFE) sequence, with a 25 ms TR, an ultrashort TE, and a 30° flip angle. The acquisition matrix was 256 × 256 and the FoV was 256 mm. The set consisted of 160 contiguous sagittal images covering the whole brain.

##### 2.3. Data Analysis

Datasets were preprocessed and analysed using BrainVoyager QX software (Brain Innovation, Maastricht, The Netherlands).

Functional images were preprocessed to reduce artefacts as follows [27]: (i) slice scan time correction was performed using a sinc interpolation algorithm; (ii) 3D motion correction was applied using a trilinear interpolation algorithm according to which all volumes were spatially aligned to the first volume by rigid body transformations, and the roto-translation information was saved for subsequent elaborations; (iii) spatial smoothing was performed using a Gaussian kernel of 8 mm FWHM. Several nuisance covariations were regressed out from the time courses to control for the effects of physiological processes, such as fluctuations related to cardiac and respiratory cycles and motion [28–30]. Specifically, we included nine additional covariations from white matter (WM), global signal (GS) [31], and cerebrospinal fluid (CSF), as well as six motion parameters. Subsequently, time courses were temporally filtered in order to keep only frequencies between 0.008 and 0.08 Hz and normalized.

Following the preprocessing, we implemented some steps to improve intersubject analysis of the anatomical location of brain activations. For each subject the functional scans were coregistered with a relatively high-resolution structural scan. This coregistration was done using both the slice positioning as stored in the raw data’s headers and fine adjustments calculated comparing the intensity values of the data sets. After this we transformed each subject’s 3D structural data into Talairach space [32]. This was obtained by translating and rotating the cerebrum on the plain passing through the anterior and the posterior commissure; then, the borders of the cerebrum were identified. The coregistration matrix of anatomical and functional data consisted of the parameters of rotation and translation during the coregistration step and the parameters of Talairach normalization. Finally, by applying the anatomical-functional coregistration matrix we transformed into Talairach space the functional time course of each subject and created the volume time course.

We applied a fuzzy -mean algorithm to the time courses of all the insular voxels and clustered these voxels on the basis of their temporal similarity. As is typical of fuzzy clustering techniques, a certain percentage of voxels can be nonunivocally attributed to the parcels. The percentage of nonunivocally attributed voxels (the fuzziness coefficient) is an arbitrary parameter. In line with other studies [33], we chose 20% of nonunivocally attributed voxels.

The fuzzy clustering technique parcels out a subset of voxels in “clusters” of activation [34]. Signal time courses of all voxels were* z*-standardized. We subsequently confronted the voxel’s time courses () with each other and derived a representative cluster of time courses (cluster centroids) (). On the basis of this unsupervised method, and starting from the original fMRI time series, we got a predefined number of spatial modes, which were composed of a spatial map and an associated centroid time course. Accordingly, a voxel is assigned to a cluster with reference to the similarity (e.g., by correlation) of its time course to the cluster centroid. This similarity is determined in a fuzzy way, which means that a voxel is not uniquely assigned to one cluster (hard clustering) because the similarity is expressed by the “membership” of voxel to cluster .

Centroids and memberships are both updated in an iterative procedure [35], which terminates when successive iterations do not further change memberships. Cluster centres are determined via classical cluster-algorithm distance measures and are expressed as follows:where is the distance between a voxel and a cluster centre and is the coefficient that determines the fuzziness of the procedure; “tunes out” the noise in the data and lies between 1 (smallest fuzziness) and infinity. The most commonly used of the several distance measures of are the Euclidean distance, , and the Mahalanobis distance, [36], which are defined as is the covariance matrix of cluster . The Mahalanobis distance takes in the elliptical shape of the cluster (i.e., it weights the differences by the range of variability, described by , in the direction of the voxel instead of treating all voxels equally when calculating the distance to the cluster centre ). The Euclidean distance assumes a spherical shape, without taking into account the shape of clustering, which corresponds to a covariance matrix with 1 s on the main diagonal and 0 s elsewhere.

Calculation starts from an initial set of membership values for the data set in the following matrix form: with and a matrix of randomly chosen cluster centres with initial . Subsequently, the new centroids and memberships are calculated using (2). When further iterations do not cause significant change to memberships and centroids, the procedure stops. With this procedure the following function is minimized:

This formula calculates the within-class variance over all clusters, . In other words, a user-defined threshold for change in is fixed when convergence is reached. The criterion of convergence is based on the first local optimum of (5) [37]. Each time series is transformed into its -score in order to avoid a classification of the voxels on signal amplitude rather than on signal shape. Then, the principal component analysis (PCA) is performed to reduce data dimensionality. The number of PCA components was calculated to retain 95% of the variance.

Optimal number of clusters: The a priori definition of the number of clusters and the fuzziness coefficient is often debated in the literature [38]. Usually, the optimal number of classes is unknown in fuzzy clustering. A number of cluster-validity indices have been proposed to estimate the optimal number of clusters in an unsupervised manner (for a review see [39]). These indices aim to identify compact and well-separated clusters. In this study, we used the silhouette validation method [40], which consists in considering the silhouette coefficient of each element:where is the average dissimilarity of the -point to all points in the same cluster and is the minimum of the average dissimilarity of the -point to all points in the other cluster.

Unlike Cauda et al. [8], to perform the clustering procedure time courses were concatenated across all subjects; this step, as has been pointed out by others [2], makes it possible to obtain a higher clustering dimensionality (i.e., more clusters).

Due to the fuzziness coefficient employed, 20% of voxels were classified as nonunequivocal. We considered these voxels as border voxels, with a time course that showed transitional characteristics between contiguous clusters.

The final step of this procedure was to place a spherical ROI with a radius of 3 mm in the local maxima of each cluster (i.e., the area of maximal similarity between voxel time courses). See Figure 1 for a graphical representation of the method.