# SPANOL (SPectral ANalysis of Lobes): A Spectral Clustering Framework for Individual and Group Parcellation of Cortical Surfaces in Lobes

Julien Lefèvre, Antonietta Pepe, Jennifer Muscato, François de Guio, Nadine

Girard, Guillaume M Auzias, David Germanaud

► To cite this version:

Julien Lefèvre, Antonietta Pepe, Jennifer Muscato, François de Guio, Nadine Girard, et al.. SPANOL (SPectral ANalysis of Lobes): A Spectral Clustering Framework for Individual and Group Parcellation of Cortical Surfaces in Lobes. Frontiers in Neuroscience, 2018, 12, 10.3389/fnins.2018.00354. hal-01868581

HAL Id: hal-01868581

https://amu.hal.science/hal-01868581

Submitted on 24 Jul 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# SPANOL (SPpectral ANalysis of Lobes): A Spectral Clustering Framework for Individual and Group Parcellation of Cortical Surfaces in Lobes

Julien Lefèvre<sup>1,2*</sup>, Antonietta Pepe<sup>1,2</sup>, Jennifer Muscato<sup>1</sup>, Francois De Guio<sup>3</sup>,
Nadine Girard<sup>4,5</sup>, Guillaume Auzias<sup>2</sup> and David Germanaud<sup>6,7,8</sup>

**OPEN ACCESS**

**Edited by:**

Pedro Antonio Valdes-Sosa,
Clinical Hospital of Chengdu Brain
Science Institute, China

**Reviewed by:**

Herve Lombaert,
École de Technologie Supérieure,
Canada

Jingxin Nie,
South China Normal University, China

**\*Correspondence:**

Julien Lefèvre
julien.lefevre@univ-amu.fr

**Specialty section:**

This article was submitted to
Brain Imaging Methods,
a section of the journal
Frontiers in Neuroscience

**Received:** 15 March 2018

**Accepted:** 07 May 2018

**Published:** 31 May 2018

**Citation:**

Lefèvre J, Pepe A, Muscato J, De Guio F, Girard N, Auzias G and Germanaud D (2018) SPANOL (SPpectral ANalysis of Lobes): A Spectral Clustering Framework for Individual and Group Parcellation of Cortical Surfaces in Lobes. *Front. Neurosci.* 12:354.
doi: 10.3389/fnins.2018.00354

Understanding the link between structure, function and development in the brain is a key topic in neuroimaging that benefits from the tremendous progress of multi-modal MRI and its computational analysis. It implies, *inter alia*, to be able to parcellate the brain volume or cortical surface into biologically relevant regions. These parcellations may be inferred from existing atlases (e.g., Desikan) or sets of rules, as would do a neuroanatomist for lobes, but also directly driven from the data (e.g., functional or structural connectivity) with minimum a priori. In the present work, we aimed at using the intrinsic geometric information contained in the eigenfunctions of Laplace-Beltrami Operator to obtain parcellations of the cortical surface based only on its description by triangular meshes. We proposed a framework adapted from spectral clustering, which is general in scope and suitable for the co-parcellation of a group of subjects. We applied it to a dataset of 62 adults, optimized it and revealed a striking agreement between parcels produced by this unsupervised clustering and Freesurfer lobes (Desikan atlas), which cannot be explained by chance. Constituting the first reported attempt of spectral-based fully unsupervised segmentation of neuroanatomical regions such as lobes, spectral analysis of lobes (Spanol) could conveniently be fitted into a multimodal pipeline to ease, optimize or speed-up lobar or sub-lobar segmentation. In addition, we showed promising results of Spanol on smoother brains and notably on a dataset of 15 fetuses, with an interest for both the understanding of cortical ontogeny and the applicative field of perinatal computational neuroanatomy.

**Keywords:** cortical parcellation, spectral clustering, Laplace-Beltrami operator, fetal MRI, group-wise analysis

<footnote n="1"><sup>1</sup> Centre National de la Recherche Scientifique, ENSAM, LSIS, Aix Marseille University, University of Toulon, Marseille, France,</footnote>

<footnote n="2"><sup>2</sup> Centre National de la Recherche Scientifique, Institut de Neurosciences de la Timone, Aix Marseille University, Marseille, France, <sup>3</sup> Institut National de la Santé et de la Recherche Médicale, University Paris Diderot, Sorbonne Paris Cité, UMR-S 1161, Paris, France, <sup>4</sup> Centre National de la Recherche Scientifique, CRMBM, Aix Marseille University, Marseille, France,</footnote>

<footnote n="5"><sup>5</sup> APHM, Hopital de la Timone, Service de Neuroradiologie, Marseille, France, <sup>6</sup> CEA, Neurospin, UNIACT, Equipe Neuropédiatrie, Gif-sur-Yvette, France, <sup>7</sup> Institut National de la Santé et de la Recherche Médicale, Université Sorbonne Paris Cité, CEA, UMR 1129, Paris, France, <sup>8</sup> APHP, Hopital Robert-Debré, DHU Protect, Service de Neurologie Pédiatrique et des Maladies Métaboliques, Université Paris Diderot, Paris, France</footnote>
# 1. INTRODUCTION

The existence of brain regions associated to specific functions (Friston, 2002) as a dominant paradigm in neuroanatomy is a legacy from the XIXth century. This localizationist view has proved more complex and several concurrent segmentations of the brain have been proposed based on microarchitecture first (Zilles and Amunts, 2010), then on neuroimaging studies, in particular with MRI in the past 20 years. Brain parcellations have received considerable interest not only because of their deemed underlying biological reality but also because of their advantages in fMRI studies, their ability to reduce the high dimensionality of the data, and their use in object-based strategies to overcome the shortcomings of spatial normalization (Thirion et al., 2006). Data-driven MRI based parcellations are most frequently sought, for instance with fMRI data (Thirion et al., 2014) or anatomical connectivity (Lefranc et al., 2016; Parisot et al., 2016) or both (Glasser et al., 2016). Yet, brain parcellations based only on the anatomy have been mostly built from information provided by experts and manual labeling. For instance the information of manually delineated cortical sulci has been used to provide volumic (Lohmann and von Cramon, 2000; Tzourio-Mazoyer et al., 2002) or surfacic regions (Cachia et al., 2003; Klein and Tourville, 2012), with different strategies to transfer the labeled information by using registration toward an atlas (Tzourio-Mazoyer et al., 2002) or by working more at the level of the single subject (Cachia et al., 2003; Klein and Tourville, 2012). Even model-based parcellations, like the one proposed in Auzias et al. (2016), require a preliminary labeling of some cortical sulci at the subject level. There are, however, recent attempts of unsupervised parcellations procedures, using spectral analysis (Germanaud et al., 2012) or 'sulcal pits' concepts (Auzias et al., 2015).

A full unsupervised segmentation of the cortical surface into anatomically meaningful entities remains a very challenging task and there might even be no a priori reasons to hope to obtain one, at least from a non-comprehensive input dataset. Nevertheless preliminary results of two previous studies of our group (Lefèvre et al., 2014; Pepe et al., 2015) suggested an intriguing relationships between parcellations obtained by spectral clustering and what has been known for a long time as brain lobes (Gratiolet, 1854). Brain lobes are a coarse but well-accepted anatomo-functional segmentation of the brain into 5 parts: frontal, parietal, temporal, occipital and insular. This segmentation is not absolute since for instance a 6th limbic part is sometimes proposed. Besides, if unambiguous landmarks ground some lobar limits (e.g., central sulcus as frontal-parietal lateral limit or parieto-occipital sulcus as a parietal-occipital internal limit), others are rather ill-defined (e.g., in the occipital-temporal continuum).

Spectral clustering has been largely studied, both theoretically and empirically, mostly in the machine learning community (Ng et al., 2002; Von Luxburg, 2007) but also in computer graphics for mesh segmentation purposes (Liu and Zhang, 2004). Spectral clustering requires a laplacian matrix (in the discrete or graph settings) which corresponds to the Laplace-Beltrami Operator (L.B.O.), often considered as the “swiss army

knife” of the geometry processing. Recently this mathematical operator has been more and more popular in the field of neuroimaging, with applications to the description of gyration pattern (Germanaud et al., 2012; Rabiei et al., 2015), to anatomo-functional variability (Lombaert et al., 2015), in shape registration (Lombaert et al., 2013b; Lefèvre and Auzias, 2015), and in shape classification (Lai et al., 2009; Wachinger et al., 2015).

In this article our contributions are threefold:

* We describe a *general framework* to parcellate a group of cortical surfaces in an arbitrary number of connected regions based on a spectral clustering algorithm adapted to triangular meshes. Only the geometry of the brain is considered and a limited supervised information (cingulate pole labeling) can be injected.

* When restricted to 6 regions, the resulting parcellations offer striking *analogies with brain lobes* that can be precisely quantified. The value 6 is optimal regarding a statistical procedure to test the independence of two parcellations.

* We test the robustness of spectral clustering to surface perturbation and in particular to unfolding. Our approach is also tested on fetal cortical meshes, providing unsupervised segmentations of the immature brain that are consistent with the adult stage. Those results suggest that folding pattern itself play a minor role in the definition of spectral lobes, that could be mainly determined by the global shape of the brain.

# 2. METHODS

## 2.1. Spectral Clustering

Spectral clustering has become a very popular method in machine learning to segment a set of points in different groups based only on pairwise similarities between those points. A dimension reduction step is performed on a matrix derived from the similarities by extracting its *K* first eigenvectors. Then, a classic clustering algorithm such as a *K*-means is performed on these eigenfunctions (Ng et al., 2002; Von Luxburg, 2007).

Graph laplacians are often chosen to define the similarity matrix because of suitable mathematical properties (non-negativity and semi-definiteness) that make them diagonalizable. In our application to cortical surfaces, it is crucial to use a matrix that makes sense in terms of shape and allows fast computations. The Laplace-Beltrami Operator and its discrete representation with Finite Element Methods (F.E.M.) are convenient because they satisfy the required mathematical properties, no parameters tuning is necessary, and the derived matrices are sparse which guarantees efficient computation of eigenvectors thanks to Krylov subspace methods as implemented in ARPACK library. Besides it is interesting to note that, following Belkin and Niyogi (2008), the Laplace-Beltrami Operator can be approximated by a point cloud Laplace Operator which corresponds exactly to an affinity matrix with a Gaussian kernel. Therefore there is an asymptotic equivalence between spectral clustering with Laplace-Beltrami Operator and the traditional one with affinity matrix.

The resulting eigenvectors correspond to an approximation of the *Fourier modes* of the ideal cortical surface. They reflect the intrinsic geometry of the brain shape, independently of how it is
embedded in the 3D space. From a more physical point of view, Fourier modes correspond to vibration modes of the surface, following a close analogy to the well-known Chladni plates (Chladni, 1787). Mathematically, Fourier modes correspond to real valued functions $\phi_i$ defined on a genus-0 surface $S$ that satisfy:

$$-\Delta_S \phi_i = \lambda_i \phi_i. \qquad (1)$$

where $\lambda_0 = 0 \le \lambda_1 \le \lambda_i...$ are eigenvalues that can be interpreted as spatial frequencies. $\Delta_S$ is the Laplace Beltrami Operator of the surface $S$ (see Berger, 2003 for more details).

Visually Fourier modes reveal spatial oscillations whose number increases with $i$. More precisely, the number of connected regions with constant sign (also called nodal domains) is bounded between 2 and $i + 1$, because of the celebrated Courant Nodal theorem. Such eigenfunctions/Fourier modes are represented in the Figure 8.

On a triangular mesh approximating $S$ with $N$ vertices, it is possible to apply the framework of Finite Element Method. A discretized Fourier mode and its corresponding frequency can then be obtained as a vector $U = (u_i)_{i=1..N}$ and a scalar $\lambda$ such that:

$$G U = \lambda M U, \qquad (2)$$

where the expression of the two matrices, in the case of first order F.E.M., are given in Lefèvre and Mangin (2010).

The two steps of the spectral clustering adapted to the Laplace-Beltrami Operator are summarized on the left block of Figure 1 and in Algorithm 1.

<caption>**Algorithm 1** Cortical surface segmentation</caption>

Require: integer $K$, mesh defined by
$n_V \times 3$ array of points $V$
$n_T \times 3$ array of triangles $T$
1: Compute sparse matrices $G$ and $M$ of size $n_V \times n_V$
2: Find $K$ first eigenvectors $X(\cdot, 1), \dots, X(\cdot, K)$ of $GX = \lambda MX$
3: Apply $K$-means clustering to the $n_V \times K$ array $X$ to get an
$n_V \times 1$ array labels
4: return labels

## 2.2. Specific Tunings for Spectral Clustering of Lobes

In the previous part we exposed the general mechanism to apply spectral clustering on general surfaces. We propose here three variations in the context of cortical surface segmentation.

### 2.2.1. Removing of Non-cortical Vertices in the Clustering

First of all we want to incorporate in the framework some supervised information that would not be efficiently extracted from the data by the spectral analysis. Namely, we consider in our experiments labeled vertices corresponding to the cingular pole, as defined in Auzias et al. (2013), to be excluded from the spectral clustering. This choice is legitimated by the fact that the surface of the cingular pole cannot be considered as homologous

to the rest of the mesh, not being cortical, but rather an artificial filling of the hole created by the segmentation procedure when dissociating the two hemispheres. Furthermore, in line with the underlying hypothesis that Fourier modes may be informative for surface parcellation, it would make sense, if the clustering process had trouble dealing with this artificial surface of cingulate pole, to treat it differently.

From a computational point of view we adopt a simple heuristic: (a) we compute first the eigenfunctions of the entire mesh, (b) we exclude the regions for which we have the desired information (e.g., cingular pole), (c) we run the spectral clustering. In the following we will use the term constrained spectral clustering when the cingular pole is excluded from the segmentation process, and unconstrained spectral clustering in the other case.

### 2.2.2. Vary the Number of Clusters and Eigenfunctions

We evaluate the qualitative influence of a given Fourier mode on a resulting segmentation. Explicitly, with an anatomical model that includes 6 regions (5 lobes and the cingular pole), the unconstrained clustering requires 6 eigenvectors (including the trivial one) to define 6 clusters and at first sight, the constrained clustering would deal with only 5 eigenvectors to define the 5 unsupervised clusters. Thus, as a test for an add-on effect of a single eigenvector, we will consider only four cases: for the unconstrained approach, the number of eigenvectors is 6 or 7, and 5 or 6 for the constrained clustering.

### 2.2.3. Effect of Smoothing

The sensitivity of the spectral clustering to cortical surface deformations may be an important issue both for theoretical (e.g., which aspect of cortical shape impacts the clustering) and practical (e.g., stability of the segmentation along development or in case of focal lesion) reasons. We tested this sensitivity by deforming the cortical surface and computing again eigenfunctions and spectral clusters on the deformed surface. For general graphs, this issue of relating clustering error to graph perturbation has been addressed in the litterature of machine learning (Huang et al., 2009). But theoretical results are often limited to two clusters and assume a small perturbation. We adopted an empirical approach and explored a family of perturbed surfaces obtained by the mean curvature flow (Huisken, 1984). This geometric flow is able to smooth a folded cortical surface and it has been shown experimentally that it produces no singularities in the case of cortical surfaces (Lefèvre et al., 2013). Moreover, this process is able to mimic the backward emergence of cortical folding pattern (Lefèvre et al., 2013). Hence, the combination of mean curvature flow with spectral clustering may be seen as a first step of investigation to ensure the stability of a cortical parcellation along the brain development and its dependency to the folding level.

## 2.3. From Subject Level to Group Level Segmentation

When spectral clustering is performed for different surfaces independently there is no guarantee (1) that clusters exhibit a good reproducibility at the group level and, even if segmentation
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Flowchart of individual spectral clustering and inter-individual strategy for relabeling. The cluster matching step is achieved via the Munkres algorithm.</figcaption></figure>

are consistent, (2) that corresponding parcels are assigned to a same label. Starting from our preliminary work (Lefèvre et al., 2014) we extended the analysis of point (1) by inspecting extremal parcellations to illustrate the variability. To satisfy point (2) we proposed two different strategies:

*   Matching of individual parcellations

This relabeling strategy is illustrated on the right block of **Figure 1**. To achieve consistent labelings, the individual segmentations obtained by the spectral clustering are mapped onto a common spherical template $S^2$ where a relabeling is performed. Note that any spherical mapping approach could be suitable if it transforms the global shape of the brain in a consistent way across subjects. Even *ad hoc* methods based on affine registration might work in practice but the choice of using a spherical template is supported by the evaluation of the reproducibility (see later). In this paper we used the Freesurfer method which controls area distortions (Fischl et al., 1999), a property that will be useful in section 2.5.

For each surface *s* we have a mapping $m_s : S_s \to S^2$. A segmentation map $f_s : S_s \to \{1, ..., K\}$ can then be extended on the spherical template as a composition $g_s = f_s \circ m_s^{-1}$. Given two surfaces *s* and *s'* we can obtain a $K \times K$ table whose value at cell $(k, l)$ is given by the cardinality of the set $g_s^{-1}(k) \cap g_{s'}^{-1}(l)$. It is simply a quantification of the overlap between two clusters given by labels *k* and *l* for different surfaces.

Based on the previous table, it is possible to find a relabeling (i.e., a permutation $\sigma_s':\{1, ..., K\} \to \{1, ..., K\}$) that maximizes the degree of overlap between the two segmentations by taking surface *p* as a reference. This step is achieved via the Munkres algorithm that solves the assignment problem Kuhn (1955).

*   Group spectral clustering

Another possibility consists in applying a global clustering to all the individual eigenvectors pooled together.

Formally, let us consider surfaces $S_1, ..., S_S$ with $N_1, ..., N_S$ vertices respectively, while $X_1, ..., X_S$ are the corresponding

arrays of *K* first eigenvectors. We can build a big vector **X** in $\mathbb{R}^{N \times K}$, where $N = N_1 + ... + N_S$, as a concatenation of all the eigenvectors of the *S* surfaces. The *K*-means clustering is applied to the vector **X** which yields a big vector of labels that produces individual segmentations. This step is summarized in Algorithm 2 and in **Figure 2**.

Given an eigenfunction $\Phi$ associated to an eigenvalue $\lambda$, note that $-\Phi$ is also an eigenfunction. To overcome this sign ambiguity issue, we used a simple correlation of two eigenfunctions resampled on the common sphere and switched the sign of one eigenfunction if the correlation was negative.

**Algorithm 2 Group spectral clustering**

**Require:** integer *K*, *S* eigenvectors obtained in Algorithm 1 and stored in
*S* × 1 cell array EV

1: **Build** X=[EV{1};...;EV{S}]
N=[length(EV{1});...;length(EV{S})]

2: Apply K-means clustering to the sum(N) × K array X to get
a
sum(N) × 1 array group labels

3: **for** s=1: S **do**

4: Compute i1=sum(N(1:s-1))+1 and
i2=sum(N(1:s))

5: Build labels{s}=group labels(i1:i2)

6: **end for**

7: **return** labels

## 2.4. Evaluation Metrics

We propose here two different metrics to evaluate the quality of parcellations across a group of surfaces. The first one assesses the relationship between cluster boundaries and some anatomical
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Flowchart of group spectral clustering. Individual eigenvectors are computed for each surface and the sign ambiguity is corrected. Concatenation of eigenvectors yields a large array on which a clustering is performed. Last, the labels are back-projected on the individual surfaces.</figcaption></figure>

landmarks that have been identified by neuroanatomical experts (sulcal lines). The second metric quantifies the distance between spectral clusters and reference clusters (obtained by concatenating smaller parcels obtained with the Desikan atlas in the Freesurfer software Desikan et al., 2006).

### 2.4.1. Spectral Boundaries and Landmarks

We consider *M* anatomical landmarks that will be compared to cluster boundaries (*M* = 2 in practice). In particular, Central Sulcus (C.S.) and Parieto-Occipital Sulcus (P.O.S.) are very important landmarks since they are reliable borders between lobes. For each discrete sulcal line *L*<sub>m</sub> we can define a mean distance to the boundaries of the segmentation map *B* of a given subject:

$$D_m = \frac{1}{|\mathcal{L}_m|} \sum_{i \in \mathcal{L}_m} d(i, \mathcal{B}) \quad (3)$$

where *d*(., .) is the geodesic distance from a point to a set of points. In practice we use the Fast Marching Algorithm and the Matlab implementation of Peyré and Cohen (2008) to obtain a distance map from *B* to any vertex of the mesh. The mean distance *D*<sub>m</sub> is also expressed as a percentage of the largest geodesic distance of the cortical mesh as approximated by extremal points of the Fiedler vector (Lefèvre et al., 2012).

### 2.4.2. Comparison of Two Parcellations

One can compare two segmentations of a mesh by considering them as partitions of the *n* vertices of the template sphere and computing a distance between these partitions. Note that the two parcellation maps may have different number of labels and we will denote them as *K* and *L* in the following.

We have used the *rand distance* which is frequently employed to obtain partition distances or partition similarities (*rand index*). It has the interesting property to be a mathematical distance and to have a very fast computation time. The *rand distance* is based on the number *a* of pairs of vertices that are in the same cluster in the two partitions, and the number *b* of pairs of vertices that are in different clusters in both partitions. The values of *a* and *b* can be obtained from the contingency table *T*<sub>*k*,*l*</sub> between the two segmentations (e.g., following Denceud and Guénoche, 2006):

$$a = \frac{1}{2} \sum_{k=1}^{K} \sum_{l=1}^{L} T_{k,l} (T_{k,l} - 1) \quad (4)$$

$$b = \frac{1}{2} \sum_{k=1}^{K} \sum_{l=1}^{L} (n - T_{k,l})(n - T_{k,l} - 1) \quad (5)$$

where *T*<sub>*k*,*l*</sub> is the number of vertices that have a label *k* (*k* between 1 and *K*) in the first segmentation and a label *l* in the second segmentation (between 1 and *L*).

Then the formula for the *rand distance* is simply:

$$1 - \frac{a + b}{n(n - 1)/2} \quad (6)$$

To evaluate the overlap between each Freesurfer lobe and the corresponding spectral cluster we have chosen *Dice coefficients*, which are commonly used in neuroimaging but do not correspond to mathematical distances. The Dice between two sets *A* and *B* is given by:

$$\frac{2|A \cap B|}{|A| + |B|} \quad (7)$$
0 corresponds to no overlap, while 1 implies equality of A and B.

### 2.4.3. Consensus Parcellation

We aim at obtaining consensus maps for visualization purposes and we have used fast and simple methods.

First it is possible to define an intuitive average segmentation in K clusters by a voting strategy. At each vertex *i* of the sphere we consider the label that is present the most which can also be mathematically written as:

$$ \mathcal{F}(i) = \arg \max_{k \in [1, K]} \left| \left\{ s = 1..S/g_s(i) = k \right\} \right| \quad (8) $$

## 2.5. Statistical Evaluation

We propose a statistical test to determine if a brain parcellation (a mapping $f_s$) obtained by spectral clustering has a statistical relationship with a segmentation map obtained by an expert system, e.g., a Freesurfer parcellation.

For a given subject *s*, an expert system provides a mapping $r_s: M_s \rightarrow \{1, ..., L\}$. Then a positive value can be obtained by evaluating Dist($r_s$, $f_s$), a distance between the ground truth $r_s$ and the reference segmentation $f_s$. This value by itself does not reflect something obvious to interpret and we need to compare it to a distribution of values under a null hypothesis of independence between the reference and the automatic segmentation. This distribution can be generated by randomly rotating the maps $g_s$ defined on the sphere, provided that the spherical parameterization $m_s$ preserves areas as much as possible to avoid biases in the spatial distribution of random parcellations on the original mesh $M_s$. Eventually the distribution yields a *p*-value.

To obtain random rotations, we generate 3 × 3 matrices with coefficients following independent standard normal distributions, apply a QR decomposition, and keep the orthogonal matrices Q which are uniformly distributed (Blaser and Fryzlewicz, 2016). It is important to note that our procedure is an improvement of the one in Alexander-Bloch et al. (2017) where 3 angles from uniform distributions along each axes were chosen and introduce a bias as pointed in Blaser and Fryzlewicz (2016).

We can sum up the procedure in the following algorithm and in Figure 3.

# 3. RESULTS

## 3.1. Datasets and Pre-processings

First we tested our method on a dataset of 62 left hemispheres previously used in Auzias et al. (2013). Triangular meshes were obtained from T1-MRI through the morphologist pipeline of BrainVisa software. Central and Parieto-Occipital sulci were delineated semi-automatically by using Surfpaint (Le Troter et al., 2011). In parallel, Freesurfer was used to obtain (1) a spherical parameterization of each mesh and (2) to generate a parcellation of each surface in 35 regions based on the atlas of Desikan et al. (2006). The 35 regions were then merged in 4 lobes (frontal,

**Algorithm 3 Significance of Rand Distance**

Require: Parcellation map labels_a, reference map labels_r
spherical mesh sphM, integer Nrot

1: Compute $d=Dist labels_a, labels_r)$

2: Initialize $rot_d$ of size $Nrot \times 1$

3: **for** n=1:**Nrot do**

4: Generate a random rotation Q

5: $rlabels_a=Interpolate$ the map labels_a from sphM to Q*sphM

6: $rot_d(n)=Dist(rlabels_a, labels_r)$

7: **end for**

8: $pval$= Number of indices n in $rot_d$ such as $rot_d(n) < d$

9: **return** $pval$

parietal, temporal, occipital), insula and cingulate pole following the lobe mapping proposed on Freesurfer wiki<fnref n="1" />.

We also tested the constrained spectral clustering on a dataset of 15 fetuses previously used in Lefèvre et al. (2015). The gestational age ranges from 21 weeks to 34 weeks. T2-weighted images were acquired on axial, coronal, and sagittal planes with a half Fourier acquisition single shot turbo spin echo (HASTE) sequence. Preprocessings consist in high resolution image reconstruction, image segmentation, and mesh generation using an adaptation of BrainVisa processes. Major sulci, when present, were manually traced with Surfpaint.

Next the results are divided in five sections. In sections 3.2, 3.3, 3.4, 3.5 results involve the first dataset. Section 3.6 illustrates the stability of the spectral partitioning with respect to surface deformation and how it can be a powerful tool to parcellate developing brains as revealed by the second dataset.

## 3.2. Descriptive Analysis of Spectral Clusters

In Figure 4, we show consensus parcellations obtained with the unconstrained (top) and constrained (bottom) approaches ($K = 6$ in both cases) in parallel with a consensus parcellation in lobes obtained with Freesurfer (middle).

The comparison of unconstrained clustering and FreeSurfer reference (first two rows) showed that spectral clusters were strongly reminiscent of anatomical lobes, even with a plain unrefined strategy. Indeed there were several differences between segmentations: (a) the insular lobe was not segmented and was divided by frontal, parietal and temporal lobes, (b) the frontal lobe was divided in two parts (light and dark blue) suggestive of prefrontal and precentral subdivisions, (c) the cingulate pole was poorly delimited and even included in a large mesial region (in red) which extended on a short portion of the external face. This parcel had little anatomical correlate and clearly impaired mesial segmentation. In this view, the introduction of a constrained cingulate pole resulted in several qualitative improvements on the

<footnote n="1">¹Cortical Parcellation/Lobe Mapping:
https://surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation</footnote>
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | Distribution of Rand distances obtained through random rotations of a parcellation on the sphere. From an initial segmentation projected onto a regularly sampled sphere, random rotations are applied to generate random parcellations onto the sphere and consequently on the original surface. Rand distances are evaluated for each random parcellation which yields an empirical distribution.</figcaption></figure>

consensus map (third row of **Figure 4**) that were also present at the individual level in **Figure 6**. By definition of the method, the cingulate pole was perfectly segmented, but it resulted in a global improvement of mesial segmentation, particularly the mesial part of the parietal lobe that was then much more consistent with the one by Freesurfer (second row) while the parietal-occipital boundary remained almost unchanged. Interestingly, the better accordance with Freesurfer segmentation extended on the external face where the temporal-occipital boundary was shifted more posteriorly. Other lobes did not show any important changes.

Striking similarities between spectral segmentation and reference are present on two boundaries, between frontal and parietal lobes (dark blue, yellow), and between parietal and occipital lobes (gray, yellow). **Figure 5** illustrates the overlap between the distribution of clustering boundaries (in green) and the distribution of central sulcus and parieto-occipital sulcus (in red), both represented on an average brain obtained with Freesurfer. This figure also offers a first look on the group variability of our approach: All the boundaries are reproducible except an isolated individual variation in the mesial part of the parietal lobe.

## 3.3. Quantitative Validation of Optimal Spectral Segmentation

After this visual description, we proposed some quantitative results to emphasize striking similarities between spectral regions and brain lobes. We will now consider mostly the constrained strategy with $K = 6$. In Supplementary Information we added other quantitative results regarding the different clustering strategies (constrained/unconstrained).

### 3.3.1. Region Boundaries and Sulci

First we obtained distances between two well-defined boundaries of lobes and corresponding spectral boundaries. Based on our

definition 3, we computed two distances $D_1$ and $D_2$ for the Central Sulcus (C.S.) and the Parieto-Occipital Sulcus (P.O.S.) of each subject. The numeric values provided a first view on the variability of the segmentations as represented at the group level on **Figure 5**. In absolute values we obtained a median distance of $D_1 = 13.3$ mm for the central sulcus and $D_2 = 9.2$ mm for the parieto-occipital sulcus. Those distances are also compared to the distances obtained with Freesurfer lobes in **Table 1**. Freesurfer is very accurate on the central sulcus with a median distance of the order of magnitude of the spatial sampling and a bit less for parietal-occipital sulcus that has a more tricky and variable pattern.

When compared to the extremal dimensions of the brain, the median distances with Spanol remain small, with a mean value of 5.7% (from 1% to 11.5%) for the C.S. and 3.9% (from 0.9% to 12.5%) for the P.O.S.

### 3.3.2. Distances Between Automatic and Freesurfer Segmentation

Secondly we evaluated another metric to quantify the adequacy between spectral clusters and anatomical information. The rand distance was computed between each automatic segmentation and the corresponding one obtained by Freesurfer. The values offered another look on the variability of the spectral clusters. **Figure 6** shows the segmentations with minimal, median and maximal distances. Quantitatively the values ranged between 0.126 and 0.153 in the constrained case (**Figure 6**).

### 3.3.3. Significance of Spectral Segmentations

To go further we obtained, for each subject, 500 random rotations of Freesurfer lobes on the sphere to generate random partitions that share similar area with the initial Freesurfer segmentation. Those random maps were then used to compute random distances with the different spectral segmentations at the individual level. We then compared the partition distance to the
<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | (A) Consensus segmentation with the unconstrained spectral approach. (B) Consensus segmentation of Freesurfer lobes. (C) Consensus segmentation with the constrained spectral approach.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**FIGURE 5** | In green, distribution of the clustering boundaries represented on an average brain. In red, distribution of central sulcus (Left) and parieto-occipital sulcus (Right). The overlap between the two spatial distribution is encoded by a mix of the two colors.</figcaption></figure>

randomly generated distribution of distances. We counted how many subjects among the 62 had a *p*-value below a significance level of 0.01 as obtained from the statistical procedure in 2.5. Constraining the cingulate pole for *K* = 6 improved the ratio of subjects showing a statistical relationship with brain lobes (from 98.4 to 100%).

### 3.3.4. Optimality Regarding the Number of Clusters

We computed partition distances between constrained spectral approaches with *K* labels and *K* eigenvectors when *K* varies

between 3 and 10, given that *K* = 2 is pointless (two regions with the cingulate pole and its complement part). On **Figure 7** (left), we observed the distribution of rand distances and the median values (red) ranged between 0.15 and 0.3 with a clear minimum for *K* = 6. On **Figure 7** (right), we added information obtained through our random rotations procedure. We displayed both the distribution of individual *p*-values (boxplot) and the ratio of subjects for which the rand distance was significant with respect to the random distribution. This ratio was 1 for *K* = 6 and 0.38 for *K* = 7 (24 subjects). For all the other values, the
<figure><img src="image_8.png" /><figcaption>**FIGURE 6** | Three spectral segmentations with constrained cingulate pole (first row) for which the distance with Freesurfer (second row) is the smallest (A), median (B), and the largest (C).</figcaption></figure>

<table><caption>**TABLE 1** | Median distances associated to central sulcus and parieto-occipital sulcus for Spanol and FreeSurfer segmentations.</caption><thead><tr><th></th><th>Spanol</th><th>FreeSurfer</th></tr></thead><tbody><tr><td>Central</td><td>13.3 mm (2.2–28)</td><td>1.4 mm (0.8–2.5)</td></tr><tr><td>Parieto-Occipital</td><td>9.2 mm (2–27.9)</td><td>4.3 mm (1.6–24.6)</td></tr></tbody></table>

ratio was below 0.3 suggesting poor association between spectral segmentations and Freesurfer lobes. On Figure 2SI we showed all the consensus segmentation obtained when K varied between 3 and 10.

## 3.4. Group Analysis of Spectral Clustering

In this part we further analyzed the ability of the spectral clustering approaches to delineate automatically brain lobes. We evaluated first the overlaps between each Freesurfer lobe and the corresponding spectral region and then by comparing the two strategies described in section 2.3. Those steps necessitated before a consistent matching between the regions, at the subject and group level, by using the hungarian algorithm.

### 3.4.1. Lobe by Lobe Statistics

Beyond the global rand distances shown in S.I., we computed Dice coefficients between each Freesurfer lobe and the corresponding one in two spectral approaches (unconstrained, constrained, with 6 eigenvectors). The corresponding frontal lobe was obtained by concatenating two spectral regions (light and dark blue in most of the previous figures). The Insula was not included since it had no corresponding region after the grouping strategy for the frontal lobe. The results were summarized in Table 2. For all lobes except for the occipital one (almost no change), the mean Dice coefficients were improved by constraining the cingulate pole, particularly for the parietal lobe.

### 3.4.2. Comparison of the Two Group Strategies

In Table 3 we showed the 6 Dice coefficients between the regions obtained with the individual spectral clustering and the joint group clustering. We remarked that dice values were large (above 0.93) and close to the perfect overlap (= 1). For all lobes except

the temporal one, the overlap was improved when the constraint on the cingular pole was added.

For the subject with minimal overlap the two lowest Dice coefficients were 0.88 for both temporal and parietal lobe (see SI for the visual comparison).

## 3.5. Features of the Spectral Clustering

For each subject, the 6 first non trivial eigenfunctions were projected onto a brain template by using Freesurfer. On Figure 8 we displayed the mean pattern of each eigenvector used in the spectral clustering. For each eigenfunction the variability is illustrated by showing the nodal lines across the 62 subjects.

It is remarkable to observe that the individual nodal lines are distributed along consistent directions. The variability of those nodal lines is simply represented by a shift transverse to the main direction. In a dual perspective, one observes for each eigenfunction a clear oscillating pattern with several connected regions sharing the same sign (so called nodal domains). For each subject, there are respectively 2, 2, 2, 2, 3, 3 nodal domains for the 6 first eigenfunctions in increasing order. Those results illustrate the perfect stability of eigenfunction from a topological point of view.

## 3.6. Robustness of the Constrained Spectral Segmentation

### 3.6.1. Effect of Surface Smoothing

For each subject we applied the mean curvature flow to smooth the folded geometry of the cortical meshes. For several iterations of the process ($t = 1, 10, 20, 50, 100, 200, 300, 500, 750, 1,000$) we ran the constrained spectral clustering procedure with 6 eigenvectors and computed the partition distance between the resulting segmentation and the initial segmentation. The variability of this distance was illustrated in Figure 9 (left). We observed a regular increase of the mean value till an interval [0.04 – 0.09] and a spreading of the distribution as well. Compared to the rand distances found between constrained spectral parcellations and Freesurfer lobes (around 0.15 on Figure 7 or Table 2SI) which reflected good agreements, one can therefore consider that the perturbation of spectral parcels during the smoothing process is small.
<figure><img src="image_9.png" /><figcaption>**FIGURE 7** | (Left) Rand distances between constrained spectral segmentation with *K* regions and Freesurfer segmentation with 6 lobes. (Right) *p*-values for each subject represented as boxplots and proportion of subjects (blue bars) for which the *p*-value is < 0.01.</figcaption></figure>

<table><caption>**TABLE 2** | Dice coefficients between Freesurfer lobes and each of the 5 corresponding regions obtained with the spectral segmentation.</caption><thead><tr><th>Cingulate</th><th>Frontal</th><th>Parietal</th><th>Temporal</th><th>Occipital</th></tr></thead><tbody><tr><td>0.39 ± 0.04</td><td>0.89 ± 0.01</td><td>0.66 ± 0.05</td><td>0.82 ± 0.02</td><td><strong>0.82 ± 0.03</strong></td></tr><tr><td><strong>0.91 ± 0.02</strong></td><td><strong>0.94 ± 0.01</strong></td><td><strong>0.83 ± 0.03</strong></td><td><strong>0.87 ± 0.01</strong></td><td>0.81 ± 0.03</td></tr></tbody></table>

Note that the corresponding frontal lobe is obtained by concatenating two spectral regions. The Insula is not included since it has no corresponding region after the grouping strategy in frontal lobe. First row: unconstrained method with 6 eigenvectors. Second row: constrained method with 6 eigenvectors. The bold value indicate the largest value between the two rows.

<table><caption>**TABLE 3** | Dice coefficients for each of the 6 spectral regions obtained with the individual and the group strategy.</caption><thead><tr><th>Front</th><th>PreFront</th><th>Occip</th><th>Pariet</th><th>Tempo</th><th>Mesial</th></tr></thead><tbody><tr><td>0.95 ± 0.03</td><td>0.93 ± 0.04</td><td>0.96 ± 0.03</td><td>0.94 ± 0.03</td><td><strong>0.97 ± 0.02</strong></td><td>0.97 ± 0.02</td></tr><tr><td><strong>0.96 ± 0.02</strong></td><td><strong>0.98 ± 0.01</strong></td><td><strong>0.97 ± 0.02</strong></td><td>0.94 ± 0.03</td><td>0.96 ± 0.03</td><td><strong>1.00 ± 0.00</strong></td></tr></tbody></table>

First row: unconstrained method with 6 eigenvectors. Second row: constrained method with 6 eigenvectors. The bold value indicate the largest value between the two rows.

In **Figure 9** (right) we observed the constrained spectral clustering for a surface whose rand distance correspond to the median value (0.06) of the distance distributions for 1000 iterations. The most visible change can be seen on the mesial face with a slight bifurcation of the yellow/blue boundary that has a counterpart on the external face close to the central sulcus.

### 3.6.2. Segmentation of the Immature Brain

In **Figure 10** we showed the constrained spectral segmentation during early development on 15 fetuses whose gestational age ranges from 21 weeks to 34 weeks. We observed an informative consistency between the different regions across ages. Moreover the position of the parcells was qualitatively similar to what was observed for adult brains. But more precisely the analogous of frontal/parietal boundary crossed the central sulcus which was not the case for adult brains as revealed in particular by **Figure 5**. It is also possible to quantify the distance between C.S. (respectively P.O.S) and its corresponding boundary, except for

the two youngest subjects (21 and 24 weeks gestational age). For the C.S. the absolute values range from 0.7 mm (28 w GA) to 6.3 mm (25 w GA) while for the P.O.S it was from 1.3 mm (25 w GA) to 6.1 mm (33 w GA). Regarding the percentage of the largest distance, the mean value is 2.1% (from 0.9 to 8.1%) for C.S., and 3.44% (from 0.0 to 6.2%) for the P.O.S. In both cases the correlation of those values with age were found non significant ($R = 0.46$ and $R = -0.07$ respectively) when correcting for the increasing global size of the cortical surface, as measured by its longest geodesic (Lefèvre et al., 2012).

## 4. DISCUSSION AND CONCLUSION

The major applicative result of our general Spanol approach is the possibility to segment a cortical surface in a few parcels that have strong similarities with brain lobes. These similarities can be precisely quantified by partition distances and a dedicated statistical procedure, which has been proposed to determine how far the resulting parcels stood from randomly segmented lobe-sized parcels. The segmentation process can be entirely unsupervised, using only the geometrical properties of the hemispheric mesh shape. However, we improved dramatically the global overlap with the four main lobes (frontal, parietal, temporal and occipital) obtained by Freesurfer software by adding limited constraints only, excluding points that belong to non-cortical regions of the mesh, namely the cingulate pole. This could be seen as a semi-supervised clustering at the level of the whole mesh, but remains strictly unsupervised regarding the cortical level. We would insist on two main neurological relevancies of Spanol, and further discuss interesting methodological considerations raised by its development.

### 4.1. Practical Aspects

First, to our knowledge it is the first attempt to extract in a fast, simple, reproducible, and almost unsupervised manner a parcellation of a cortical surface in lobes. This approach does not require any atlas and can be applied as it is to triangular meshes, at the individual level, in a purely intrinsic way. The computation time is limited to a few seconds for each surface.
<figure><img src="image_10.png" /></figure>

**FIGURE 8** | Center: Mean pattern of the 6 first non trivial eigenfunctions. Positive values are in red, negative values in dark blue and zero values in green with the nodal line in black. Around each eigenfunction is represented the variability of nodal lines across the 62 subjects.

<figure><img src="image_11.png" /></figure>

**FIGURE 9** | (Left) Distribution of rand distances between initial constrained segmentations and segmentations for increasing number of iterations in the smoothing process. (Right) Illustration of the smoothing process for a surface whose rand distance correspond to the median value (0.06) of the distance distributions for 1,000 iterations. 3 iterations are considered (0, 100, and 1,000).

<figure><img src="image_12.png" /></figure>

**FIGURE 10** | Constrained spectral segmentations of fetal brains with preserved proportions and gestational age in weeks. The central sulcus, when present, is superimposed in red.
Our unsupervised spectral clustering approach identifies boundaries that may be less accurate for C.S. and P.O.S than the ones obtained by FreeSurfer in a pure supervised manner. The interpretation of the larger discrepancy with our spectral approach clearly depends on the expectancies one can have for a new totally unsupervised approach by comparison to a long-optimized multi-input one. The main difference is on the central boundary, for which the spectral approach generates a kind of global posterior shifting, even if the limit can be really accurate in some subject. This interesting point requires further investigation that we are currently conducting.

Although some local improvements can be expected according to experimental needs, they may be achieved with complementary segmentation steps, either supervised or not, following Spanol backbone. For instance, the strict respect of bottom sulci lines along some lobar boundaries (e.g., for C.S. or P.O.S) would be a possible finishing stage, in the same spirit as Klein and Tourville (2012). It could be also interesting to consider this segmentation problem in a semi-supervised manner following for instance (Lu and Carreira-Perpinán, 2008). The exclusion of otherwise segmented Insula before spectral clustering has been also tested without disturbing the relevance of lobar boundaries. Indeed, the implementation of Spanol could allow regional analysis of morphometric parameters, in the spirit of Toro et al. (2009), to be rapidly performed without requiring complicated pipelines and potentially skewing normalization steps.

## 4.2. Link With Neurodevelopmental Questions

Second and more fundamental, spectral lobes are obtained from the intrinsic geometrical informations contained in the low-frequency eigenfunctions only. Yet, the lobar limits given by Spanol are not only consistent with the classical sulcal ones (C.S., P.O.S) but also for the ill-defined one such as the lateral parietal occipital or the temporal occipital. Even the segmentation of a prefrontal sub-lobar region seems rather consistent. While classical anatomy relies on almost *ad hoc* landmarks such as the preoccipital notch, or locally defined virtual lines (e.g., occipital temporal limit), spectral analysis kind of intergrade the whole shape into objectively defined boundaries. Indeed, the fact that lobes can be intrinsically defined from the shape of the brain is *per se* an intriguing result. Interestingly, not all the geometrical information is required, but only the one of the low-frequency eigenfunctions that give mostly information on the global shape of the brain (Seo and Chung, 2011; Germanaud et al., 2012). It may suggest a strong relationship between global shape of the cortical surface and the appearance of great primary sulci that defined main lobar limits (C.S., P.O.S), but also with functional aspects associated with lobes. Anyway, this dependency upon global shape only would explain the stability of spectral segmentation of lobes along development regardless of the full expansion of the main sulci. This claim has been empirically verified by smoothing the cortical surfaces to decrease the influence of sulci and gyri in the geometry and by quantifying the stability of spectral lobes during this deformation.

Moreover the good reproducibility of spectral lobes on a small population of fetuses with various folding complexity is another argument in favor of a very early determination of brain lobes, even from a very smooth cortical surface. This question is also closely intertwined with the emergence of primary folds during fetal development. Among various explanations of this phenomenon, the mecanical hypothesis (Tallinen et al., 2016) is seducing in our case because of oscillatory properties of the surface as proxies of its ability to fold and determine the position of primary folds. Indeed, the question of the relationship between shape and function is all the more interesting since one inverts the classical paradigma and addresses the fact that the determinant of shape may impact the functional outcome (Foubet and Toro, 2015). So far, it is a very general question, not restricted to the lobar one.

From a more practical point of view, Spanol offers a promising way to define automatically brain lobes in fetuses or newborns because of this nice continuity in the ontogeny. It would remain to test more systematically our method in larger longitudinal databases with complementary information provided by cytoarchitecture, connectivity or functions.

## 4.3. Final Remarks

Spanol provides a segmentation of cortical surface by using a *K*-means clustering of Laplace-Beltrami operator eigenfunctions. In other terms, only few low frequency descriptors of the brain geometry are required (optimum with 6 in our work) to provide a given number of relevant connected regions on the cortical surface. This method has been used at the individual level and at the group level by concatenating all the eigenfunctions of each subject of the group. A major advantage of the second approach is the ability to obtain directly a consistent labeling of the resulting parcels. In the individual case, it requires to solve the assignment problem to match parcels of different subjects. From a computational point of view the group spectral clustering with 62 meshes of 100 – 150k vertices takes approximately 1 min on a standard laptop which is comparable to 62 individual *K* means (*K* = 6 in our case). The results showed strong similarities between the segmentations obtained at the individual or group level as measured by the rand distance. This result legitimates the group spectral clustering on our data as a way to segment individual meshes in a consistent, reliable and fast manner. It also suggests a very reproducible pattern in the low frequency eigenfunctions of cortical surfaces which has already been pointed out in several works (Germanaud et al., 2012; Lombaert et al., 2013b; Lefèvre and Auzias, 2015). More precisely in Lombaert et al. (2013a) the authors observe that between 5 and 15 eigenvectors perform equally well in terms of brain regions overlap when used in their spectral registration algorithm. Besides it could be interesting to explore whether a group spectral clustering procedure like our's could be applied successfully to a larger number of regions. Indeed, in the field of computer graphics the co-segmentations of shapes with spectral approaches have been shown efficient but often tested for limited number of components (Sidi et al., 2011). Conversely studies with other kinds of neuroimaging data (fMRI) suggested that spectral
clustering approaches would not be satisfactory for large values of K (Thirion et al., 2014).

A novelty in our approach concerns the statistical procedure to test whether a given brain segmentation could be considered as randomly sampled from a null distribution of brain parcellations sharing similar properties. We escaped the intractable (and biologically irrealistic) combinatoric of testing all possible partitions of a triangular mesh by proposing random rotations of a segmentation map onto a reference sphere and computing distances with a reference segmentation. This approach is able to find statistical associations between spatial partitions independently of the partition distance used and of the disparity between the number of clusters in the reference and the automatic segmentations. This procedure takes approximately a few minutes for each subject for 500 random samples. With this statistical tests, we are able to determine which subjects have a spectral parcellation that is significantly close to a reference one. At the end we can find the number of clusters which maximizes the number of significant subjects.

Regarding the supervised information injected in fixing the cingulate pole, we have adopted a strategy where eigenfunctions of the original surface are used in the spectral clustering while excluding the points of the constraint. Something more mathematically natural would have been to consider the eigenfunctions of the Laplace-Beltrami Operator with boundary conditions on the borders of the cingular pole. But our choice was motivated by purely pragmatic considerations regarding unsatisfactory segmentations in the second case. Why the sound mathematical framework of Neuman or Dirichlet boundary conditions has failed and why our *ad-hoc* treatment of boundaries in the cingulate pole is more successful, remains an interesting question that might even be more general than in our specific context of neuroanatomy.

As closing remarks, our work illustrates, from the applicative point of view of brain anatomy, one of the fascinating properties of the Laplace-Beltrami Operator. For a long

time low-frequency patterns of its first eigenfunctions have been well described by theorems like Courant's nodal one, which estimates the number of spatial oscillations. Despite a few theoretical results, spectral clustering remains a rather empirical way to combine the geometrical information contained in the eigenfunctions. Still, our results show its ability to unveil structure in datasets such as those in neuroimaging. Besides we rediscovered after (Jin et al., 2005) that varying slightly the number of eigenvectors can improve clustering results. This discrepancy with the idealized view of spectral clustering offers also stimulating theoretical perspectives in the future.

## ETHICS STATEMENT

All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Ethical Committee of Aix-Marseille University.

## AUTHOR CONTRIBUTIONS

JL and DG: designed research; JL, AP, JM: conducted analyses; FD, NG: provided and processed fetal data; JL, GA, and DG: wrote the paper. All authors took part in the scientific discussion at multiple stages of the study.

## ACKNOWLEDGMENTS

This work has been funded by the French Agency for Research (ANR-12-JS03-001-01, Modegy).

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2018.00354/full#supplementary-material

## REFERENCES

* Alexander-Bloch, A., Vandekar, S. N., Shinohara, R. T., Liu, S., Satterthwaite, T. D., Glahn, D. C., et al. (2017). The correspondence problem: which brain maps are significantly similar? *bioRxiv* 203083. doi: 10.1101/203083

* Auzias, G., Brun, L., Deruelle, C., and Coulon, O. (2015). Deep sulcal landmarks: algorithmic and conceptual improvements in the definition and extraction of sulcal pits. *Neuroimage* 111, 12–25. doi: 10.1016/j.neuroimage.2015.02.008

* Auzias, G., Coulon, O., and Brovelli, A. (2016). Marsatlas: a cortical parcellation atlas for functional mapping. *Hum. Brain Mapp.* 37, 1573–1592. doi: 10.1002/hbm.23121

* Auzias, G., Lefevre, J., Le Trotter, A., Fisher, C., Perrot, M., Régis, J., et al. (2013). Model-driven harmonic parameterization of the cortical surface: hip-hop. *IEEE Trans. Med. Imaging* 32, 873–887. doi: 10.1109/TMI.2013.2241651

* Belkin, M., and Niyogi, P. (2008). Towards a theoretical foundation for Laplacian-based manifold methods. *J. Comput. Syst. Sci.* 74, 1289–1308. doi: 10.1016/j.jcss.2007.08.006

* Berger, M. (2003). *A Panoramic View of Riemannian Geometry*. Berlin; Heidelberg: Springer Verlag.

* Blaser, R., and Fryzlewicz, P. (2016). Random rotation ensembles. *J. Mach. Learn. Res.* 17, 1–26. Available online at: http://www.jmlr.org/papers/v17/blaser16a.html

* Cachia, A., Mangin, J.-F., Rivière, D., Papadopoulos-Orfanos, D., Kherif, F., Bloch, I., et al. (2003). A generic framework for the parcellation of the cortical surface into gyri using geodesic voronoi diagrams. *Med. Image Anal.* 7, 403–416. doi: 10.1016/S1361-8415(03)00031-8

* Chladni, E. F. F. (1787). *Entdeckungen über die Theorie des Klanges*. Leipzig: Zentralantiquariat der DDR.

* Denœud, L., and Guénoche, A. (2006). “Comparison of distance indices between partitions,” in *Data Science and Classification* (Ljubljana: Springer Berlin Heidelberg), 21–28.

* Desikan, R. S., Ségonne, F., Fischl, B., Quinn, B. T., Dickerson, B. C., Blacker, D., et al. (2006). An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest. *Neuroimage* 31, 968–980. doi: 10.1016/j.neuroimage.2006.01.021

* Fischl, B., Sereno, M. I., and Dale, A. M. (1999). Cortical surface-based analysis: II: inflation, flattening, and a surface-based coordinate system. *Neuroimage* 9, 195–207. doi: 10.1006/nimg.1998.0396

* Foubet, O., and Toro, R. (2015). Mechanical morphogenesis and the development of neocortical organisation. *bioRxiv* 021311.
Friston, K. (2002). Beyond phrenology: what can neuroimaging tell us about distributed circuitry? *Annu. Rev. Neurosci.* 25, 221–250. doi: 10.1146/annurev.neuro.25.112701.142846

Germanaud, D., Lefèvre, J., Toro, R., Fischer, C., and Dubois, J. (2012). Larger is twistier: spectral analysis of gyration (spangy) applied to adult brain size polymorphism. *NeuroImage* 63, 1257–1272. doi: 10.1016/j.neuroimage.2012.07.053

Glasser, M. F., Smith, S. M., Marcus, D. S., Andersson, J. L., Auerbach, E. J., Behrens, T. E., et al. (2016). The human connectome project's neuroimaging approach. *Nat. Neurosci.* 19, 1175–1187. doi: 10.1038/nn.4361

Gratiolet, L. P. (1854). *Mémoire sur les Plis cérébraux de l'Homme et des Primates: Mit Einem Atlas (4 pp. XIV pl.) in fol.* 33i. Paris: A. Bertrand.

Huang, L., Yan, D., Taft, N., and Jordan, M. I. (2009). “Spectral clustering with perturbed data,” in *Advances in Neural Information Processing Systems*, eds D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou (Vancouver, Curran Associates Inc.), 705–712.

Huisken, G. (1984). Flow by mean curvature of convex surfaces into spheres. *J. Diff. Geometry* 20, 237–266. doi: 10.4310/jdg/1214438998

Jin, R., Kang, F., and Ding, C. H. (2005). “A probabilistic approach for optimizing spectral clustering,” in *Advances in Neural Information Processing Systems* (Vancouver, CA), 571–578.

Klein, A., and Tourville, J. (2012). 101 labeled brain images and a consistent human cortical labeling protocol. *Front. Neurosci.* 6:171. doi: 10.3389/fnins.2012.00171

Kuhn, H. W. (1955). The hungarian method for the assignment problem. *Naval Res. Logist.* Q. 2, 83–97. doi: 10.1002/nav.3800020109

Lai, R., Shi, Y., Dinov, I., Chan, T., and Toga, A. (2009). “Laplace-beltrami nodal counts: a new signature for 3d shape analysis,” in *IEEE International Symposium on Biomedical Imaging: From Nano to Macro, ISBI'09* (Boston, MA), 694–697.

Le Troter, A., Rivière, D., and Coulon, O. (2011). “An interactive sulcal fundi editor in brainvisa,” in *Proceedings of the Organization for Human Brain Mapping Conference* (Quebec City, QC), 61.

Lefèvre, J., and Auzias, G. (2015). “Spherical parameterization for genus zero surfaces using laplace-beltrami eigenfunctions,” in *Proceedings of the Second International Conference, GSI 2015, October 28-30, 2015*, Vol. 9389. (Palaiseau).

Lefèvre, J., Auzias, G., and Germanaud, D. (2014). “Brain lobes revealed by spectral clustering,” in *Pattern Recognition (ICPR), 2014 22nd International Conference on* (Stockholm: IEEE), 562–567.

Lefèvre, J., Germanaud, D., Dubois, J., Rousseau, F., de Macedo Santos, I., Angleys, H., et al. (2015). Are developmental trajectories of cortical folding comparable between cross-sectional datasets of fetuses and preterm newborns? *Cereb. Cortex* 26, 3023–3035. doi: 10.1093/cercor/bhv123

Lefèvre, J., Germanaud, D., Fischer, C., Toro, R., Riviere, D., and Coulon, O. (2012). “Fast surface-based measurements using first eigenfunction of the laplace-beltrami operator: interest for sulcal description,” in *Biomedical Imaging (ISBI), 2012 9th IEEE International Symposium on* (Barcelona: IEEE), 1527–1530.

Lefèvre, J., Intwali, V., Hertz-Pannier, L., Hüppi, P. S., Mangin, J.-F., Dubois, J., et al. (2013). “Surface smoothing: a way back in early brain morphogenesis,” in *International Conference on Medical Image Computing and Computer-Assisted Intervention* (Nagoya: Springer), 590–597.

Lefèvre, J., and Mangin, J. (2010). A reaction-diffusion model of human brain development. *PLoS Comput. Biol.* 6:e1000749. doi: 10.1371/journal.pcbi.1000749

Lefranc, S., Roca, P., Perrot, M., Poupon, C., Le Bihan, D., Mangin, J.-F., et al. (2016). Groupwise connectivity-based parcellation of the whole human cortical surface using watershed-driven dimension reduction. *Med. Image Anal.* 30, 11–29. doi: 10.1016/j.media.2016.01.003

Liu, R., and Zhang, H. (2004). “Segmentation of 3d meshes through spectral clustering,” in *Computer Graphics and Applications, 2004. PG 2004. Proceedings. 12th Pacific Conference on* (Seoul: IEEE), 298–305.

Lohmann, G., and von Cramon, D. Y. (2000). Automatic labelling of the human cortical surface using sulcal basins. *Med. Image Anal.* 4, 179–188. doi: 10.1016/S1361-8415(00)00024-4

Lombaert, H., Arcaro, M., and Ayache, N. (2015). “Brain transfer: spectral analysis of cortical surfaces and functional maps,” in *International Conference on Information Processing in Medical Imaging* (Sabhal Mor Ostaig: Springer), 474–487.

Lombaert, H., Grady, L., Polimeni, J. R., and Cheriet, F. (2013a). Focus: feature oriented correspondence using spectral regularization–a method for

precise surface matching. *IEEE Trans. Patt. Anal. Mach. Intell.* 35, 2143–2160. doi: 10.1109/TPAMI.2012.276

Lombaert, H., Sporring, J., and Siddiqi, K. (2013b). “Diffeomorphic spectral matching of cortical surfaces,” in *International Conference on Information Processing in Medical Imaging* (Asilomar: Springer), 376–389.

Lu, Z., and Carreira-Perpinhán, M. A. (2008). “Constrained spectral clustering through affinity propagation,” in *Computer Vision and Pattern Recognition, 2008. CVPR 2008. IEEE Conference on* (Anchorage, AK: IEEE), 1–8.

Ng, A. Y., Jordan, M. I., and Weiss, Y. (2002). On spectral clustering: analysis and an algorithm. *Adv. Neural Inform. Process. Syst.* 2, 849–856. doi: 10.1.1.19.8100

Parisot, S., Arslan, S., Passerat-Palmbach, J., Wells, W. M., and Rueckert, D. (2016). Group-wise parcellation of the cortex through multi-scale spectral clustering. *NeuroImage* 136, 68–83. doi: 10.1016/j.neuroimage.2016.05.035

Pepe, A., Auzias, G., Guio, F. D., Rousseau, F., Germanaud, D., Mangin, J.-F., et al. (2015). “Spectral clustering based parcellation of FETAL brain MRI,” in *12th IEEE International Symposium on Biomedical Imaging, ISBI 2015*, April 16-19, 2015 (Brooklyn, NY), 152–155.

Peyré, G., and Cohen, L. D. (2008). Heuristically driven front propagation for fast geodesic extraction. *Int. J. Comput. Vis. Biomech.* 1, 55–67. doi: 10.1007/11567646_15

Rabiei, H., Richard, F., Roth, M., Anton, J.-L., Coulon, O., and Lefèvre, J. (2015). “The graph windowed fourier transform: a tool to quantify the gyration of the cerebral cortex,” in *Workshop on Spectral Analysis in Medical Imaging (SAMI)*. (Munich).

Seo, S., and Chung, M. K. (2011). “Laplace-beltrami eigenfunction expansion of cortical manifolds,” in *Biomedical Imaging: From Nano to Macro, 2011 IEEE International Symposium on* (Chicago, IL: IEEE), 372–375.

Sidi, O., van Kaick, O., Kleiman, Y., Zhang, H., and Cohen-Or, D. (2011). “Unsupervised co-segmentation of a set of shapes via descriptor-space spectral clustering,” in *ACM Transactions on Graphics (TOG)*, Vol. 30 (New York, NY: ACM), 126.

Tallinen, T., Chung, J. Y., Rousseau, F., Girard, N., Lefèvre, J., and Mahadevan, L. (2016). On the growth and form of cortical convolutions. *Nat. Phys.* 12, 588–593. doi: 10.1038/nphys3632

Thirion, B., Flandin, G., Pinel, P., Roche, A., Ciuciu, P., and Poline, J.-B. (2006). Dealing with the shortcomings of spatial normalization: multi-subject parcellation of fMRI datasets. *Hum. Brain Mapp.* 27, 678–693. doi: 10.1002/hbm.20210

Thirion, B., Varoquaux, G., Dohmatob, E., and Poline, J.-B. (2014). Which fMRI clustering gives good brain parcellations? *Front. Neurosci.* 8:167. doi: 10.3389/fnins.2014.00167

Toro, R., Chupin, M., Garnero, L., Leonard, G., Perron, M., Pike, B., et al. (2009). Brain volumes and val66met polymorphism in the bdnf gene: local or global effects? *Brain Struct. Funct.* 213, 501–509. doi: 10.1007/s00429-009-0203-y

Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., et al. (2002). Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. *Neuroimage* 15, 273–289. doi: 10.1006/nimg.2001.0978

Von Luxburg, U. (2007). A tutorial on spectral clustering. *Stat. Comput.* 17, 395–416. doi: 10.1007/s11222-007-9033-z

Wachinger, C., Golland, P., Kremen, W., Fischl, B., Reuter, M., Initiative, A. D. N., et al. (2015). Brainprint: a discriminative characterization of brain morphology. *NeuroImage* 109, 232–248. doi: 10.1016/j.neuroimage.2015.01.032

Zilles, K., and Amunts, K. (2010). Centenary of brodmann's map. Conception and fate. *Nat. Rev. Neurosci.* 11, 139–145. doi: 10.1038/nrn2776

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Lefèvre, Pepe, Muscato, De Guio, Girard, Auzias and Germanaud. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.