1 Introduction

The morphological complexity of the mammalian cerebral cortex has fascinated scientists for generations, with cortices across and within species exhibiting a large diversity of shapes and sizes. Such diversity is not arbitrary, however. The mammalian brain folds into stereotypical, hierarchically-organized structures such as lobes and major gyri. In fact, qualitative and quantitative regularities in cortical scaling have often been suggested and observed [1, 2, 3, 4]. More specifically, through modelling the mechanism of cortical folding from a statistical physics approach, we have previously derived a theoretical scaling law relating pial surface area A_t, exposed surface area A_e [the exposed surface area can be thought of as the surface area of a piece of cling film wrapped around the brain. Mathematically, for the remaining paper it is the convex hull of the brain surface], and average cortical thickness T:

This scaling law, relating powers of cortical thickness and surface area metrics, was shown to be valid across mammalian species [5] and within the human species [6], as well as to the structures and substructures of individual brains [7, 8]. Notably, across all these cases the dimensionless offset k is shown to be near invariant. However, this universality, presumed to be arising from universal physical principles and evolutionarily conserved biomechanics, says little about what that cortical shape actually is, beyond a constraint binding three morphometric parameters. In this paper, we take Eqn. (1) as an empirical starting point to create a new and hierarchical way of expressing cortical shape. Specifically, we introduce a coarse-graining procedure that renders the cortex at different spatial scales, or resolutions. We show that coarse-grained primate cortices at each spatial scale can be understood as approximations of the same universal self-similar archetypal form, of which the observed scaling law (Eqn. 1) can then be shown to be a direct consequence.

Besides revealing a symmetry in nature hidden under much apparent complexity, our results indicate a conservation of morphological relationships across evolution. We will show that these results further provide us with a new and powerful tool to express and analyse cortical morphology. As an example, we will calculate the effects of human ageing across spatial scales and show that the effects are highly scale-dependent.

1.1 Mathematical background

The universal scaling law (Eqn. 1) can be rewritten in a suggestive way

where the is a fundamental area element that defines the threshold between gyrencephaly (folded cortex) and lissencephaly (smooth cortex) when A_t = A_e = A₀. For a constant k the value of A₀ is a multiple of T², indicating that cortical thickness determines the size of the smallest possible gyri and sulci.

This re-writing highlights a new perspective, or interpretation of the scaling law: it now suggests a relationship between intrinsic and extrinsic measures of cortical size [given the folded laminar structure of the cortex, areas are the more natural way of measuring its “size”]: A_t and A_e respectively, measured in units of A₀. This is reminiscent of fractal scaling [9], where a complex shape reveals ever smaller levels of self-similar detail as it is probed in ever smaller scales (or equivalently, higher resolutions), represented here by A₀. The scaling, or power exponent between the measured intrinsic and extrinsic sizes is the so-called fractal dimension.

Although actual fractals are mathematical abstractions, they can often be defined as the limit of iterative processes. Many structures in nature, and in particular biology [10, 11, 12, 13, 14, 15, 16], are good approximations of a fractal. Equation (2) is suggestive, but not proof, that cortices are among these forms, with a fractal dimension of 1.25 × 2 = 2.5 (the factor 2 being the topological dimension of areas). Indeed, fractal scaling for various aspects of cortical morphology has often been postulated [17, 18], with a number of recent papers making use of MRI data [19, 20, 21]. Most recently published estimates of fractal dimension for the whole cortex are indeed close to 2.5 [22, 23, 24].

Here, for the first time, we propose to directly construct morphologically plausible realisations of cortices at any specified spatial scale, or resolution. This is achieved through a coarse-graining method that removes morphological details smaller than a specified scale while preserving surface integrity. For example, at a set scale of 3 mm, sulcal walls that are less than 3 mm apart would be removed, and the neighbouring gyri would be fused. This method is a new systematic way of obtaining shape properties from the cortex in terms of a sequence of morphometric measurements as spatial scale varies. By examining how areas scale across coarse-grained versions of actual primate cortices, we will be able to directly verify cortical self-similarity.

2 Method

2.1 Coarse-graining method

As a starting point for a coarse-graining method, we suggest to turn to a well-established method that measures fractal dimension of objects: the so-called box-counting algorithm [25, 26]. Briefly, this algorithm fills the object of interest (the cortex in our case) with boxes, or voxels of increasingly larger sizes and counts the number of boxes in the object as a function of box size. As the box size increases, the number of boxes decreases; and in a log-log plot, the slope of this relationship indicates the fractal dimension of the object. In our case, this method would not only provide us with the fractal dimension of the cortex, but, with increasing box size, the filled cortex would also contain less and less detail of the folded cortex. Intuitively, with increasing box size, the smaller details below the resolution of a single box would disappear first, and increasingly larger details will follow – precisely what we require from a coarse-graining method. We therefore propose to expand the traditional box-counting method beyond its use to measure fractal dimension, but to analyse the reconstructed cortices as different realisations of the original cortex at the specified spatial scale.

Concretely, our proposed method requires the bounding pial and white matter surfaces of the cortical ribbon as input. We obtained these surfaces based on reconstructions from magnetic resonance imaging data in 11 different primate species. Algorithmically, we then segment the space between the original pial and white matter surfaces into a 3D grid of boxes of the desired scale λ, where each box is a cube of dimensions λ × λ × λ. We also term the 3D grid of cubes “voxelisation”, as it effectively captures the cerebral cortex as voxels in 3D space (Fig. 1B bottom row). At any given scale, or voxel size, this process effectively erases morphological features (folds) that are smaller than the cube size. Visually, increasing the voxel size appears as if the cortex is “melting” and “thickening” (Fig. 1, and videos: https://bit.ly/3CDoqZQ).

Figure 1:

A more technical and detailed description and discussion of the algorithm is provided in Suppl. S1. Note this method has also no direct dependency on the original MR image resolution, as the inputs are smooth grey and white matter surface meshes reconstructed from the images using strong (bio-)physical assumptions and therefore containing more fine-grained spatial information than the raw images (see also Suppl. S2).

2.2 Rescaling coarse-grained outputs for analysis

Morphological properties, such as cortical thicknesses measured in our “melted” brains are to be understood as a thickness relative to the size of the brain. Therefore, to analyse the scaling behaviour of the different coarse-grained realisations of the same brain, we apply an isometric rescaling process that leaves all dimensionless shape properties unaffected (more details in Suppl. S3.1). Conceptually, this process fixes the voxel size, and instead resizes the surfaces relative to the voxel size, which ensures that we can compare the coarse-grained realisations to the original cortices, and test if the former, like the latter, also scale according to Eqn. (1). Resizing, or more precisely, shrinking the cortical surface is mathematically equivalent to increasing the box size in our coarse-graining method. Both achieved an erasure of folding details below a certain threshold. After rescaling, as an example, the cortical thickness also shrinks with increasing levels of coarse-graining, and never exceeds the thickness measured at native scale.

2.3 Independent morphological measures of shape

To better characterise the coarse-grained cortices in terms of their similarity in offset, we use a previously introduced [27] set of independent measures, K, I and S, that summarise the morphometry of the cortex in a natural and statistically robust way. In this framework, isometrically scaled copies of the same morphometry all map onto a line along the I = log A_t + log A_e + log T² direction, which is perpendicular to a K × S plane that fully summarises their shape. K = is the direction defined by the offset k of the scaling law Eqn. (1), while direction captures the remaining information about shape, and can be regarded as a simple measure of morphological complexity. Ref [27] provides a detailed derivation and demonstration of the superior sensitivity and specificity of these new morphometric measures. The advantage of using this framework here is that we can assess the offset K (and shape term S) without interference by isometric size effects, including any re-scaling procedures.

2.4 Data and processing

With the exception of the marmoset data, all other cortical surface reconstructions were based on healthy individual brains.

2.4.1 Human data

To study healthy human adults, we used the Human Connectome Project (HCP) MRI data, available at https://db.humanconnectome.org/ [28], obtained using a 3T Siemens Skyra scanner with 0.7 mm isotropic voxel size. We used the HCP minimally pre-processed FreeSurfer data output, which provided the pial and white matter surface meshes we required. We selected five random subjects in the age category 22–25 y.o. and show one example subject (103414) in the main text, and the remaining subjects in Suppl. S3.

To study the alterations associated with human ageing, we used T1 and T2 weighted MRI brain scans from The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) dataset (available at http://www.mrc-cbu.cam.ac.uk/datasets/camcan/ [29, 30]). Cam-CAN used a 3T Siemens TIM Trio System with 1 mm isotropic voxel size (for more details see [29, 30]). From the Cam-CAN dataset we retained 644 subjects that successfully completed preprocessing (with Freesurfer recon-all) without errors. From these subjects we selected all subjects between the ages of 17–25 inclusive (forming the 20 y.o. cohort, n=27); we also selected all subjects between the ages of 77–85 inclusive (forming the 80 y.o. cohort, n=86).

To confirm the ageing results, we also obtained an independent dataset from the Nathan Kline Institute (NKI)/Rockland sample [31] (http://fcon_1000.projects.nitrc.org/indi/pro/nki.html) using the same procedure as described for the CamCAN dataset.

The MR images of both CamCAN and NKI datasets were first preprocessed by the FreeSurfer 6.0 pipeline recon-all, which extracts the grey-white matter boundary as well as the pial surface. These boundaries were then quality checked by visual inspection for particularly the young and old cohorts and manually corrected where needed.

For all three datasets, we obtained the pial and white matter surfaces for further analysis. In the current work, the analysis is always hemisphere based, as in our previous work [5, 6]. We did not perform a more regionalised analysis, which is also possible [7, 8]. Future work using the principle demonstrated here can be directly extended to derive regionalised measures across scales.

2.4.2 Non-human primate data

Macaque

Rhesus Macaque MRI scans were carried out at the Newcastle University Comparative Biology Centre. Macaques were trained to be scanned while awake and sat in a primate chair. Both T1 weighted MP-RAGE and T2 weighted RARE sequences were acquired, using a vertical MRI scanner (Biospec 4.7 Tesla, Bruker Biospin, Ettlingen, Germany).

Scans were processed using a custom macaque MRI pipeline, incorporating ANTs, SPM, FreeSurfer and FSL. Briefly, this involved the creation of precursor mask in SPM, denoising (ANTs DenoiseImage) and debiasing (ANTs N4BiasField-Correction, and Human Connectome BiasField-Correction script), creation of a final mask in SPM and then processing in FreeSurfer (using a modified version of the standard FreeSurfer processing pipeline). Reconstructed pial and white matter surfaces were visually quality controlled in conjunction with the MR images.

Marmoset

The marmoset MRI structural scan was collected as part of the development of the NIH marmoset brain atlas [32, 33]. Data was collected ex vivo from a 4.5-year-old male marmoset using a T2-star weighted 3D FLASH sequence using a horizontal MRI scanner (Biospec 7 Tesla, Bruker Biospin, Ettlingen, Germany).

A total of ten scans were collected and averaged into one final image. In combination with scans from other modalities, cortical boundaries were manually delineated on each coronal slice. Boundaries were then refined through comparisons with other atlases. Volumetric data was then converted to surfaces using a custom pipeline involving an intermediate generation of high-resolution mesh data [34], decimation [35], and remeshing. Reconstructed pial and white matter surfaces were visually quality controlled in conjunction with the MR images.

Other non-human primates (NHPs)

The remaining NHP MRIs and subsequent brain surface extraction is detailed in [36, 37], and provided to the authors in a processed format. Briefly, a range of specialised scanners were used to acquire optimal images for each species. FreeSurfer 6.0 with some modifications was used for surface reconstruction, complemented by FSL, ANTS, and Matlab. All surfaces were visually inspected for accuracy and consistency across datasets.

2.4.3 Comparative neuroanatomy data

The comparative neuroanatomy dataset for different mammalian species is the same as previously published [5]. Note that for this dataset, we only had numerical values for the total and exposed surface area, as well a average cortical thickness estimates. We did not perform any analysis across scales in this dataset (hence surfaces were not required), but only used it as a reference dataset.

2.5 Statistical analyses

Briefly, linear regression is used in either a mixed-effect model to capture effects across individuals and species, or in simple fixed-effect settings to estimate regression slopes to obtain fractal dimension.

In the final part of Results, we analyse the effect between a group of 20-year-olds and 80-year-olds. Effect size is calculated as Cohen’s D between the two groups.

Throughout the paper, statistical significance is not a crucial argument, and we report p-values only for reference and completeness.

More details can be seen in the analysis code, where the reader can directly reproduce all main result figures.

2.6 Code and data availability

The code for coarse-graining has been integrated into our MATLAB toolbox Cortical Folding Analysis Tools: https://github.com/cnnp-lab/CorticalFoldingAnalysisTools, include a graphical user interface in the latest version. Users will also see the latest updates in this repository.

The data and analysis code underpinning this paper will be published on GitHub upon final decision: https://github.com/cnnp-lab/2024_Folding_scales/ to allow reproduction of our main results.

3 Results

3.1 All primate brains follow the same scaling law across spatial scales

We have analysed cortices of 11 different primate species: aotus, cebus, chimpanzee, colobus, galago, human, lagothrix, lophocebus, macaque, marmoset, and pithecia, and various cohorts of human subjects. We applied our coarse-graining procedure to their pial and white matter surfaces, and empirically determined (i) that all species followed a power law (linear regression R² > 0.999 for all species); (ii) the slope of said power law is α = 1.255 on a group level (CI: [1.254 1.256]) using linear mixed effect modelling (see Fig. 2 for visualisation, and Suppl. S3 for a detailed breakdown by species); and importantly, (iii) all species also show a similar offset log k ≈ -0.65263, with a standard deviation of intercept across species estimated at 0.02 from linear mixed effect modelling in log k.

Figure 2:

Taken separately, the scaling for each species is proof that their cortices are self-similar with the same scaling: they each approximate a fractal with fractal dimension d_f = 2.5. Considering all species together, different species also overlap substantially (similar offset), and only differ from each other in the range of scales over which the approximation is valid (see Suppl. S3.2 and S4). Thus, as Fig. 2 illustrates, the data supports a universal scaling law across primate species and spatial scales:

with k = 0.2277.

3.2 Primate brains at different spatial scales are morphometrically similar to each other and other mammalian species

To better characterise the coarse-grained cortices in terms of their similarity in offset, we use a set of independent morphometric measures, K, I and S, that summarise the morphometry of the cortex in a natural and statistically robust way. We can therefore assess the offset K (and shape term S) without interference by the isometric size or rescaling.

We can measure K and S for any object, but a fuller expression is captured by the trajectory of said object as a function of coarse-graining in the K × S plane. This is a very convenient and informative way of summarising an object: self-similar objects correspond to straight trajectories as the K × S plane is in log-log space. In particular, objects without any folds or protrusion (i.e., convex, such as the box with finite thickness in Fig. 3 A) correspond to the line , as A_e = A_t for all levels of coarse-graining. Horizontal trajectories (constant K) represent fractal objects with fractal dimension d_f = 2.5 (Fig. 3 B). And finally, in the K × S plane, a group of objects can said to be “universal” when their trajectories overlap, so that they can all be regarded as coarse-grained versions of one another (Fig. 3 C).

Figure 3:

Primate cortices (Fig. 3 D) display a nearly invariant K in all cases. But, over all levels of coarse-graining, K also remains near-invariant in all trajectories as S decreases, resulting in a set of horizontal lines that largely overlap with each other and other mammalian species. The variance in K across scales and all 11 species is < 0.01, which is at least an order of magnitude lower than the variance in S. Primate brains therefore have all three characteristics of self-similarity, fractality (with d_f = 2.5), and universality (invariant K for all scales and species) at the same time.

Thus, coarse-grained primate cortices are morphometrically similar to, and in terms of the universal law, ‘as valid as’ actual existing mammalian cortices. [Note, of course, that the coarse-grained brain surfaces are an output of our algorithm alone and not to be directly/naively likened to actual brain surfaces, e.g. in terms of the location or shape of the folds. Our comparisons here between coarse-grained brains and actual brains is purely on the level of morphometrics across the whole cortex.] In contrast, we tested various non-brain objects, and while e.g. the walnut, and bell pepper form (partially) straight lines, they vary in both K and S (see Suppl. S5). These objects may have a fractal regime, but their fractal dimension is not 2.5, nor are they similar to primate or mammalian brains in terms of K or S. Furthermore, supplementary S6 underscores the algorithmic and statistical robustness of these results using multiple realisations of the coarse-graining procedure on the same object. In this framework, our main result can thus be expressed simply: for all the cortices we analysed, and for none of the non-cortices, coarse-graining will leave K largely unaffected, while morphological complexity S will decrease.

3.3 Morphometric measures as functions of scale reveal scale-specific effects of ageing

In the final part of our work, we show how our algorithm and the associated new understanding of brain morphology may become useful in applications. As an example, we will focus on how the ageing process affects human cortical morphology across scales. In Fig. 4 A, we compare the total surface area A_t(λ) as a function of scale λ for a young (20 year old) vs. an old (80 year old) group of human brains. The difference in A_t between the groups takes a U shape, and the strongest effect is seen at approximately 2 millimeter (greatest effect size of -8.635 seen at scale 2.188 mm), where older subjects have higher A_t. For scales over ~ 5 mm, and under ~ 0.5 mm the differences become relatively small, suggesting the ageing process has less effect on the largest and smallest cortical morphological features. Finally, we reproduced these results in an independent dataset in Suppl. S7.2.

Figure 4:

In this particular example, the scale-dependency of morphological measures can be visually and intuitively understood by looking at the reconstructed surfaces at each scale: Fig. 4 B shows some coronal slices of the cortical surface. At scale 0.27 mm, the gyri in the younger subjects are densely packed, but the older subjects show the expected widening between gyral walls and decrease in gyral surface area at the crown (see e.g. [38, 26] for recent investigations and references therein). At 1.86 mm, the younger cortices have already partially “melted”, erasing most small sulci between the densely packed gyri. In the older humans, however, the gyri are less dense, the sulci more open, and thus, at 1.86 mm, most gyri and sulci have not been erased yet. At scale 7.94 mm, both young and old brains have “melted” down to similarly near-lissencephalic cortices. For a more detailed multiscale investigation over the entire human lifespan, please refer to our new preprint [39]

More broadly, one can regard the melting process as a way of determining how cortical area is allocated across different scales. As the cortex ‘melts’, the contributions to the total area from features smaller than the cut-off scale are eliminated. In the example in Fig. 4, for instance, we can say about half (10^5-4.7 = 10^0.3 ≈ 2) the total area in the 80 y.o. cortices is present in features smaller than 4 mm.

4 Discussion

We have devised a new way of expressing the morphology of the mammalian cerebral cortex, as the flow in the values of morphometric measures over a range of spatial scales. This was achieved by coarse-graining cortical surfaces, erasing morphological features smaller than the specified scale, while preserving surface integrity. After applying this method to the cortices of 11 primate species, we have shown that all these diverse cortices are approximations of the same archetypal self-similar (fractal) shape. Most of their morphological diversity can be ascribed to the species-specific ranges of spatial scales over which each approximation is valid, with the smallest scale being an invariant multiple of cortical thickness. This was a proof-by-construction of fractality, a step beyond the usual box-counting approach, which also yields scale-dependent morphometrics. As a proof-of-principle we showed that healthy human ageing has highly scale-dependent effects in a range of morphometrics.

4.1 Advantages and advances

Compared to previous literature, we can summarise our main contribution and advance as follows: (i) We are showing for the first time that representative primate species follow the exact same fractal scaling - as opposed to previous work showing that they have a similar fractal dimension [40, 41], i.e. slope, but not necessarily the same offset, as previous methods had no consistent way of comparing offsets. (ii) Previous work could also not show direct agreement in morphometrics between the coarse-grained brains of primate species and other non-primate mammalian species. (iii) Demonstrating in proof-of-principle that multiscale morphometrics, in practice, can have much larger effect sizes for classification applications. This moves beyond our previous work where we only showed the scaling law across [4] and within species [6], but all on one (native) scale with comparable effect sizes for classification applications [6].

In simple terms: we know that objects can have the same fractal dimension, but differ greatly in a range of other shape properties. However, we demonstrate here, that representative primate brains and mammalian brain indeed share a range of other key shape properties, on top of agreeing in fractal dimension. This suggests a universal blueprint for mammalian brain shape and a common set of mechanisms governing cortical folding. As a practical additional outcome of our study, we could show that our novel method of deriving multiscale metrics can differentiate subtle morphological changes much better (4 times the effect size) than the metrics we have been using so far at a single native scale.

Expressing cortical morphology as a function of scale is more detailed than a list of summary morphometric measures, and more informative than the mere listing of every sulcus and gyrus. We propose this new syntax as the basis for a more rigorous characterisation of brain morphology and morphological changes. A clear advantage is that some biological processes may only act on a specific spatial scale, leaving other scales untouched (ageing in our example). By disambiguation of the spatial scale, it allows for an extra dimension of understanding, and tracking of biological processes. In the future, this approach can also be extended to cortical development and to various degenerative [6] and congenital [27] neuropathic conditions, especially if combined with a regionalised version of this method applied to specific cortical regions [7] or local patches [8].

4.2 Implications of universality

Empirically, the main result of this paper is the demonstration of a universal self-similar scaling (Eqn. 3) for primate, and presumably mammalian cortices. There are two aspects of this universality: first that for each and every cortex the value of K remains the same for all scales as one removes substructures smaller than a varying length scale (or equivalently, that the fractal dimension is almost exactly 2.5 in all cases). Second, that the value of K for the cortices of different species is approximately the same, as previously observed (Eqn. 1) across species [5] and individuals [6]. One could imagine a set of objects for which one but not the other aspect of the universality in K holds true. However, the fact that both universalities hold true is significant. It suggests the existence of a single highly conserved mechanism for cortical folding, operating on all length scales self-similarly with only a few morphological degrees of freedom. It also hints at the possibility of deriving cortical scaling from some variational principle. Finally, this dual universality is also a more stringent test for existing and future models of cortical gyrification mechanisms at relevant scales, and one that moreover is applicable to individual cortices. For example, any models that explicitly simulate a cortical surface as an output could be directly coarse-grained with our method and the morphological trajectories can be compared with those of actual human and primate cortices. The simulated cortices would only be ‘valid’ in terms of the dual universality, if it also produces the same morphological trajectories. [Note, we do not suggest to directly compare coarse-grained brain surfaces with actual biological brain surfaces. As we noted earlier, the coarse-grained brain surfaces are an output of our algorithm alone and are not to be directly/naively likened to actual brain surfaces, e.g. in terms of the location or shape of the folds. Our comparisons here between coarse-grained brains and actual brains is purely on the level of morphometrics across the whole cortex.]

The scaling itself does not imply or favour any particular proposed gyrification model (ours [5] included), and all results in this paper are agnostic about this choice. Indeed, our previously proposed model [5] for cortical gyrification is very simple, assuming only a self-avoiding cortex of finite thickness experiencing pressures (e.g. exerted by white matter pulling, or by CSF pressure).The offset K, or ‘tension term’, precisely relates to these pressures, leading us to speculate that subtle changes in K correlate with changes in white matter property [6, 27]. In the same vein of speculation, the scale-dependence of K shown in this work might therefore be related to different types of white matter that span different length scales, such as superficial vs. deep white matter, or U-fibres vs. major tracts. However, there are also challenges to the axonal tension hypothesis [42]. Indeed, white matter tension differentials in the developed brain may not explain location of folds, but instead white matter tension may contribute to a whole-brain scale ‘pressure’ during development that drives the folding process overall. Aside from speculations about the biological interpretation, the simplicity of the highlighted scaling law parallels many complex phenomena in nature that display simple and universal scaling that can be derived from first principles [43, 44, 45]. In addition, recent results suggest simplicity and symmetry are generically favoured on statisticalensemble grounds by evolution [46]. Our model correctly predicts the scaling law (Eqn. 1), but a more complete explanation for cortical gyrification is probably far more complex [47] than can be accounted by such a simple model.

One specific example of said complexity are the exact patterns, locations, depth, and features of gyri and sulci. We know such patterns to be, for example, variable but also somewhat heritable in humans, whilst in macaques such patterns are relatively preserved across the species. Our work does not explain any of these observations, nor are the coarse-grained versions of human brains supposed to exactly resemble the location/pattern/features of gyri and sulci of other primates. The similarity we highlighted here are on the level of summary metrics, and our goal was to highlight the universality in such metrics to point towards highly conserved quantities and mechanisms.

4.3 Biological plausibility and implications

The observation that with increasing voxel sizes, the coarse-grained cortices tend to be smoother and thicker is particularly interesting: the scaling law in Eq. 3 can be understood as thicker cortices (T) form larger folds (or are smoother i.e. less surface area A_t) when brain size is kept constant (A_e). This way of understanding has also been vividly illustrated by using the analogy of forming paper balls with papers of varying thickness in Mota et al.[4]: to achieve the same size of a paper ball (A_e), the one that uses thicker paper (T) will show larger folds (or is smoother i.e. less surface area A_t) than the one using thinner paper. [Note that these observations should not be interpreted in absolute terms, smoother brains are of course not always thicker (measured in millimeters) than folded brains, neither in the real-world, nor in our outputs. These statements are formed under the assumption of the same isometric size of the cortex, or the paper ball.] The scaling law can therefore be understood as a physically and biologically plausible statement and our algorithm yields results in line with the scaling law.

More broadly, the interaction between brain development and evolution may also benefit from a scale-specific understanding. This may be important in elucidating what in cortical morphology is selected for by evolution, what is determined by physics; what is specified by genes, and what is emergent. For example, one can estimate the number of structural features at each scale (A_t(λ) as multiples of A₀(λ)), and it will be interesting to correlate this number to other quantifiers of cortical structure, such as number of neurons [4], number of functional areas [48] or the number of cortical columns [49], possibly over different stages of development. Generally, the larger the cortical feature (i.e. from gyri to functional areas to lobes to hemispheres), the earlier during development it appears [50, 51], and more broadly it is conserved over kinship [52] and phylogeny [53, 54]. It thus seems likely that comparative neuroanatomical methods [55, 56] may be directly used to identify and contrast structures in coarse-grained cortices of more highly gyrified species with their analogues in less-gyrified species. One could then perhaps specify when evolution conserves and when it invents old and new cortical features.

From an application perspective, our final result illustrates clearly that the surface area difference between older and younger subjects at “native” scale (i.e. original freesurfer surfaces) is negligible (effect size smaller than 2 standard deviations). However, in our analysis across scales, there is a clear optimal scale at ~ 2 mm where the effect size is maximised between older and younger subjects (effect size is -8 standard deviations). For most classification applications in biology and medicine, the increased effect size and hence separability of groups in the scale-dependent morphometrics represent a huge advance over the native scale.

4.4 Outlook

Our work here was limited to summary descriptors of entire cortical hemispheres, but future work will explore extensions of these methods to lobes and cortical areas, similarly to [7, 8]. This will generate precise characterisations of the morphological differences between phylae and across developmental stages, and perhaps pinpoint the time and location of morphological changes leading to congenital and neurodegenerative conditions. Ultimately, we hope this new framework for expressing and analysing cortical morphology, besides revealing a hitherto hidden regularity of nature, can become a powerful tool to characterise and compare cortices of different species and individuals, across development and ageing, and across health and disease.

Acknowledgements

We thank members of the Computational Neurology, Neuroscience & Psychiatry Lab (www.cnnp-lab.com) for discussions on the analysis and manuscript, and Dirk Jan Ardesch and Martijn van den Heuvel for helpful discussions and NHP brain surface data. P.N.T. and Y.W. are both supported by UKRI Future Leaders Fellowships (MR/T04294X/1, MR/V026569/1); Y.W. and K.L are further supported by the EPSRC (EP/Y016009/1, EP/L015358/1). B. Mota is supported by Fundaçao Serrapilheira Institute (grant Serra-1709-16981) and CNPq (PQ 2017 312837/2017-8).

Supplementary

S1 Coarse-graining algorithm

The core coarse-graining algorithm underpinning our analyses takes pial and white matter surfaces from the cortical surface reconstruction as inputs and at any specified spatial scale λ “voxelises” the gray matter ribbon at the specified resolution. To achieve this, we set up a 3D voxel grid, where each grid voxel is of size λ × λ × λ.

In detail, we assign all voxels in the grid with at least four corners inside the original pial surface to the pial voxelization. This process allows the exposed surface to remain approximately constant with increasing voxel sizes. A constant exposed surface is desirable, as we only want to gradually “melt” and fuse the gyri, but not grow the bounding/exposed surface as well. We want the extrinsic area to remain approximately constant as we decrease the instrinsic area via coarse-graining; it is like generating iterates of a Koch curve in reverse, from more to less detailed, by increasing the length of smallest line segment.

We then assign voxels with all eight corners inside the original white matter surface to the white matter voxelization. This is to ensure integrity of the white matter, as otherwise white matter voxels in gyri may become detached from the core white matter, and thus artificially increase white matter surface area. Indeed the main results of the paper are not very sensitive to this decision using all eight corners, vs. e.g. only four corners, as we do not directly use white matter surface area for the scaling law measurements. However, we still maintained this choice in case future work wants to make use of the white matter voxelisations or derivative measures.

Finally, subtracting the white matter voxelisation from the pial voxelisation, we obtain the grey matter voxelisation (Fig. 1B bottom row). The pial surface voxelisation is designed such that neighbouring gyri can fuse if their separation is smaller than the scale of interest, whilst not growing the cortex outwards. The white matter surface voxelisation allows the walls of each gyrus to thicken and fuse inwards. Visually, this process looks as if the cortex is “thickening” inwards and smoothing on the outside.

From the voxelised cortex at each scale A, we can then obtain coarse-grained pial and white matter surfaces, and extract estimates of global morphometric measures such as the average cortical thickness T(λ), the total cortical surface area A_t(λ), and the exposed surface area A_e(λ). For each λ, we derived an outer isosurface equal to 0.5 for the pial voxelisation (voxel values are 1 within the voxelised pial surface and 0 otherwise). This isosurface is then defined as the coarse-grained pial surface at this scale. The surface area of this isosurface is used as A_t(λ). The exposed surface area A_e(λ) is subsequently derived from the convex hull of the pial isosurface. Finally, the cortical thickness T(λ) is estimated as , where V_G(λ), which is the estimated grey matter volume, derived from the number of grey matter voxels, is multiplied by the voxel volume. These global morphometric measures are summary statistics of the brain at each particular scale, capturing information about both their intrinsic geometry (A_t(λ), T(λ)) and extrinsic geometry (A_e(λ)). In this manner, for each cortex, we obtain not just one set of summary morphometric measures, but rather an set measures as functions of λ. A detailed walk-through of the coarse-graining algorithm and estimation of morphometric measures is provided on Github: https://github.com/cnnp-lab/CorticalFoldingAnalysisTools/blob/master/Scales/fastEstimateScale.m.

Note that although this process is inspired by the box-counting algorithm, it is different from box-counting and related convolution-based algorithms: we do not simply apply successive convolutions with an increasing kernel size (or equivalent), which effectively would achieve a uniform and spatially isometric “smearing” of the original cortical ribbon to a given scale, rather than a targeted erasure of surface details smaller than said scale. As a result, this method yields well-defined and well-behaved white and gray matter surfaces. Thus, one can apply all the usual analytical tools to these realisations that are applicable to actual cortices. Note that our approach is more comparable in its principles to the calculation of the outer smoothed pial surface in FreeSurfer [57], which utilises a dilation and erosion convolution to effectively erase details below a certain scale. Of course, our proposed procedure is not the only conceivable way to erase morphological details below a given scale; and we are actively working on related algorithms that are also computationally cheaper. Nevertheless, the current version requires no fine-tuning, is computationally feasible and conceptually simple, thus making it a natural choice for introducing the methodology and approach.

Given how broadly it has been verified, we expect the observed universality in cortical selfsimilar scaling to be robust to the details of the coarse-graining algorithm. It would be very informative to test this proposition in future. For example, an alternative method inspired by [58] could be implemented, eschewing voxelization and dealing only with the flow of nested surfaces with self- and mutual-avoidance explicitly implemented. Going in the opposite direction, more detailed models for the mechanisms of cortical folding (see e.g. [59]) can be regarded as a type of reverse melting, and could perhaps be implemented and tested in a similar fashion as fine-graining procedures.

S2 Influence of original image resolution

In this supplementary section, we want to demonstrate the relatively weak effect of the original image resolution on our analysis outputs. To this end, we used five example HCP subjects, who were scanned at 0.7mm isotropic image resolution, and downsampled their images to 1mm isotropic images. We then proceeded with our analysis using two freesurfer outputs. (1) the HCP freesurfer pipeline output optimised for the 0.7mm resolution, and (2) a standard freesurfer pipeline output on the 1mm downsampled images. We proceeded with these two sets of surfaces in our analysis and show the resulting morphology measures in Fig. S2.1.

We observe a relatively weak difference between these two sets of inputs/surfaces. Both sets largely follow the same trajectory across scales. Especially in the exposed area, the within and between subject differences are noticeably larger than between image resolutions. In total area and thickness, some small but systematic differences are seen in all subjects between the scales of 1–2mm. These are most likely differences in the cortical morphology reconstruction in the freesurfer surfaces using the two different resolution images as input. However, in none of these resulting morphology measures do we see an artifact specifically at 0.7 or 1mm. Any differences between image resolutions only result in subtle changes in the freesurfer meshes, which are smooth. Our analysis method therefore has no direct dependency on the image resolution, as our inputs are these smooth freesurfer meshes.

Figure S2.1:

S3 Scaling properties by species

S3.1 Obtaining the scaling law

From the coarse-grain procedure (Fig. 1), we obtain surface meshes for the pial and white matter surface at each spatial scale. From those, we derive exposed area, total pial surface area, and average cortical thickness as described in Methods. However, the coarse-graining procedure, by itself, barely changes the exposed surface area, and only changes the total surface area minimally (Fig. S3.1 A). As the voxel size changes at each scale, it only starts to affect the exposed area at very large scales, where effectively the voxels are no longer a good description of the shape of the skull. Thus, if we scatter the raw data points from the coarse-graining procedure in the plane of the scaling law, there is barely any variance in the data (Fig. S3.1 B). Even after zooming in, we see an almost vertical line (Fig. S3.1 B). This is expected and is not evidence for or against our hypothesis.

Instead, we need to transform the data into a perspective that makes sense to be be seen in the scaling law plane. Instead of thinking of the coarse-graining procedure as a process that re-renders the cortex at increasing voxel sizes, we can instead think of the procedure as rescaling the original cortical mesh into increasingly smaller sizes, and re-rendering it with the same fixed voxel size. The analogy using Britain’s coastline is that instead of measuring the coastline with increasingly smaller rulers, we resize the map of the coastline to increasingly smaller sizes, but keep the size of the ruler the same. Both procedures are equivalent and produce the same fractal dimension.

Figure S3.1:

To achieve this procedure, we re-scaled our measurements of A_e, T and A_t by the voxel size (λ) and a fixed factor (l_r):

As we used isometric cubes as voxels, the voxel size refers to the length of a single side. E.g. λ = 1 if we used a isometric 1 × 1 × 1mm³ voxel.

l_r is a fixed factor for each cortical hemisphere, and does not change with λ. We use it to systematically shift all the data points within a range, such that the re-scaled quantities are not larger than those from the original cortical meshes. One can easily verify that l_r will not change the slope or offset of any scaling law, but simply represents a constant shift to all data points. In our data, some of the re-scaled quantities would indeed be larger than those from the original cortical meshes, as the voxel size we choose is limited at the smaller end only by computational resources. In other words, we can use very small voxel sizes (relative to the mesh), which after re-sizing would yield very large values of , and T^t. To avoid this, we chose l_r simply as the ratio of the I (isometric term) of the mesh at the smallest scale we used relative to the original mesh, divided by the λ of the smallest scale:

where λ_s is the smallest voxel size used for a particular cortex, I_s is the corresponding I for this cortex at the smallest voxel size.

Finally, I_o is the I term for the original cortical mesh. Indeed, the ratio of is always close to, but larger than one in our dataset. Thus, we can see in Fig. 2 A that most traces start very close to the original data point, indicating that our finest scale is reconstructing the original surfaces well.

Note that the re-scaling is isometric, meaning that is only affects I , but does not change the data in the K × S plane (Fig. 3). Rescaling by a factor proportional to λ² can therefore be understood as distributing the data points along the I axis, while l_r can be understood as fixing the position in the I axis relative to the I of the original mesh.

Fig. 2 A is produced with these rescaled quantities as described above. The only final step in producing Fig. 2 A, and in calculating the associated slopes, is the removal of artifactual data points where , which can occur at very large voxel sizes relative to the cortical mesh.

The algorithmic implementation in MATLAB can be found on Github: https://github.com/cnnp-lab/CorticalFoldingAnalysisTools/blob/master/Scales/fastEstimateScale.m, as part of our Cortical Folding Analysis MATLAB package https://github.com/cnnp-lab/CorticalFoldingAnalysisTools/, which also has been recently updated with a graphical user interface.

S3.2 Species-specific details

In the following, we will show the detailed data for each species in terms of their scaling behaviour in Fig. S3.2. Videos for all species, showing the pial surface at each scale can be found under https://bit.ly/3CDoqZQ for review purposes. Final versions of all underlying data, analysis code, and videos will be published on Zenodo and eLife upon acceptance of the paper.

Given the demonstrated overlap between species, and if their cortices are all approximations of the same form, how can one tell apart their cortices? The answer is that all approximations have a range of validity, which varies between cortices. For a gyrified cortex of area A_t [when no dependence on λ is indicated then we are referring to the values for the original cortex], the coarse-graining will remove details of ever-increasing scale, lowering A_t(λ) until attaining lissencephaly for A_e(λ_lys) = A_t(λ_lys). Indeed, it follows from rewriting Eqn. (3) in the form of Eqn. (2) that a given cortex’ shape will be comprised of self-similar structures with areas ranging from A₀ = g^-5A_t at their smallest to A_t at their largest, where g is the gyrification index . Whenever this self-similar scaling is valid, A₀ acquires a further interpretation as the typical size of the smallest structures in a cortex: patches smaller than that must be approximately smooth. Consequently, estimates the number of morphological features in each cortex (and adds a new interpretation for the gyrification index). For example, we estimate the human in our dataset to have about 105 morphological features in each cortical hemisphere, and the galago to have about 4 such morphological features. The corresponding N_structures(A) estimates how this number changes over coarse-graining. Suppl. S4 later provides a more detailed discussion on this topic.

Figure S3.2:

S4 Morphometric relations and the range of validity of the fractal approximation

There are many geometrical quantities that describe aspects of entire cortical hemispheres. But most information contained in such summary morphometric measures is captured, exactly or to good approximations, by the three used in Eqn (1): the pial (or total) surface area A_t, exposed surface area A_e, and average cortical thickness T. Other quantities, like the gyrification index , or the Grey Matter and Total volumes, can be expressed as products of power laws of these quantities; and the logarithm of products of power laws are linear combinations of the logarithms of the constituting variables, with their exponents as coefficients. Thus, the morphology of a given cortex can be fairly summarized as point in a 3-dimensional morphometric space with components given by log A_t, log A_e and log T² (the last exponent guarantees all axes have the same dimension of log[Area]).

This logarithmic structure guarantees that the quantities denoted by products of power laws of the original variables correspond to vectors in the morphometric space: the increase and decrease of their values occur along these vectors, and the planes perpendicular to each vector denote all configurations with the same value for the associated variable. For example, G = log g = log A_t - log A_e, the point in the morphometric space associated with the gyrification ratio value , defines a displacement along a vector with coefficients {1, -1, 0}.

The quantities K, S and I used in Sec. 3.2 likewise correspond to displacements along specific directions in morphometric space, which are all orthogonal to one another. These can be normalized, so that the sum of their squared coefficients equals unity. In this case, expressing cortical shape in the morphometric space in terms of the original variables or the new variables amounts simply to a rotation: a change of orthonormal base.

Figs. 2 and 3 can be then regarded as particular 2D snapshots, taken from a certain direction, of a 3D set of points.

Of particular note, those vectors with a sum of coefficients equal to zero correspond to dimensionless variables. Any such vector is perpendicular to the I direction, and lie on the K × S plane.

In this framework, we can picture the step-wise coarse-graining of a cortex, as described in Sec. 3.1, as a trajectory in this morphometric space. Empirically, Sec. 3.1, as codified by Eqn. 3, finds that the trajectories in morphometric space for all studied primates are linear, and largely overlap along a line of constant K. For cortical morphology, the first fact implies scale-invariance, the second implies universality.

Aging, development, and the progression of neurodegenerative conditions likewise each correspond to different morphometric trajectories.

Now consider a morphometric trajectory, corresponding to the coarse-graining of a cortex up to the point of lissencephaly, i.e., to the spatial scale λ_lys such that A_t(λ_lys) = A_e(λ_lys). In the idealized case (but very close to empirical, as seen in Sec. 3.1), all points along the trajectory should follow Eqn (3), starting at the original values for the morphometric measures and ending at the intersection of the K = constant coarse-graining trajectory and the g = 1 line that defines the limit for lissencephaly (points with g < 1 are geometrically possible, but are largely avoided by the trajectories of actual cortices, although not by those of all objects; see S5.1). This trajectory spans a range of values for each morphometric measure: for dimensionless morphometric measures, such span is fully specified by Eqn (3). For all others, one needs also to specify the voxel rescaling. The scheme described in Sec.S3.1 simply guarantees that, for a constant , the fundamental area element (and thus also the average thickness T) is kept constant for all realisations of the coarse-graining process. We can regard these as approximations of the same original cortex, all drawn with the same resolution, but in ever smaller isometric sizes.

This choice of spatial scale then allows our identification of A₀ as the typical size of the smallest morphological features in a given cortex. Equivalently, it is also the area of the largest possible lissencephalic cortex for a given average cortical thickness, just at the cusp of the onset of gyren-cephaly.

This same choice also enable us to easily compute the expected span of the various morphometric measures over coarse-graining, assuming the universal scaling: each will range between its original value (expressed as a function of the original values for A_t, A_e and T), and its value at spatial scale λ_lys, where A_t(λ_lys) = A_e(λ_lys) = A₀: we simply replace A_t and A_e by A₀, and T by . The ratios between the initial and final values of all morphometric measures will be given as powers of g, imbuing the gyrification index with a new significance.

We list below the exspans for all such morphometric measures (S = log s, K = log k and I = log v_I)

S5 Validation with non-brain objects

We chose to apply our coarse-graining procedure to a range of objects to validate our algorithm, and we included both negative and positive controls.

As a first positive control, we used a simple box with finite thickness (i.e. we simulated a inner “white matter” surface and an outer “grey matter” surface, both being a cube, one positioned inside the other). A schematic is shown in Fig. S5.1. We know based on theoretical consideration that this box must be aligning exactly with the line, as its A_t must be the same as its A_e at every stage of coarse-graining, which also implies a “fractal dimension” of 1. This is exactly the case in our plot of K against S, as all the grey data points from the coarse-graining algorithm align on (thick black line).

Figure S5.1:

As negative controls, we included three non-brain objects, a bell pepper, two walnut halves, and a coarse outline of a bust/figurine known as Laurana. All three of these objects show a self-similar, or fractal, regime (corresponding to partially straight trajectories in K × S) in our coarse-graining procedure, but their trajectories in K × S space are not flat (Fig. S5.1), meaning their fractal dimensions are distinct from 2.5. Their trajectories also do not overlap with each other, indicating that they are fundamentally different shapes.

The human brain is included here as a reference and shows a clear flat trajectory (Fig. S5.1) as presented in the main text. Also worth noting, after the human trajectory intersects the line, it veers off to follow the line closely, indicating that the A_e and A_t remain the same with further steps of coarse-graining. Effectively, the human brain transitions to be lissencephalic convex structures once their fractal regime ends. For simplicity, we excluded these lissencephalic data points of extreme coarse-graining from the results in the main text.

S6 Validation with randomising grid for coarse-graining

To assess robustness of the coarse-graining algorithm, we ran 30 different realisations of the algorithm, but with a small random shift of the grid position relative to the surface meshes. The shift is chosen to be within the radius of λ. This allows subtle changes in the voxelisation at the boundary of grey matter and white matter. Interestingly, over five different human individuals, we observed subtle changes in the values S over the 30 realisations, and to a lesser degree K (Fig. S6.1). Additionally, the variation between individuals, especially in K at smaller spatial scales, is far greater than the variation introduced by the jittered realisations. These results suggest that the coarse-graining algorithm is extremely robust towards subtle changes at the grey or white matter boundary, especially far away from the lissencepahlic limit.

Figure S6.1:

S7 Ageing process

S7.1 All morphometric variables

In the main text, we used the ageing process as an example for scale-dependent biological process. For completeness, here we also show all the morphometric variables across scales in the 20 y.o. and 80 y.o. cohort in Fig. S7.1, and the corresponding effect sizes in Fig. S7.2. As expected, A_e shows very little difference between the two cohorts. I, S, and T show some effects at scales smaller than 4 mm, but the effect decreases monotonously for higher scales. K demonstrates a more complex scale-dependent effect: larger K in younger subjects at small scales of 0.25 mm (in agreement with previous native scale analyses [6]), smaller K in younger subjects at scales of approx. 2 mm, and a return to larger K in younger subjects at large scales of more than approx. 5 mm. Note, however, that despite these large effect sizes, the actual change of values of K is within the range of variation expected across species (seen in the main text). The range of variation of K is approx 0.04 here, and at least an order of magnitude smaller than the range of variation in S (approx 1.5).

S7.2 Confirmation in independent dataset

To confirm that our observed ageing effect was not driven by sample or data-specific properties, we also analysed an independent dataset (NKI) with the same methods. Here, we show the equivalent figures to Fig. S7.1 and Fig. S7.2 for the NKI data in Fig. S7.3 and Fig. S7.4. The same qualitative patterns can be seen in all plots, including the scale-specific effects in K and A_e at around 2 mm.

Figure S7.1:

Figure S7.2:

Figure S7.3:

Figure S7.4: