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, and average cortical thickness T:

This relationship 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, elegant and hierarchical way of expressing cortical shape. Specifically, we introduce a coarse-graining procedure that will render cortical shapes at different spatial scales, or resolutions. We show that the shapes of primate cerebral cortices at each spatial scale can be understood as approximations of the same universal self-similar archetypal form.

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

2 Coarse-graining procedure

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 and lissencephaly 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¹: At and Ae 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 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 is 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].

Here, for the first time, we show how to directly construct morphologically plausible realisations of cortices at any specified spatial scale, or resolution. This is achieved through a coarse-graining procedure 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 could be fused. This method is a new systematic way of expressing cortical shape, in terms of a flow of morphometric measurements as spatial scale varies. By examining how areas scale across coarse-grained realisations of actual primate cortices, we will be able to directly verify cortical self-similarity.

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 (see Suppl. S1). Algorithmically, we then segment the space between the original pial and white matter surfaces into a 3D grid of cubes of the desired scale λ, where each cube is of dimensions λ × λ × λ. We also term the 3D grid of cubes “voxelisation”, as it effectively captures the shape of 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” (Fig. 1, and videos: https://bit.ly/3CDoqZQ). A detailed description and discussion of the algorithm in provided in Suppl. S2.

Finally, 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 accordingly, which ensures that we can compare the coarse-grained realisations to the original cortices, and test if the former, like the latter, also scale as in Eqn. (1).

Figure 1

3 Results

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

We have analysed cortices of 11 different primate species (see Suppl. S1): 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.

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.

Figure 2

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 previously introduced [24] set of independent morphometric measures, K, I and S, that summarize the shape of the cortex in a natural and statistically robust way. In this framework, isometrically scaled copies of the same shape 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 summarizes their shape. 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 [24] 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.

One can directly compute the values of K and S for any cortex-like object, thus ascribing to it a point in the K × S plane. But a fuller expression of a particular shape is captured by the trajectory of said object as a function of coarse-graining. This is a very convenient and informative way of summarising shape: Self-similar objects correspond to straight trajectories. 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; and to horizontal trajectories (K = constant) if d_f = 2.5 (Fig. 3 B). A shape is said to be universal for a set of objects when their trajectories overlap, so that they can all be regarded as coarse-grained versions of one another (Fig. 3 C).

From this perspective, the regularity observed in primate cortices is remarkable (Fig. 3 D). The value of K is nearly invariant in all cases, in parallel to what was found in [4]. 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 cortical shapes are morphometrically similar to and in terms of the universal law ‘as valid as’ actual existing mammalian cortices. In contrast, we tested various non-brain shapes, and while e.g. the walnut, and bell pepper shapes form (partially) straight lines, they vary in both K and S (see Suppl. S5). These shapes may have a fractal regime, but their fractal dimension is not 2.5, nor is their shape 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 shape. In this framework, our main result can thus be expressed simply: For all cortices we analyzed, and for none of the non-cortices, coarse-graining will leave K largely unaffected, while morphological complexity S will decrease.

Figure 3

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 K(λ) as a function of scale λ for a young (20 year old) vs. an old (80 year old) group of human brains. Strikingly, the difference in K between the groups takes a U shape, where the strongest positive effect is seen at sub-millimeter scales, but the effect reverses and turns negative at millimeter scales (older subjects have higher K), peaking at about 2 mm. For scales over 5 mm the differences become negligible, suggesting the aging process has little effect on the larger cortical morphological features. Finally, we reproduced these results in an independent dataset in Suppl. S7.2.

In this particular example, the dramatic scale-dependency of the effect of K can be visually and intuitively understood by looking at changes in total surface area (Fig. 4 B, also see Suppl. S7.1): Fig. 4 C 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 (see e.g. [25, 26] for recent investigations and references therein) widening between gyral walls and decrease in gyral surface area at the crown. 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 very similar near-lissencephalic shapes.

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.

Figure 4

4 Discussion

We have devised a new way of expressing the shape 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 shapes are all approximations of 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.

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 sulci and gyri. We propose this new syntax as the basis for a more rigorous characterization 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 [24] neuropathic conditions, especially if combined with a regionalized version of this method applied to specific cortical regions [7] or local patches [8].

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 shapes 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, and one that moreover is applicable to individual cortices.

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, which is both a strength and a weakness. There are many complex phenomena in nature that display simple and universal scaling that can be derived from first principles [27, 28, 29]. In addition, recent results suggest simplicity and symmetry are generically favoured on statistical-ensemble grounds by evolution [30]. Our model correctly predicts the scaling law (Eqn. 1), but a more complete explanation for cortical gyrification is probably far more complex [31] 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, 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.

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 neuron numbers [4], number of functional areas [32] or the number of cortical columns [33], 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 [34, 35], and more broadly it is conserved over kinship [36] and phylogeny [37, 38]. It thus seems likely that comparative neuroanatomical methods [39, 40] 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 structures.

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 characterizations 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 analyzing cortical morphology, besides revealing a hitherto hidden regularity of nature, can become a powerful tool to characterize and compare cortices of different species and individuals, across development and aging, 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ção Serrapilheira Institute (grant Serra-1709-16981) and CNPq (PQ 2017 312837/2017-8).

Supplementary

S1 Primate species and their cortical surface reconstruction

S1.1 Human data

To study healthy human adults, we used the Human Connectome Project (HCP) MRI data, available at https://db.humanconnectome.org/ [41], obtained using a 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/ [42, 43]). Cam-CAN used a 3T Siemens TIM Trio System with 1 mm isotropic voxel size (for more details see [42, 43]). From the Cam-CAN dataset we retained 644 subjects that successfully completed preprocessing (with Freesurfer recon-all) without errors. From these subjects we randomly selected 5 females and 5 males between the ages of 19-21 (forming the 20 y.o. cohort); we also randomly selected 5 females and 5 males between the ages of 79-81 (forming the 80 y.o. cohort).

To confirm the ageing results, we also obtained an independent dataset from the Nathan Kline Institute (NKI)/Rockland sample [44] (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 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.

S1.2 Non-human primate data

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 N4BiasFieldCorrection, and Human Connectome BiasFieldCorrection 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.

The marmoset MRI structural scan was collected as part of the development of the NIH marmoset brain atlas [45, 46]. 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 [47], decimation [48], and remeshing. Reconstructed pial and white matter surfaces were visually quality controlled in conjunction with the MR images.

The remaining NHP MRI and subsequent brain surface extraction is detailed in [49, 50], 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.

S1.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.

S2 Coarse-graining algorithm

The core coarse-graining algorithm underpinning our analyses takes pial and white matter surfaces from the cortical surface reconstruction (Suppl. S1) 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. We then assign voxels with all eight corners inside the original white matter surface to the white matter voxelization. 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 gryi 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 λ, 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 coarsegrained 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 shape 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/.

Note that this process is substantially different from existing approaches measuring the fractal dimension of cortex. We do not simply apply successive convolution with an increasing kernel size (or equivalent), which effectively would achieve a uniform and spatially isometric “smearing” of the original cortical shape 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 [51]. Of course, our proposed procedure is not the only conceivable way to erase shape details below a given scale. However, it requires no fine-tuning, is computationally feasible and conceptually simple, thus making it a natural choice.

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 preposition in future. For example, an alternative method inspired by [52] 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. [53]) can be regarded as a type of reverse melting, and could perhaps be implemented and tested in a similar fashion as fine-graining procedures.

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.

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.

Figure S3.1:

λ = 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. I.e. we can use very small voxel sizes (relative to the mesh), which after re-sizing would yield very large values of . 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 Zenodo: TBA

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 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 ², the coarsegraining 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(λ) 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 shapes; 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 gyrencephaly.

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 At 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 shapes

We chose to apply our coarse-graining procedure to a range of shapes 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 shapes, a bell pepper, two walnut halves, and a coarse outline of a bust/figurine known as Laurana. All three of these shapes 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 chimpanzee 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 chimpanzee 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, these brains transition 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.

Figure S7.1:

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 (see Section S1). 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.2:

Figure S7.3:

Figure S7.4: