Quantitative methods

The normal curvature on a 3D surface in some direction is the inverse of the radius of the circle that best approximates a surface normal slice in that direction³⁵. The normal curvature for a smooth surface can be represented by the Weingarten matrix, i.e. the second fundamental tensor II, which is defined in terms of the directional derivatives of the surface normal:

where (u, v) are the directions of an orthonormal coordinate system in the tangent frame (the sign convention used here produces positive curvatures for convex surfaces with outward-facing normals).

In this study, we compute the curvature based on Rusinkiewicz estimation²², which may be thought of as an extension of common methods, such as the curvature presented in Knutsen et al.³⁶ and used by subsequent authors like Garcia et al. ³⁷, for the purpose of estimating per-vertex normals by averaging adjacent per-face normals. This algorithm uses the “Voronoi area” weighting which can produce more accurate normal estimates of curvature than other weighting methods for triangles of varying sizes and shapes. In this algorithm, the per-face (per-triangle) curvature tensor is first computed by its three well-defined directions (the edges) together with the differences in normals in those directions (computed from the per-vertex normals). Then, the algorithm performs a coordinate system transformation for converting the curvature tensor to the vertex coordinate frame. Eventually, a per-face coefficient is applied to allow to weight the face curvature around each vertex.

Mean curvature of a vertex is defined by the average of the two principal curvatures (the maximal and minimal curvatures) of the vertex, and the principal curvatures are the eigenvalues of the vertex normal curvature tensor computed by Rusinkiewicz estimation:

where K1 and K2 are the eigenvalues and (u′,v′) are the principal directions, which are the directions in which the normal curvature reaches its minimum and maximum. Since surface curvature is useful to describe spatial variations in folding, thus for the overall folding complexity comparison, we first compute a dimensionless mean curvature by multiplying the mean curvature by the square root of the surface area (Karea), where K is the mean curvature. Then we calculate the average (across all vertices on the mesh surface) of the absolute value of dimensionless mean curvatures (at a vertex on the mesh surface) for each simulated surface. In the remainder of the manuscript, we simply use the term curvature for the sake of clarity.

Curvature-based features do not provide complete description of the folding patterns. In order to describe globally the folding complexity by considering the depth and wideness of the cortical folding, we also use the surface-based three-dimensional gyrification index (3D GI). It is a global measurement which is defined as the ratio of the cortical surface area to the area of its smooth “convex hull” (the minimum surface area needed to completely enclose the brain)²³:

To get the convex hull, we scale the initial smooth surface in three dimensions so that the three-dimensional lengths of the convex hull are equal to those of the simulated cortical surface.

Sulcal depth can be used as a quantitative marker of cortical morphology³⁸. Several approaches have been proposed to compute the sulcal depth^39–41 but a well-defined computation of depth remains an open question. In this work, we make use of an intuitive approach to calculate the sulcal depth by using the distance between the deformed mesh surface and the corresponding convex hull. Specifically, for each surface vertex of the deformed mesh, we find the intersection point on the convex hull by using the vector determined by the corresponding vertex of the initial mesh and this vertex and the method of traversing all triangles on the convex hull. Then we compute the distance between each surface vertex of the deformed mesh and its corresponding intersection point on the convex hull.

For the purpose of describing and comparing the direction of the folds on the simulated surfaces, we calculate the angle between the gradient of Fiedler vectors²⁴ and the principal directions of curvatures²⁵, which helps to understand whether the folds are isotropic. The Fiedler vector is the first non-constant eigenfunction of Laplace-Beltrami operator, represented by ϕ1 in Eq. ⁸²⁴. The Laplace-Beltrami operator is defined as ΔM=div·∇M, where M is a Riemannian manifold. The eigenvalues of -ΔM are λ0=0≤λ1≤… and ϕ0, ϕ1, … are associated orthonormal basis of eigenfunctions, which satisfy

The Fiedler vector allows to describe the longitudinal extension of surfaces^24,42–44. The Fiedler’s extrema are the most distant points²⁴, and its contour lines are slices in the elongation axis. The gradient of the Fiedler vector that is perpendicular to the contour lines gives the direction of elongation. The principal directions of curvatures are the corresponding eigenvectors of the principal curvatures (the eigenvalues of the Weingarten matrix). Based on the local scalar product between the gradient of the Fiedler vector and the principal directions of curvatures, we can obtain the folds angle.

In order to compare quantitatively the uniformity of the angular distribution of folds, we use the Kullback–Leibler (KL) divergence. The KL divergence, also called relative entropy, is used to measure how one probability distribution is different from a second reference probability distribution. For two discrete probability distributions P and Q defined on the same probability space, the KL divergence from P to Q is defined to be

For angular uniformity calculations, P corresponds to the fold angular distribution on the folded surface, Q represents the theoretically uniform distribution of fold angles.