Hamed Rabiei

Frédéric Richard

Olivier Coulon

Julien Lefèvre

Local Spectral Analysis of the Cerebral Cortex: New Gyrification Indices

Keywords: Brain folding, gyrification index, spectral analysis, mean curvature, windowed Fourier transform, allometric relation

Gyrification index (GI) is an appropriate measure to quantify the complexity of the cerebral cortex. There is, however, no universal agreement on the notion of surface complexity and there are various methods in literature that evaluate different aspects of cortical folding. In this paper, we give two intuitive interpretations on folding quantification based on the magnitude and variation of the mean curvature of the cortical surface. We then present a local spectral analysis of the mean curvature to introduce two local gyrification indices that satisfy our interpretations. For this purpose, the graph windowed Fourier transform is extended to the framework of surfaces discretized with triangular meshes. An adaptive window function is also proposed to deal with the intersubject cortical size variability. The intrinsic nature of the method allows us to compute the degree of folding at different spatial scales. Our experiments show that while more classical surface area-based GIs may fail at differentiating deep folds from very convoluted ones, our spectral GIs overcome this issue. The method is applied to the cortical surfaces of 124 healthy adult subjects of OASIS database and average gyrification maps are computed and compared with other GI definitions. In order to illustrate the capacity of our method to capture and quantify important aspects of gyrification, we study the relationship between brain volume and cortical complexity, and design a scaling analysis with a power law model. Results indicate an allometric relation and confirm the well-known observations that larger brains are more folded. We also perform the scaling analysis at the vertex level to investigate how the degree of folding varies locally with the brain volume. Results reveal that in our healthy adult brain database, cortical regions which are the least folded on average show an increased folding complexity when brain size increases.

I. INTRODUCTION

The process of folds appearing on the cerebral cortex is called gyrification. It occurs mostly during the second half of fetal life but continues changing slightly the shape of the folds during post-natal life [START_REF] Armstrong | The ontogeny of human gyrification[END_REF]- [START_REF] Lefèvre | Are developmental trajectories of cortical folding comparable between crosssectional datasets of fetuses and preterm newborns?[END_REF]. Morphometric parameters of the cerebral cortex such as volume, surface area, sulcal depth, curvature and gyrification index (GI) are commonly used to quantify this gyrification process. Among these, GI is a direct measure that attempts to quantify the degree of folding of the cerebral cortex. It provides valuable information about changes occurring on a brain surface during development, aging, and disease [START_REF] Lefèvre | Are developmental trajectories of cortical folding comparable between crosssectional datasets of fetuses and preterm newborns?[END_REF]- [START_REF] Shimony | Comparison of cortical folding measures for evaluation of developing human brain[END_REF]. In this paper, we propose two new GIs that take into account two different aspects of folding quantification.

Existing GI definitions in the literature, from a general point of view, may be categorized in two classes: surface area (perimeter)-based methods [START_REF] Toro | Brain size and folding of the human cerebral cortex[END_REF]- [START_REF] Li | Mapping longitudinal development of local cortical gyrification in infants from birth to 2 years of age[END_REF], [START_REF] Zilles | The human pattern of gyrification in the cerebral cortex[END_REF]- [START_REF] Su | Geometric computation of human gyrification indexes from magnetic resonance images[END_REF] and curvaturebased methods [START_REF] Shimony | Comparison of cortical folding measures for evaluation of developing human brain[END_REF], [START_REF] Luders | A curvature-based approach to estimate local gyrification on the cortical surface[END_REF]- [START_REF] Shishegar | Quantifying gyrification using Laplace Beltrami eigenfunction level-sets[END_REF].

Methods in the first category compute the gyrification index as a ratio between the local area (perimeter) of the cortical surface and that of a reference surface. For instance, Toro et al. [START_REF] Toro | Brain size and folding of the human cerebral cortex[END_REF] defined a local GI as a ratio between the area of the pial surface contained in a spherical neighbourhood of each cortical point and the area of the great disc of the sphere. Schaer et al. [START_REF] Schaer | A surface-based approach to quantify local cortical gyrification[END_REF] proposed a ratio between the area of a region of interest, determined by intersection of a sphere with the convex hull Copyright (c) 2010 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to pubs-permissions@ieee.org.

Corresponding author: H. Rabiei (email: hamed.rabiei@univ-amu.fr).

of the cortical surface, and that of the corresponding patch on the pial surface as a local GI. The locality of results of those methods depends on the radius of the spherical neighbourhood.

In order to deal with the size variability in a database (e.g. fetal brains), the radius would have to be adapted for each subject but no intrinsic strategy exists to choose an appropriate one.

To tackle this issue, Lefèvre et al. [START_REF] Lefèvre | Are developmental trajectories of cortical folding comparable between crosssectional datasets of fetuses and preterm newborns?[END_REF] proposed to adapt the radius r so that it equals a fraction of the brain length in the rostro-caudal direction. Another drawback of those methods is that if there are two cortical regions with equal areas inside the sphere, one contains deep folds and the other one is formed by shallower but more oscillating folds, these methods give equal GI value to both regions. This means that they may not distinguish between deep folds and oscillating folds with equal surface areas. For instance, Schaer's and Toro's methods both give high GI values to deep folds like the insula and the central sulcus. In a more general way, both methods tend to produce similar GI maps; see [5, Fig. 5a-c] and [START_REF] Schaer | A surface-based approach to quantify local cortical gyrification[END_REF]Fig. 4]. Methods in the second category rely on the curvature and its derivatives e.g. the mean curvature. The mean curvature is a geometric tool that measures locally how a surface is deviated from being flat. It defines a function on the cerebral cortex and assigns positive values to points on gyri and negative values to points on sulci [START_REF] Meyer | Discrete differentialgeometry operators for triangulated 2-manifolds[END_REF]. The mean curvature map, however, is too local to deliver a helpful insight into the surface folding [START_REF] Luders | A curvature-based approach to estimate local gyrification on the cortical surface[END_REF]. Moreover, as Shimony et al. [START_REF] Shimony | Comparison of cortical folding measures for evaluation of developing human brain[END_REF] have recently shown, the curvature and its derivatives (e.g. the mean curvature, Gaussian curvature, shape index etc.) by themselves may not be able to discriminate between normal and aberrant cortical development. Nevertheless, the curvature contains useful information about the surface bending. Luders et al. [START_REF] Luders | A curvature-based approach to estimate local gyrification on the cortical surface[END_REF] defined a local GI by smoothing the magnitude of the mean curvature. Kim et al. [START_REF] Kim | Development of cortical shape in the human brain from 6 to 24 months of age via a novel measure of shape complexity[END_REF] has recently proposed a new GI by quantification of the shape index variance in a local region. Due to the heat kernel smoothing procedure used in these methods, they are likely to miss some folding features at fine scales.

Recently, the spectral analysis of mean curvature based on the spectrum of Laplace-Beltrami operator has been used to study cortical folding. By binning the global frequency distribution in several frequency bands and computing the dominant and the determinant bands, Germanaud et al. [START_REF] Germanaud | Larger is twistier: spectral analysis of gyrification (SPANGY) applied to adult brain size polymorphism[END_REF] derived two metrics of folding that can discriminate primary, secondary and tertiary folds. In a work by Shishegar et al. [START_REF] Shishegar | Quantifying gyrification using Laplace Beltrami eigenfunction level-sets[END_REF], a new GI has been proposed by computing the weighted differential mean curvature along the level sets of the first nontrivial Laplace-Beltrami eigenfunction. This method focuses mainly on sulcal bending and assigns the lowest values to points on gyri.

The above-mentioned methods have different implicit interpretations of the concept of "surface complexity" that underlies GIs and sometimes their findings lack consistency [START_REF] Awate | Cerebral cortical folding analysis with multivariate modeling and testing: Studies on gender differences and neonatal development[END_REF]. In this paper, we propose an explicit interpretation of this notion in an intuitive way that relies on the surface bending properties :

In a neighbourhood around each point of the cerebral cortex, the surface complexity is quantified by I. the magnitude of sulcal/gyral bending, or II. the spatial frequency of sulcal/gyral bends.

Due to the natural link between the surface bending and the mean curvature, the above interpretations mean that in a highly folded region, (I) the magnitude or (II) the variations of the mean curvature are higher in comparison to other regions.

Given those interpretations, we propose two novel local GIs, called spectral gyrification index (sGI) and weighted spectral gyrification index (wGI) based on a local spectral analysis of the mean curvature. These GIs are scale invariant and are computed directly on the cortical surface with neither of them requiring a reference surface nor a smoothing procedure. While sGI gives information about the magnitude of the mean curvature in a neighbourhood around each point on the cortical surface (Interpretation I), wGI takes into account the spatial frequency of folds in that neighbourhood (Interpretation II). The intrinsic nature of the method enables us to compute the surface complexity at different spatial scales. Moreover, in contrast to the surface area-based GIs, the proposed indices are able to disentangle the effect of depth on folding quantification. Finally, the issue of inconsistent analysis arising from the inter-subject brain size variability is also addressed by introducing an adaptive neighbourhood. To our best knowledge, the GIs proposed by Li et al. [START_REF] Li | Mapping longitudinal development of local cortical gyrification in infants from birth to 2 years of age[END_REF] and Kim et al. [START_REF] Kim | Development of cortical shape in the human brain from 6 to 24 months of age via a novel measure of shape complexity[END_REF] are the only ones that tried to tackle this issue before.

In the following, we present in Section II a local spectral analysis method based on the graph windowed Fourier transform of the mean curvature. This analysis is used to define the two gyrification indices, sGI and wGI, together with their properties. The method is applied to some synthetic surfaces and real data in Section III. Finally, we discuss the results and the specificities of our method, followed by a conclusion.

II. METHOD

In this section, we use a local spectral analysis of the cortical surface mesh to define the mesh windowed Fourier transform of the mean curvature function and we show that it can be used to define two new gyrification indices. Spectral methods for surface analysis rely mainly on eigenvalues and/or eigenfunctions of an operator defined on the surface. For example, for a Riemannian surface, the eigenfunctions of the Laplace-Beltrami operator serve as bases for Fourier transforms. The idea has been extended to triangular meshes modelling surfaces to analyse their structural properties; see [START_REF] Zhang | Spectral methods for mesh processing and analysis[END_REF] for a comprehensive survey on this topic. Here, we use such eigenfunctions to define a mesh windowed Fourier transform and provide an explicit access to the local frequency components of the mean curvature of the cortical surface.

A. Mesh Fourier Transform

For a compact Riemannian manifold S as a surface in IR 3 , one can consider a set of square integrable functions defined on the surface:

L 2 (S) = {u ∶ S → IR| ∫ S u 2 < ∞}.
The Laplace-Beltrami operator ∆, associated with the surface S, is defined as a generalization of the Laplacian operator in Euclidean space to Riemannian manifolds. The spectrum of this operator {(λ i , u i ) ∈ IR + × L 2 (S), i = 1, 2, . . .} is generated by solving the differential eigenvalue problem of ∆: ∆u i = -λ i u i in which λ i and u i are called the ith-eigenvalue and eigenfunction of ∆ [START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF]. The spectral theory based on the Laplace-Beltrami spectrum can be used to obtain a new representation of the space L 2 in the so called spectral domain. In practice, the differential eigenvalue problem of ∆ is usually solved on a triangulated surface by numerical methods.

Formally, let G = {V, E} be a triangular mesh modelling the surface S where V is the set of vertices, V = {P 1 , P 2 , . . . , P N }, |V | = N < ∞ and E is the set of edges. By using the linear finite element method (FEM) [START_REF] Reuter | Discrete Laplace-Beltrami operators for shape analysis and segmentation[END_REF], the abovementioned differential eigenvalue problem is discretized to the following algebraic generalized eigenvalue problem

Aχ = λBχ, ( 1
)
where A and B are N × N sparse matrices with the following elements:

A(i, j) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ cot αij +cot βij 2 if (i, j) ∈ E, -∑ k∈N (i) A(i, k) if i = j, 0 o.w. (2)
and

B(i, j) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ |t1|+|t2| 2 if (i, j) ∈ E, ∑ k∈N (i) |t k | 6 if i = j, 0 o.w. ( 3
)
where α ij and β ij are the angles opposite to the edge P i P j in two triangles t The bases of the Fourier transform in continuous domain are complex exponential functions which are the eigenfunctions of the Laplace-Beltrami operator. Inspired by this fact, the eigenvectors of the discretized Laplace-Beltrami operator serve as Fourier atoms on the triangulation [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF]. Given a function f defined on the vertices of triangulation, Fourier transform coefficients of f are given by the set { f (l) ∶= ⟨f, χ l ⟩ B , l = 1, 2, . . . , N }. This set is a representation of function f in the spectral domain and gives the frequency distribution of this function. In this setting, the Parseval's identity is

⟨f, g⟩ B = ⟨ f , ĝ⟩, (5)
where ⟨., .⟩ denotes the Euclidean inner product. It yields ∥f ∥ B = ∥ f ∥ 2 where ∥.∥ B is the norm induced by the B-inner product.

B. Local Mesh Fourier Transform

The mesh Fourier transform described in Section II-A gives a global frequency distribution of a function f . In other words, it is not able to provide information about the frequencies of f in a spatial neighbourhood. In continuous domain, the windowed Fourier transform was introduced in order to find the local frequency distribution of a function [START_REF] Gabor | Theory of communication[END_REF]. This transform gives local information about a function simultaneously in spatial and frequency domains.

Shuman et al. [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF] have recently extended this transform to graph setting which enables us to do a "vertex-frequency analysis". The general idea of this method is to localize a function defined on the vertices of a graph around a vertex by a translated window function and then, compute the graph Fourier transform of this localized function. In this paper, we extend the method to the mesh framework by using the mesh Fourier transform that is aware of the surface geometry.

Since a mesh is a special kind of a graph, embedded in a surface, one may argue that the graph spectral theory tools can be applied on a mesh without any adaptation. In general, it is true but unlike the graph Laplacian in which only the connectivity of vertices are considered, the geometric (FEM) Laplacian takes into account geometric properties of the surface. Equations ( 2) and (3) show how the geometry of neighbouring triangles on the mesh contributes to the definition of Laplacian operator. More discussions and comparisons between the graph Laplacian and the geometric Laplacian can be found in [START_REF] Zhang | Spectral methods for mesh processing and analysis[END_REF], [START_REF] Reuter | Discrete Laplace-Beltrami operators for shape analysis and segmentation[END_REF], [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF], [START_REF] Levy | Laplace-Beltrami Eigenfunctions Towards an Algorithm That "Understands" Geometry[END_REF]- [START_REF] Tan | Spectral Laplace-Beltrami wavelets with applications in medical images[END_REF].

1) Window function and translation operator

Let f ∶ V → IR be a function defined on the vertices of a triangulation (e.g. the mean curvature). To localize this function around a specific vertex, we need a window function with local support and a translation operator to move the window function to that specific vertex. Following [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF], we consider the window function

ĝ(l) = C exp(-τ λ l ), (6)
defined in the spectral domain. In this formula, τ is a parameter which determines the size of the window, λ l is the l-th Laplace-Beltrami eigenvalue and C is chosen such that ∥ĝ∥ 2 = 1.

The window size parameter τ sets a locality trade-off between the frequency and spatial domains [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF], [START_REF] Kaiser | Windowed Fourier Transforms[END_REF]. By increasing the window size, we will have a wider window in the spatial domain and the function f will be localized in a larger neighbourhood around each vertex. On the other hand, more local frequency distribution of the function in that neighbourhood is obtained. The spread of window function in spatial and frequency domains is measured by the area of the Heisenberg box [31, Section 4.2], [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF]Section 6.6]. It is proved that the Gaussian function is the unique window that minimizes the area of the Heisenberg box [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]Theorem 2.5]. Since λ l is proportional to the square of the spatial frequency [START_REF] Grebenkov | Geometrical structure of Laplacian eigenfunctions[END_REF], [START_REF] Golbabai | Hybrid shape parameter strategy for the RBF approximation of vibrating systems[END_REF], i.e. λ l ∝ ω 2 l , the window function ( 6) corresponds to a Gaussian function in the frequency domain ĝ(l) ∝ exp(-τ ω 2 l ). The fact that there is no canonical origin and direction on a triangulated mesh makes it difficult to define a translation operator on this mesh. Inspired by the properties of the generalized Fourier transform in continuous domain and the mechanism proposed in [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF] for graph setting, the translation operator is defined on a triangulation as (see Appendix A for the details of derivation):

T i ∶ IR N → IR N i = 1, 2, . . . , N, (T i g)(n) ∶= √ N N ∑ l=1 N ∑ m=1 (ĝ(l)B(i, m)χ l (m)χ l (n)). (7)
The translation operator T i shifts the center of the window function to a vertex P i . Now, by multiplying the function f by the translated window function T i g, it is localized around the vertex P i :

fi (n) = (T i g)(n)f (n), n = 1, 2, . . . , N. (8)

2) Mesh windowed Fourier transform

The mesh windowed Fourier transform coefficients of a function f ∈ IR N are defined as the modulation of the localized function fi by Fourier atoms {χ k , k = 1, 2, . . . , N }:

Sf (i, k) ∶= ⟨ fi , χ k ⟩ B , ( 9
)
where i = 1, 2, . . . , N is the index of vertex and k is the index of frequency. This gives us a frequency spectrum {|Sf (i, k)| 2 , k = 1, 2, . . . , N } for every vertex P i of the mesh which can be seen as the frequency distribution of the function f in a local neighbourhood around the vertex P i .

3) Adaptive window function

It is easily seen from Eqs. ( 1)-( 3) that if the surface area is scaled by a factor q 2 , the eigenvalues are scaled by 1/q 2 . In this case, due to the definition of the window function [START_REF] Schaer | A surface-based approach to quantify local cortical gyrification[END_REF], the relative spread of window (i.e. the ratio between the region area covered by the window and the surface area) is affected by the size of the surface. In other words, given a fixed window size τ for the original and the scaled surfaces, the window covers relatively larger space on the smaller surface and vice versa. This leads to an inconsistent large and small scale spectral analysis for small and large surfaces, respectively (see Fig. S1a of supplementary figures).

To keep the relative spread of window constant automatically across surfaces (Fig. S1b), we introduce an adaptive window function in which the surface area is incorporated:

ĝ(l) = C exp(-τ |S|λ l ), (10)
where |S| denotes the surface area of surface S. Note that by this definition, a dimensionless parameter is obtained inside the exponential. As we will see in Section II-C1, the adaptive window function also plays an important role to derive scale invariant gyrification indices. In Fig. 1, the spread of the adaptive window function (10) with 3 different window size parameters, τ = 2e -4, 1e -3 and 5e -3 is plotted on a cortical surface around a specific vertex shown as a yellow point on the precentral gyrus. As τ increases, the spread of window in spatial domain increases as well. While the narrow window, τ = 2e -4, covers a part of a gyrus and/or sulcus, the medium window, τ = 1e -3, covers several folds and the wide window, τ = 5e -3, covers a big portion of the cortical surface equivalent to a lobe. Tuning the window size parameter τ changes the spatial scale of the analysis.

C. Gyrification Indices

To define gyrification indices that fulfill our interpretations on the surface complexity, stated in Introduction section, we take the mean curvature as a function defined on the vertices of a mesh. By applying the mesh windowed Fourier transform to this function, we will have a frequency spectrum at each vertex P i of the mesh which consists of the so called "frequency powers" |Sf (i, k)| 2 , k = 1, 2, . . . , N . The summation of the frequency powers is called the total power (TP) of the frequency spectrum [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. We proved that the TP of the frequency spectrum of the vertex P i equals to the norm of the localized mean curvature around this vertex (see Proposition 1 in Appendix B):

N ∑ k=1 |Sf (i, k)| 2 = ∥ fi ∥ 2 B . ( 11
)
Due to our first interpretation, in a more folded region, the magnitude of the mean curvature increases. This interpretation, together with Eq. ( 11), leads to the definition of the first GI:

I. Spectral Gyrification Index (sGI):

sGI(i, S) = N ∑ k=1 |Sf (i, k)| 2 . ( 12
)
In this equation, sGI(i, S) denotes the spectral gyrification index of vertex P i of surface S.

On the other hand, based on our second interpretation, the variation of the mean curvature increases in more folded regions. Since the high variation of a function is encoded in the high frequency band of its frequency spectrum, we give larger weights to higher frequency powers to take into account this interpretation. It brings us the second GI: II. Weighted Spectral Gyrification Index (wGI):

wGI(i, S) = N ∑ k=1 ( λ k λ 2 ) 2 |Sf (i, k)| 2 . ( 13
)
In this equation, the weights are the normalized eigenvalues of the Laplace-Beltrami operator that contain information about the shape of surface [START_REF] Reuter | LaplaceBeltrami spectra as Shape-DNA of surfaces and solids[END_REF], [START_REF] Wachinger | BrainPrint: A discriminative characterization of brain morphology[END_REF]. The normalization by the first nonzero eigenvalue λ 2 removes the effect of the size of surface on weighting [START_REF] Lefèvre | Fast surface-based measurements using first eigenfunction of the laplace-beltrami operator: Interest for sulcal description[END_REF]. In this definition, both Laplace-Beltrami eigenvalues and eigenvectors are involved. We proved that

wGI(i, S) = 1 λ 2 2 ∥L fi ∥ 2 B , ( 14
)
where L = B -1 A is the Laplace-Beltrami operator (see Proposition 2 in Appendix B). Since this operator measures how much a function differs at a point from its average value at neighbour points, (L fi )(m) measures the variation of fi at vertex P m and so, wGI(i, S) sums up all of these variations of localized f around vertex P i .

1) Geometric invariance

We now provide some important properties of the proposed GIs. The Laplace-Beltrami spectrum is invariant under isometric transformations. It makes sGI and wGI isometry invariant. Moreover, as demonstrated in Proposition 3 of Appendix B, sGI and wGI are both scale invariant by their constructions.

2) Global GIs Given a surface mesh with M triangular faces S = {t 1 , t 2 , . . . , t M }, GI defines a function on the vertices of this mesh. We define the global GI of surface S as the mean value of this function:

GI(S) = 1 |S| M ∑ j=1 1 3 3 ∑ i=1 GI j (i, S)|t j |, (15)
where GI j (i, S) denotes the GI value (sGI or wGI) of the i-th vertex of the face t j and |t j | is the area of t j . It is noteworthy that the global GI is scale invariant as the local ones.

III. EXPERIMENTS AND RESULTS

A. Synthetic data

In order to illustrate the efficiency of sGI and wGI for measuring the degree of folding of a surface, and to show the effect of sulcal depth on surface area-based methods, we computed our GIs and Toro's GI [START_REF] Toro | Brain size and folding of the human cerebral cortex[END_REF] on some synthetic surfaces for which we control the degree of folding.

1) Wavy rectangle A wavy rectangle in IR 3 is created by the equation z = 2 sin(60πx 2 )/(60πx) where -0.7 ≤ x ≤ 0.7 and 0 ≤ y ≤ 1 with geodesic length equals to 4. This surface is triangulated with N = 40, 000 equidistant vertices and is shown in Fig. 2a. The intersection of this surface with the plane y = 0.5 is indicated on the surface by green points and is plotted in Fig. 2b. As we go far from the center to the left or right side of this surface, it becomes more folded i.e. the surface becomes more bended and the spatial frequency of the folds increases. The windowed Fourier transform is applied on the mean curvature of this surface and the spectrogram that consists of the frequency powers |Sf (i, k)| 2 is shown in supplemental Fig. S2a (see supplementary figures).

In Fig. 2b, the neighbourhoods of two vertices P m and P n are shown. Obviously, the surface is more folded around P m than P n while there are deeper folds around the latter. The local frequency spectrums of these vertices are plotted in the supplemental Figs. S2b and S2c that show clearly how higher frequency powers increase in regions with high oscillating folds. The values of sGI, wGI and a surface area-based GI defined by Toro et al. [START_REF] Toro | Brain size and folding of the human cerebral cortex[END_REF] for those vertices are given in Table I for comparison. While sGI and wGI give appropriately higher values to vertex P m , Toro's GI gives equal values to both vertices. It shows that surface area-based GIs may not be able to discriminate between deep folds and complex folds with the same area. sGI, wGI and Toro's GI values of the vertices along the green middle line of the wavy surface are plotted in Fig. 2c. The observations clarify that sGI and wGI discriminate efficiently deep folds from complex ones with a clear increase from the center towards the borders whereas Toro's GI shows high values on the center where there are deep folds. To elucidate more the effect of the fold depth on the surface area-based GIs, two other surfaces are constructed; see Figs. S3 and S4 of the supplementary figures. These surfaces are designed to study the effect of fold depth and oscillation frequency separately.

B. Real Data

In this section, the proposed method is applied to a real database of cortical surfaces. The gyrification maps of an individual subject as well as group average maps across all subjects are presented. The relation between the proposed GIs and the brain volume is also studied.

1) Data and preprocessing

We applied the method to 124 healthy adult subjects from the Open Access Series of Imaging Studies (OASIS) database 1 . For each subject, three or four T1 anatomical Magnetic Resonance Images (MRIs) had been acquired at inplane resolution of 1mm × 1mm, slice thickness = 1.25 mm, TR = 9.7 ms, TE = 4 ms, flip angle = 10u, TI = 20 ms, TD = 200 ms. Images of each subject were motion corrected and averaged to create a single image per subject with a high contrast-to-noise ratio. The resulting anatomical MR images were segmented using BrainVISA 2 . The white matter surface of each hemisphere was then meshed using this software which resulted in triangular meshes with spherical topology and approximately 50,000 nodes depending on the subjects. The Hip-Hop algorithm [START_REF] Auzias | Model-driven harmonic parameterization of the cortical surface: HIP-HOP[END_REF] in BrainVISA was then applied to compute the spherical interindividual correspondence between cortical surfaces.

2

) Gyrification maps

The maps of sGI and wGI of a subject from the database are shown in Fig. 3 for 3 different window sizes, τ = 2e -4, 1e -3, 5e -3 associated to very local, medium and wide windows, respectively (see Fig. 1 for window sizes). On this figure, it is shown that the window size parameter τ can be used to control the scale of observations. At τ = 2e -4, the spatial scale is fine and high values are located mostly on the ridge of complex gyri, while low values are located on the walls of regular sulci. As the window size increases, a more regional effect becomes visible, with a very smooth and low variation map at value τ = 5e -3, which gives a coarse scale observation of the gyrification.

As it has been discussed in Section II-C, the indices sGI and wGI, based on their constructions, measure complementary properties of surface folding. This is shown in Fig. 4 where the mean curvature of the same cortical surface of Fig. 3 is depicted. Two regions on this surface, R1 and R2, have been 2 http://brainvisa.info/ chosen. R1 is a sharp spike located on the postcentral gyrus and R2 is a very shallow fold located on the superior parietal lobe. The mean curvature of the region R1 is very high while that of the region R2 varies a lot between positive and negative values. The maps of sGI and wGI of R1 and R2 are shown in this figure (τ = 2e -4). As expected theoretically from Eqs. ( 11) and ( 14), sGI gives high value to R1 while wGI assigns high value to R2.

The group average of both GIs in the medium scale (window size τ = 1e -3) have been computed by using the cortical surface inter-subject matching method, Hip-Hop [START_REF] Auzias | Model-driven harmonic parameterization of the cortical surface: HIP-HOP[END_REF]. Results are depicted on the template cortical surface hiphop1383 in Fig. 5 for the left hemisphere. The average patterns of sGI and wGI represented in this figure are similar to those observable on an individual subject (Fig. 3) which shows that the spatial patterns of the proposed GIs are reproducible across subjects. As Fig. 5a shows, sGI gives higher values to vertices on the prefrontal and occipital lobes, inferior parietal lobe, inferior temporal sulcus and the medial area of the superior parietal cortex. wGI, as shown in Fig. 5b, assigns higher values to the prefrontal lobe, medial part of the occipital lobe and the posterior cingulate gyrus. Some folding is also captured by high wGI values in the insula. We also computed the average GIs of the right hemispheres and by visual inspection, we observed no remarkable difference in gyrification patterns of the left and right hemispheres in the medium scale (see Supplemental Fig. S5).

To compare with our results, the average map of Toro's GI across all subjects in the database is presented in Fig 5c .  This method gives high GI values to deep folds like the central sulcus, the insula, the superior temporal sulcus and the parietooccipital sulcus.

3) Scaling analysis

Recent studies show that larger brains are more folded [START_REF] Toro | Brain size and folding of the human cerebral cortex[END_REF], [START_REF] Germanaud | Simplified gyral pattern in severe developmental microcephalies? New insights from allometric modeling for spatial and spectral analysis of gyrification[END_REF], [START_REF] Germanaud | Larger is twistier: spectral analysis of gyrification (SPANGY) applied to adult brain size polymorphism[END_REF], [START_REF] Im | Brain size and cortical structure in the adult human brain[END_REF]. To investigate this phenomenon, the global sGI and wGI of each hemisphere are modelled by the following power law

G = kV α , ( 16
)
where G denotes the global gyrification index (sGI or wGI) of a hemisphere computed through Eq. ( 15), V is the hemispheric volume and k and α are coefficients to be determined. Since gyrification indices, sGI and wGI, are scale invariant (Proposition 3), the scaling exponent coefficient α of the power law [START_REF] Kim | Development of cortical shape in the human brain from 6 to 24 months of age via a novel measure of shape complexity[END_REF], under the isometric scaling of brain volume, should equal 0. The hemispheric data along with the fitted power law model in logarithmic scale are represented in Fig. 6. The positive exponent coefficient reveals a positive allometric scaling of gyrification indices with volume. Moreover, it confirms that the larger brains are more folded. The values of the exponent coefficient α of wGI is more than that of sGI (0.67 > 0.37 and 0.66 > 0.36 for left and right hemispheres, respectively) while the proportion of the variance of sGI explained by the volume is higher than that of wGI (0.44 > 0.36 and 0.47 > 0.34 for left and right hemispheres, respectively). The results also suggest that in the global hemispheric scale, the degree of folding of the left and right hemispheres increase with volume symmetrically.

The above-mentioned results support this hypothesis that the larger brains are more folded but they do not illustrate which cortical regions get more folded in larger brains. To address this question, we perform the allometric analysis ( 16) at the vertex level. Thanks to the inter-subject matching by Hip-Hop method [START_REF] Auzias | Model-driven harmonic parameterization of the cortical surface: HIP-HOP[END_REF], we are able to find the corresponding vertices across all subjects. Our method naturally allows us to monitor the changes of the degree of gyrification at different spatial scales. The maps of the significant exponent coefficient α (p < 0.05, corrected for multiple comparisons using False Discovery Rate (FDR) [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]) for sGI and wGI of the left hemispheres derived from the medium window (τ = 1e -3) are represented in Fig. 7. In this figure, the vertices for which α is not significant are masked by the gray color. Fig. 7 shows that in the adult population, as brain size increases across subjects, deep folds with low gyrification (represented in blue in Fig. 5a and5b), such as the central sulcus, the insula, the superior temporal sulcus, or the parieto-occipital sulcus, show the largest increase in folding complexity. We have not observed any vertex with significant negative exponent except in the right hemisphere with the most local window (τ = 2e -4). In this case, for sGI map, we observed about 0.04% of vertices with significant negative α exponent located on the anterior cingulate and the superior parietal cortices. For wGI, there are only 0.02% of vertices with significant negative α exponent located on the isthmus cingulate cortex. The αmaps of the right hemisphere (τ = 1e -3) are given in Fig. S6 (see supplementary figures).

IV. DISCUSSION AND CONCLUSION

In this study, we interpreted the cortical surface complexity in two ways that are consistent with the intuitive concept of it. The regions where α is not significant (n.s.) are masked by gray color. The higher the value of α, the higher the variation of folding with respect to volume.

Based on these assumptions and by applying the local spectral analysis of the mean curvature, we developed a method to estimate locally the degree of folding of the cerebral cortex. The method was applied to a database of 124 healthy adult brains and the results demonstrate that the primary folds like the central sulcus and the insula are less folded than other regions. The central sulcus is a good example since it is a very deep fold but with a high regularity and straight shape where we expect a low measure of folding whereas surface area-based methods usually show high gyrification values (see Fig. 5).

Through some experiments on synthetic surfaces, we clarified that surface area-based methods such as Toro's [5] and Schaer's [START_REF] Schaer | A surface-based approach to quantify local cortical gyrification[END_REF] GIs may fail to discriminate between deep and complex folds. Indeed, these GIs give higher values to deep folds; see the similar gyrification maps of Toro's GI [5, Fig. 5a-c] and Schaer's GI [START_REF] Schaer | A surface-based approach to quantify local cortical gyrification[END_REF]Fig. 4]. To address this issue, Su et al. [START_REF] Su | Geometric computation of human gyrification indexes from magnetic resonance images[END_REF] proposed to weight the GI by the geodesic sulcal depth. Nevertheless, we believe that this weighting strategy is in contradiction with the notion of folding. In fact, if a cortical region with deep sulci and another region with shallow sulci have equal areas, the latter region should be more folded to keep the same area as the former one.

As an illustrative example to show how our method discriminates between deep and oscillating folds of a cortical surface, the medial face of a subject is considered and its sGI and wGI maps are shown in Fig. 8a and 8b respectively (the lateral maps are shown in the middle row of Fig. 3). Two lines with equivalent geodesic length were drawn on the medial side: line 1 in the medial prefrontal region (red line, geodesic length 35.35 mm) where both GIs show high values, and line 2 in the medial precentral region (blue line, geodesic length 34.63 mm) where both GIs have low values (Fig. 8c). Geodesic sulcal depth on the cortical surface was computed [START_REF] Im | Brain size and cortical structure in the adult human brain[END_REF] in order to get depth values at each point of both lines and produce depth curves along these lines. These depth curves, plotted in Fig. 8d and8e, show that the frequency of folds in the medial prefrontal region (Line 1) is double of that in the medial precentral region (Line 2), while the depth of folds on the medial prefrontal region (Line 1) is almost half of that in the medial precentral region. This explains the high values of our gyrification index in the medial prefrontal region: despite an apparent smoothness due to the low sulcal depth, the folding frequency is high. Another difference between our method and surface areabased methods is that the latter may be not localized enough for some applications. For example, in Fig. 4 of [START_REF] Schaer | A surface-based approach to quantify local cortical gyrification[END_REF], for a small spherical neighbourhood, the most folded region of the cortex is around the Sylvian Fissure and as the size of the neighbourhood increases, the same pattern propagates across the cortex. Therefore, it may fail to catch other folded parts of the brain, thus affecting the reliability of findings. In our method, by tuning the neighbourhood size, we get the results at different spatial scales, ranging from a very local scale (in order of a part of a sulcus/gyrus) to a more global scale (in order of a lobar cortex); see Figs. 1 and3.

Despite the differences we highlighted between our method and previous GI definitions, we do think that these previous gyrification indices are still valid but should be considered in the light of their specific interpretations of the notion of "surface complexity". Our experiments on synthetic surfaces show that how different interpretations of this notion may lead to different results. Accordingly, a message of this paper is to point out that one should pay attention to this notion and interpret results based on it. Another work on this topic has been done by Awate et al. [START_REF] Awate | Cerebral cortical folding analysis with multivariate modeling and testing: Studies on gender differences and neonatal development[END_REF] discussing about the different meanings of surface complexity and categorizing GIs based on their responses to different situations like surface scaling and increased spatial frequency of folds.

An important matter that has been mostly neglected so far is the effect of brain size on the relative neighbourhood spread. Due to inter-subject size variability, analysis of the cerebral cortices with a fixed neighbourhood size may cause inconsistency. To address this issue, we introduced a neighbourhood which is adapted for each hemisphere by its surface area. In longitudinal/developmental studies, however, a more local adaptation is needed to deal with the non-uniform changes of local surface area. To quantify the changes of folding complexity of infants within 6 months to 24 months of age, Kim et al. [START_REF] Kim | Development of cortical shape in the human brain from 6 to 24 months of age via a novel measure of shape complexity[END_REF] introduced a two-phase adaptation of neighbourhood. First, for each 6-month-old subject, a fixed neighbourhood size is adapted globally by the ratio of the subject's hemisphere surface area to the average hemisphere surface area. Then, for older subjects, the normalized neighbourhood size is scaled by the changes of local surface area. In a developmental study of local cortical gyrification in infants from birth to 2 years of age, Li et al. [START_REF] Li | Mapping longitudinal development of local cortical gyrification in infants from birth to 2 years of age[END_REF] used the so-called N -ring neighbourhood around each mesh vertex. By resampling all cortical meshes to have equal number of vertices, they took into account the local intra-and inter-subject brain size variability. Although out of the scope of this paper, a similar approach could be employed in our method.

To investigate the relationship between the volume and the degree of folding of the cerebral cortex, we considered a power law to model the GI as a response to the volume. The resulted positive exponent coefficient α indicates a positive allometric relation. It implies that if the brain volume is doubled, the GI is multiplied by factor 2 α . Moreover, it supports the well-known hypothesis in literature that larger brains are more folded. This observation is in accordance with some mechanical models of cortical folding process [START_REF] Tallinen | Surface sulci in squeezed soft solids[END_REF], [START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF] and deserves further investigations, in particular in longitudinal databases. The relative low coefficient of determination, R 2 , clarifies that there is still enough room for other covariates, beyond the volume, to explain the degree of gyrification. One interesting direction for future studies is to take into account some biological factors like age, sex, genetic conditions etc. or cognitive and behavioural ones like special skills, IQ etc. in the scaling model.

We also fitted the power law model to GIs at the vertex level to investigate how the degree of folding changes locally with the hemispheric volume. The results illustrate that the less folded cortical regions, in terms of either the magnitude of the mean curvature or its variation, like the walls of the precentral and postcentral gyri and the insula, are more convoluted in larger brains. Regarding the allometric relation between the brain volume and the cortical surface area [START_REF] Toro | Brain size and folding of the human cerebral cortex[END_REF], one may conclude that the cerebral cortex of a larger brain is twisted in less folded regions to accommodate the additional surface.

We believe that the two new gyrification indices we proposed, sGI and wGI, will be useful to investigate various aspects of cortical folding complexity and variability, such as the link between brain size and gyrification but also the search for biomarkers of developmental pathologies, or the understanding of cortical development.

APPENDIX A TRANSLATION OPERATOR

The translation operator in the mesh setting is inspired by that in continuous and graph settings [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF]. More precisely, let {ψ l , l = 1, 2, . . .} be the basis of the generalized Fourier transform in continuous domain and assume that those functions are orthonormal with respect to a function ϕ i.e. ∫ ψ l ϕψ k = δ kl . From the properties of the generalized Fourier transform, translation in spatial domain causes modulation in Fourier domain: h(x) = f (x -x 0 ) ⇔ ĥ(k) = ψ k (x 0 )ϕ(x 0 ) f (k). [START_REF] Shishegar | Quantifying gyrification using Laplace Beltrami eigenfunction level-sets[END_REF] In other words, the translation of function f to point x 0 can be given by the inverse Fourier transform of the modulated Fourier coefficients f :

(T x0 f )(x) ∶= f (x -x 0 ) = F -1 {ψ k (x 0 )ϕ(x 0 ) f (k)}, (18)
where F -1 denotes the inverse Fourier transform. In mesh setting, the Fourier basis is the set of Laplace-Beltrami eigenvectors which are orthonormal with respect to the matrix B. If g is a function defined on vertices of mesh, inspired by Eq. ( 18), the translation of g to vertex P i is defined as following

(T i g)(n) (18) ∶= √ N F -1 {B(i, ∶)χ l ĝ(l)} = √ N N ∑ l=1 N ∑ m=1 χ l (n)(B(i, m)χ l (m)ĝ(l)),
where B(i, ∶) denotes the i-th row of B. Another derivation of the translation operator through convolution in graph setting can be found in [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF].

APPENDIX B THEORETICAL COMPUTATIONS OF SGI AND WGI

Proposition 1. For gyrification index sGI at vertex P i of subject S we have sGI(i, S) = ∥ fi ∥ 2 B .

Fig. 1 .

 1 Fig. 1. The spread of the adaptive window function (10) with 3 different window sizes around the yellow vertex on a cortical surface by thresholding the values of the window functions at the same percentage of their maximal values. The red color highlights the window spread around the vertex. (a) A narrow window, with τ = 2e -4, covers about 1.5% of the surface. (b) A medium window, with τ = 1e -3, covers about 8% of the surface. (c) A wide window, with τ = 5e -3, covers about 36% of the surface.

Fig. 2 .

 2 Fig. 2. (a) Wavy rectangle in IR 3 . The "middle line" on the surface is formed by the vertices lie on the surface with the the same y-coordinate 0.5. (b) The middle line along with zoom on two neighbourhoods around vertices Pm and Pn. (c) Values of Toro's GI, sGI and wGI of the vertices located on the middle line.

Fig. 3 .

 3 Fig. 3. The map of gyrification indices (a) sGI and (b) wGI of the left hemisphere of an individual subject from our database at 3 different scales. In the first row τ = 2e -4, in the second row τ = 1e -3 and in the third row τ = 5e -3. The blue and red colors indicate the extremes of low and high degree of folding respectively.

Fig. 4 .

 4 Fig. 4. Zoom on 2 regions, R1 and R2, of the cerebral cortex. The colormap of the cortex encodes its mean curvature (the blue and red colors indicate the extremes of negative and positive values respectively). sGI and wGI (τ = 2e -4) of the regions R1 and R2 are also represented (the blue and red colors indicate the extremes of low and high degree of folding respectively).

Fig. 5 .

 5 Fig. 5. Average of gyrification indices: (a) sGI and (b) wGI derived from the medium window (τ = 1e-3), and (c) Toro's GI with spherical neighbourhood radius r = 20 across the left hemispheres of all 124 subjects on the hiphop138 template surface 3 . The blue and red colors indicate the extremes of low and high degree of folding respectively.

Fig. 6 .

 6 Fig. 6. Up: Relationship between the volume (mL) of the left/right hemispheres and global wGI (τ = 1e -3). The fitted line for the left hemisphere is y = 0.67x + 14.57, R 2 = 0.36, p < 0.001. The fitted line for the right hemisphere is y = 0.66x + 14.74, R 2 = 0.34, p < 0.001. Down: Relationship between the volume (mL) of the left/right hemispheres and global sGI (τ = 1e-3). The fitted line for the left hemisphere is y = 0.37x+3.03, R 2 = 0.44, p < 0.001. The fitted line for the right hemisphere is y = 0.36x + 3.2, R 2 = 0.47, p < 0.001.

Fig. 7 .

 7 Fig. 7. The colormaps encode the significant exponent α (p < 0.05, corrected for multiple comparisons using FDR) of vertex-wise allometric analysis (16) of gyrification indices (a) sGI and (b) wGI of the left hemisphere derived from the medium window (τ = 1e -3). The magenta and red colors indicate the extremes of low and high values of α.The regions where α is not significant (n.s.) are masked by gray color. The higher the value of α, the higher the variation of folding with respect to volume.

2 Fig. 8 .

 28 Fig. 8. Medical view of the map of (a) sGI and (b) wGI of a subject. The blue and red colors indicate the extremes of low and high degree of folding respectively. (c) Two lines with almost equal geodesic length are drawn on the medial prefrontal region (Line 1, red, geodesic length=35.35 mm) and the medial central region (Line 2, blue, geodesic length=34.63 mm). The depth maps of (d) Line 1 and (e) Line 2 show that the frequency of folds in the medial prefrontal region is double of that in the medical precentral region while the depth of folds on the medial prefrontal region is almost half of that in the medial precentral region.

  [START_REF] Armstrong | The ontogeny of human gyrification[END_REF] and t 2 sharing this edge, |t k | indicates the area of the triangle t k and N (i) denotes the index set of all vertices of the 1-ring neighbourhood of P i . The matrix B is positive definite and defines the so called B-inner product in IR N : ∀f, g ∈ IR N , ⟨f, g⟩ B = f t Bg. = λ 1 < λ 2 ≤ . . . ≤ λ N and a set of eigenvectors {χ j , j = 1, 2, . . . , N } in IR N which are orthonormal with respect to B-inner product i.e. ⟨χ i , χ j ⟩ B = δ ij where δ ij is Kronecker delta.

	(4)
	Solutions of the discrete eigenvalue problem (1) are nonneg-
	ative real eigenvalues 0

TABLE I DIFFERENT

 I GIS FOR VERTICES Pm AND Pn OF THE WAVY RECTANGLE DEPICTED IN FIG. 2B

	vertex	sGI	wGI	Toro's GI
	Pm	8.5561e2 2.8090e10	2.2268
	Pn	1.8487e2	3.8802e9	2.2271

http://www.oasis-brains.org/

The hiphop138 template is available at http://www.meca-brain.org/ softwares/hiphop138-cortical-surface-group-template/

 [START_REF] Lefèvre | Are developmental trajectories of cortical folding comparable between crosssectional datasets of fetuses and preterm newborns?[END_REF]http://labex-archimede.univ-amu.fr/ and sharing MRI databases of the brain. The OASIS project was supported by NIH grants P50 AG05681, P01 AG03991, R01 AG021910, P50 MH071616, U24 RR021382, and R01 MH56584. We also thank Antonietta Pepe and Guillaume Auzias for very helpful discussions. We are grateful to anonymous reviewers for their constructive suggestions.

Proof:

where [START_REF] Germanaud | Larger is twistier: spectral analysis of gyrification (SPANGY) applied to adult brain size polymorphism[END_REF] and [START_REF] Zhang | Spectral methods for mesh processing and analysis[END_REF] are derived by Parseval's identity ( 5) and ( 20) is based on the fact that χk = δ k (δ k is Kronecker delta).

Lemma 1. Let f ∈ IR N be a function defined on the vertices of a triangulation and L ∈ IR N ×N be the discrete Laplace-Beltrami operator i.e. L = B -1 A. Then, the Fourier coefficients of the function

Proof: The Fourier coefficients of y are as

where Eq. ( 22) is given by Eq. ( 1) and the symmetry of matrices A and B, and the definition of B-inner product (4) gives [START_REF] Reuter | Discrete Laplace-Beltrami operators for shape analysis and segmentation[END_REF].

Proposition 2. For gyrification index wGI at vertex P i of subject S we have wGI(i, S) = 1

B . Proof: Lemma 1 and the Parseval's identity give the first and the second equality, respectively:

Proposition 3. The gyrification indices sGI and wGI are scale invariant.

Proof: Assume that the surface S 2 is the scaled version of the surface S 1 by a factor q 2 i.e. |S 2 | = q 2 |S 1 |. Then the Laplace-Beltrami eigenvalues are scaled by 1/q 2 while the eigenvectors and the mean curvature are scaled by 1/q.

Thanks to the adaptive window function [START_REF] Shimony | Comparison of cortical folding measures for evaluation of developing human brain[END_REF], the translation operator and the mesh windowed Fourier coefficients Sf (i, k) remain unchanged because

It implies that sGI and wGI remain unchanged under scaling.

ACKNOWLEDGMENT

The work was supported by Labex Archimède 4 . The authors would like to thank the OASIS project for providing


Andrea Brovelli
email: andrea.brovelli@univ-amu.fr

Jean-Michel Badier

Francesca Bonini

Fabrice Bartolomei

Olivier Coulon

Guillaume Auzias

Sophie Chen

Patrick Marquis

Catherine Liegeois-Chauvel

Dynamic reconfiguration of visuomotor-related functional connectivity networks

This work was partly supported by Projets exploratoires pluridisciplinaires (PEPS) from the CNRS/Inserm/Inria "Bio-math-Info" entitled "Computational and neurophysiological bases of goal-directed and habit learning". This work was performed on a MEG platform member of France Life Imaging network (grant ANR-11-INBS-0006).