Establishing a Baseline for Human Cortical Folding Morphological Variables: A Multisite Study

Differences in the way human cerebral cortices fold have been correlated to health, disease, development, and aging. However, to obtain a deeper understanding of the mechanisms that generate such differences, it is useful to derive one's morphometric variables from the first principles. This study explores one such set of variables that arise naturally from a model for universal self-similar cortical folding that was validated on comparative neuroanatomical data. We aim to establish a baseline for these variables across the human lifespan using a heterogeneous compilation of cross-sectional datasets as the first step to extending the model to incorporate the time evolution of brain morphology. We extracted the morphological features from structural MRI of 3,650 subjects: 3,095 healthy controls (CTL) and 555 patients with Alzheimer's Disease (AD) from 9 datasets, which were harmonized with a straightforward procedure to reduce the uncertainty due to heterogeneous acquisition and processing. The unprecedented possibility of analyzing such a large number of subjects in this framework allowed us to compare CTL and AD subjects' lifespan trajectories, testing if AD is a form of accelerated aging at the brain structural level. After validating this baseline from development to aging, we estimate the variables' uncertainties and show that Alzheimer's Disease is similar to premature aging when measuring global and local degeneration. This new methodology may allow future studies to explore the structural transition between healthy and pathological aging and may be essential to generate data for the cortical folding process simulations.


INTRODUCTION
Mapping the human brain development and aging longitudinally from a unique dataset is barely impossible due to our extensive lifespan. Suppose one wants to understand the time evolution of the brain morphology through the whole lifespan. In this case, it is mandatory to combine multiple datasets and devise methods that allow a fair comparison among them. The recent advent of large heterogeneous datasets (to ensure the legit results across several populations) and data curation innovations allowed neuroscience to explore the brain's changes on an enormous scale, from functional data to structural studies. Here, we propose a novel methodology of combining structural MRI from different acquisition sites and equipment, acknowledging their heterogeneity and providing insights about cortical folding during development and aging, with a practical application of these results in Alzheimer's Disease, the most common dementia worldwide. This effort is essential to provide experimental information about the evolution span of the cortex morphology and build a cortical folding theory.
In the past years, the study of cortical folding in humans and other mammals was extensively promoted by the application of translational science, open-access databases, and data science tools which allowed the field to grow in multiple and sometimes convergent directions (Madan and Kensinger, 2016;Llinares-Benadero and Borrell, 2019;Mota et al., 2019;Essen, 2020). Due to its biological background, cortical folding measurements have been included in human brain structure analysis, discriminating disease from healthy controls (Cao et al., 2017;Lamballais et al., 2020), describing its correlation with cognition (Núñez et al., 2020), and relation to aging (Hogstrom et al., 2013;Madan, 2021). Despite the lack of consensus defining a unique theory that explains cortical folding on every scale. Mota and Herculano-Houzel (2015) proposed a cortical folding model respected by more than 50 mammals' brain hemispheres. It predicts a powerlaw relationship between cortical thickness (T), exposed (A E ), and total area (A T ) (Equation 1) where α is the self-similarity index, a universal constant, with a theoretical value of 1.25 and calculated as 1.305 (Mota and Herculano-Houzel, 2015) when considering a heterogeneous dataset of different species. It is an index of complexity, reflecting how much detail in a real fractal pattern changes with the scale at which it is measured. The constant k is dimensionless and associated with conserved viscoelastic properties of cortical matter. Following this equation, the usual geometrical variables used to describe the cortex (T, A T , and A E ) are not independent in this model. Wang et al. continued this study, validating the theory for different groups of humans (Wang et al., 2016) in specific regions of interest (Wang et al., 2019b), obtaining new evidence that the cortex is indeed a self-similar structure. Wang et al. (2021) showed that from this model one could derive a more natural set of nearly independent morphological variables, K, S, and I, that could be used to improve disease discrimination from regular and typical structural changes of the brain, such as aging.
In this framework, the morphology of each cortex is expressed as a point in a three-dimensional abstract space, with each coordinate component corresponding to the log-value of an independent morphometric variable (Wang et al., 2021). The use of log-values guarantees that linear combinations of the basis vectors correspond to power-law relations. Conversely, as long as any new variables are expressible as power-laws, different but equivalent sets of variables can, thus, be related to each other by a change of Cartesian coordinate bases. As a starting point, we use the log-values of the commonly morphometric variables, the total area log 10 A T , exposed area log 10 A E and average thickness log 10 T 2 . We then derive our new variables thusly: K = log 10 k (Equation 2) is a near-invariant quantity obtained by isolating k in Equation (1). S (Equation 3), also dimensionless, encapsulates the aspect of brain shape that changes more significantly across cortices, and I (Equation 4), represents brain isometric volume and carries the information about overall cortical size. Both are derived from planes perpendicular to K, noting that K remains unchanged when applying an isometric transformation on the brain (a direction in which all areas scale equally) resulting in I, and the perpendicular dimensionless vector S resultant from K × I 12 . Those three variables combined may help distinguish pathological events similar to age effects, such as AD. K = log 10 k = log 10 A T − 5 4 log 10 A E + 1 4 log 10 T 2 (2) S = 3 2 log 10 A T − 3 4 log 10 A E − 9 4 log 10 T 2 (3) I = log 10 A T + log 10 A E + log 10 T 2 (4) Baseline values (age-specific norm values) of those cortical morphological variables and their inherent uncertainties over the human lifespan were never defined due to methodological limitations of combining multicenter studies of brain structural images (Gronenschild et al., 2012;Fortin et al., 2018). However, both are crucial for applied studies. Uncertainties determine measuring limitations and allow a proper comparison between different cohorts. Baseline values allow clinical applications and, when necessary, simulations of healthy control data. Multi-center studies have become more common in the last few decades due to improved data sharing and the open science trend. However, combining MRI data from multiple experiments is not a trivial task. Part of these uncertainties is explained by random errors, a direct consequence of the natural variance within the species individuals, and the impossibility of replicating the exact same condition over multiple data acquisitions. Theoretically, these random errors are the same if one considers multiple experiments. Confounding components could be added when gathering MRI images with differences in acquisition parameters, acquisition equipment, and versions of the post-processing software (Dickerson et al., 2008;Gronenschild et al., 2012). Together, they add uncertainty to the data due to a systematic effect, ultimately reflecting a systematic difference in morphological variables calculated from different databases. In practice, the natural random fluctuations are convoluted with the uncertainty due to a systematic effect, enhancing the data spread.
Recently, a vast literature on primary morphological variables from MR images (Schnack et al., 2010;Pomponio et al., 2020; 1 One may understand this as similar to "Principal Component Analysis", a change of basis on the data representation, but using scaling considerations rather than covariance eigenvectors to define the principal directions. The new set of variables are combinations of the usual geometrical variables used to describe the cortex. The principal component K is given by the power-law relation, which contains the correlation between the variables. The variables S and I are orthonormal to K. 2 One can also start requiring a plane perpendicular to K, impose the dimensionless requirement, resulting in the dimensionless vector S, and find I from K × S. Bethlehem et al., 2022), and more complex ones (Fortin et al., 2018), explored the limitations of compiling images acquired with different scanners and protocols and its implication in statistical analysis results. The suggested solution for these limitations is a procedure called data harmonization, in which there is an "explicit removal of site-related effects in multi-site data" (Pomponio et al., 2020).
It is essential to add that one cannot, in general, distinguish the natural variance in healthy human anatomy (σ Natural ) and data acquisition random errors (σ Random ). Also, it is not possible to distinguish between the uncertainty due to a systematic effect of the acquisition and the processing. In other words, any measurement setup can include an additional value at its final measure that is repeated for all measurements. Here, the acquisition effects are due to differences in MR structural images acquisition protocols and their parameters, while processing errors are due to differences in the software version, pipeline parameters, and computer used for processing (Gronenschild et al., 2012). In Fortin et al. (2018) acquisition and processing effects are described as scanner and site effects 3 . However, in the case of the methodology of this manuscript, it is impossible to quantify the contribution of this systematic shift, from acquisition and processing, since we did not compare processing or acquisition parameters by controlling the other error sources. For this, one would need a diligent procedure to track all differences between acquisition and processing as well as multiple image of the same individual for each sample. This manuscript deals only with the first case, approximating the random errors due to acquisition and processing for all samples using three subsequent structural images of the subjects available only for one particular sample (refer to Section Supplementary Material Section 2). The total uncertainty associated with the variables (K, S, and I) and their components is then summarized in Equation (5), the error propagation formula with no correlation since we are treating errors from uncorrelated sources by definition.
At which σ X , is the total uncertainty for K, S, and I, σ Natural represents the so-called here Natural Fluctuation, σ Random , the random variance that can be evaluated from the repeated measures with the same subject at each Sample acquisition and processing steps, and σ Sample , for the uncertainty due to a systematic effect on each sample. In this manuscript, we developed a simplistic harmonization method that allows a multi-site analysis, determining K, S, and I baseline values, their uncertainties, and the ratio of changing rates through the years for cross-sectional images of healthy controls. Although this work focused is these novel morphological variables, we also provided the baseline of usual morphological variables such as cortical thickness (T) and Gyrification Index (GI). A robust estimation framework of a baseline time function is important to provide future clinical studies with enough information to compare brain structure trajectories in case of a pathological investigation. We then, as a clinical application and extension of de Moraes et al. (2022) and Wang et al. (2019b), verify if Alzheimer's Disease (AD) is similar to a premature accelerated aging brain in terms of the independent morphological variables, K, S, and I, by comparing their rates. Finally, we have determined the time evolution of the self-similarity dimension α proposed in the Mota and Herculano-Houzel (2015) model, interpreting all results within the framework of this theory.

MATERIALS AND METHODS
Participants and data included in this analysis were either acquired by the Instituto D'Or de Pesquisa e Ensino (IDOR), Rio de Janeiro, Brazil; used at de Moraes et al. (2022) Wang et al. (2016Wang et al. ( , 2019b, as the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu), the Human Connectome Project (HCP), Information eXtraction from Images (IXI), Nathan Kline Institute (NKI) and Open Access Series of Imaging Studies (OASIS) Wang et al. (2016Wang et al. ( , 2019b; or from open-access databases, as Ultra-high field adult lifespan (AHEAD) (Alkemade et al., 2020) and (Snoek et al., 2021). The total number of subjects is 3,095 healthy controls from 4 to 96 years old and 555 Alzheimer's Disease subjects from 56 to 92 years old. Datasets' demographics are summarized in Table 1. AHEAD and HCPr900 only include age range instead of the actual age due to local Ethics Committee rules. To overcome the lack of an exact number, we assumed the interval's mean age as a reference.
To investigate the uncertainty due to a systematic effect during image acquisition and variances in repeated measures performed in the same conditions, we included the first 50 subjects from the AOMIC ID1000 (Snoek et al., 2021), described in Table 1.
The structural images AHEAD, AOMIC ID1000, PIOP1, and PIOP2 were processed in FreeSurfer v6.0.0 (Fischl, 2012) without manual intervention at the surfaces (McCarthy et al., 2015). The original project from IDOR images is longitudinal. Those images in particular were processed with the longitudinal pipeline. The FreeSurfer localGI pipeline generates the external surface and calculates each vertex's local Gyrification Index (localGI) (Schaer et al., 2008). Values of Cortical Thickness, Total Area, Exposed Area, and Local Gyrification Index were extracted with Cortical Folding Analysis Tool (Wang et al., 2019a). We defined as ROI the whole hemisphere, frontal, temporal, occipital, and lateral lobes (based on FreeSurfer definition of lobes). The lobes' area measurements were corrected by their integrated Gaussian Curvature, removing the partition size effect and directly comparing lobes and hemisphere cortical folding, as described in Wang et al. (2019b).
Datasets' demographics, acquisition, and processing information are summarized in Supplementary Material.

Statistical Data Analysis
Multiple comparisons of means were made with ANOVA, correlations between cortical folding independent variables, and age estimated with Pearson r and post-hoc evaluations with Tukey multiple comparisons of means, including comparisons within Diagnostic groups. The statistical significance threshold was set at α = 0.05, and multiple corrections (Bonferroni) were applied when needed. The Healthy Control subjects are used as a reference in the necessary normalization and harmonization procedures. A standard linear regression was used to obtain the self-similarity index α. All analyses were done considering the whole cortical hemisphere and the lobes (refer to Supplementary Material), as most diseases imply local or non-uniform global structural damage. Linear Mixed Models (LMM), standard linear regression, and statistics derived from their results were analyzed with RStudio (v1.4.1717 and R v4.

Harmonization
The data harmonization procedure is based on the time evolution of the basic morphological variables Y = T, A T , or A E , modeled as an exponential decay with the age t age . Using Linear Mixed Models (LMM) in the log-linear scale, we are able to perform a joint linear fit for each ROI, maintaining the same angular coefficient a among all samples but different linear coefficient b j,ROI according to Equation (6). We then subtracted the sample shift, given by the linear coefficients b j,ROI ), directly obtained from the Sample Random Effects (log 10 (Variable) ∼ Age×ROI+ (1|Sample : ROI)) in Healthy Controls. From this procedure, we directly calculate the variables K, S, and I and determine for all morphological variables of interest (also T, A T , A E , GI = A t /A e ) their baseline and uncertainties. This harmonization procedure was validated as the most probable description of the data using a Bayesian approach. This simple harmonization procedure is equivalent to the approach presented in Bethlehem et al. (2022) for the special case of a purely linear function. log 10 Y j,ROI = at age + b j,ROI , for each ROI and each sample j (6)

Bayesian Model Selection
In order to show that the harmonization procedure used was the most suitable, one can compare it with other models to describe the data. For example, the angular coefficient could also be a random effect for the LMM, varying across the samples. Using Bayesian model selection (Gregory, 2005), one can show that the LMM used in this study is the most probable description of the data. Considering our dataset comprised by J samples where D (j) = {d 1(j) , d 2(j) , · · · , d N j (j) , } denotes the data of a morphological variable in the j-th sample. The degree of belief in a model M that describes a normally distributed data for each age is given by the posterior probability: The term µ(x i(j) , A (j)β ) is given by the model M proposed to describe the data, with the set {A (j)β } being the free parameters of the model defined within a range A (j)β . A prior probability P(A (j)β |I) is assigned for all free parameters reflecting the prior knowledge about their values. The basic model of this study assumes that the log of morphological variables follows µ = ax + b (j) , with the dispersion σ (j) for each sample. Uninformative priors were used for the free parameters. This integral can be easily done numerically using Python3 with PyMC3 package or similar. The advantage of this approach is that it naturally penalizes the insertion of unnecessary new parameters avoiding over-fitting. With multiple models M i , one defines the odds as O ij = p(M i |D) p(M j |D) . As i p(M i |D) = 1, the probability of the model M 1 is Frontiers in Neuroscience | www.frontiersin.org given by: The alternative models considered by this study were: (i) µ = ax + b (j) , with the same dispersion σ (j) ; (ii) µ = a (j) x + b (j) , with the dispersion σ (j) for each sample; (iii) µ = a (j) x + b (j) , with the same dispersion σ (j) for each sample; (iv) we also tested the simplest non-linear trend µ = cx 2 + ax + b (j) , with the dispersion σ (j) for each sample; we have found that the LMM used in the study is the most probable. The code implementing this approach is given by the authors.
It is important to note that all parameters calculated using the LMM approach can be obtained in this framework by a slight modification of Equation (8). Instead of marginalizing all the parameters, by leaving one out, we obtain P(A j1 |D) which gives directly the most probable value for the parameter A j1 and its credible region.

Multisite Harmonization
The analysis starts with data harmonization, a crucial step to combining information from multiple sites allowing a fair comparison among datasets. Figure 1 shows the result of the harmonization procedure and its reflection on the primary morphological variables. The LMMs are used to fit the data assuming a single angular coefficient for all samples and different linear coefficients due to the uncertainty caused by systematic effect. We considered the variables' trend with age linear based on the correlation between T, A T , and A E with Age for the Healthy Control group calculated from our data. There was a statistically significant correlation for T (Pearson r = −0.6949; R 2 = 0.48; p < 0.001), A T (Pearson r = −0.4711; R 2 = 0.22; p < 0.001), and A E (Pearson r = −40.2005; R 2 = 0.42; p < 0.001). The harmonization was validated as the most probable description of this data using the Bayesian approach.

Baseline, Rate Estimates, and Alzheimer's Disease Diagnostic Discrimination
An immediate consequence of the harmonization is defining a baseline value and its uncertainty for any morphological variable considered in this work. The LMM provides quantitative information about the trends of these variables through age. The estimated uncertainties values are summarized in Table 2, and summaries of linear mixed models are available in Supplementary Material. The so-called Natural Fluctuation (σ Natural ), representing human diversity, was considered the Residual Standard Deviation of the models. The approximation for the random component (σ Random ) of the experimental error derived from the multiple-image acquisition and processing obtained using only 50 subjects from the AOMIC ID1000 dataset is fully described in Supplementary Material. We considered the Random Effect of the Sample category (intercept) Standard Deviation as the uncertainty caused by a systematic effect component of the experimental error (σ Sample ). The σ X is the resultant uncertainty, composed of the previously described errors, and is the square root sum of each one's squared values. Thereby, σ K = 0.026 (4.9%), σ S = 0.14 (1.6%), and σ I = 0.091 (0.9%) with the corresponding fraction of their healthy controls intercept.
Considering Alzheimer's Disease, we used the baseline (Figure 2) to verify if the age trends of K, S, and I are different between the diseased and Healthy Control groups. We extracted the rates of change for each variable and diagnostic, summarized in Table 3 and displayed in Figure 3. The slopes were compared with post-hoc pairwise mean comparisons.
There is a significant difference between Control and AD slopes in K (pairwise comparison estimate −0.00090, p < 0.001), S (pairwise comparison estimate 0.0013, p = 0.004), and I (pairwise comparison estimate −0.0017, p < 0.001) as expected from Figure 3. For K, the AD curve is virtually flat as if it had reached a plateau, while for I, we see a constant reduction in brain isometric volume, with a reduced intercept. Comparing the mean values within AD and CTL, there is a significant difference for K (pairwise comparison estimate 0.042, p < 0.001, representing 7.9% of the CTL mean) value, S (pairwise comparison estimate −0.14, p < 0.001, 1.5%), and I (pairwise comparison estimate 0.11, p < 0.001, 1.0%).
In order to complete the full description of the cortical morphology within the framework established by the Mota and Herculano-Houzel's model, we separated the data from both Health controls and patients with AD per decade and fitted the self-similarity dimension α theorized to be 1.25. The result is shown in Figure 4.
As most diseases imply local or non-uniform global structural damage, we expanded the results of trajectories to the lobes. We included ROI as a fixed effect in the linear mixed model equation with harmonized variables: Variable ∼ Age × Diagnostic × ROI + (1|Sample : Diagnostic : ROI). Rates are summarized in Supplementary Material. We discriminated AD and healthy aging (CTL) K trends with pairwise comparisons at the Frontal (p adj < 0.0001), Occipital (p adj < 0.0001), Parietal (p adj < 0.0001), and Temporal lobes (p adj < 0.0001), as expected from de Moraes et al. (2022). Despite the rates being smaller for AD than CTL, the results mimic the whole brain with Alzheimer's Disease trajectories appearing as a plateau during aging. The Frontal and Parietal lobe shows the least folded pattern (K) in AD (within lobe comparison are described in Supplementary Material). Also, there are significant differences between S pairwise comparisons in CTL and AD at the Frontal (p adj = 0.00035) and Parietal lobes (p adj < 0.0001); and among I, between CTL and AD at the Frontal (p adj < 0.0001), Occipital (p adj = 0.0031), Parietal (p adj < 0.0001), and Temporal lobes (p adj < 0.0001)).

DISCUSSION
We analyzed a combination of datasets aiming to establish the K, S, and I baseline through the human lifespan. The suggested model included complexity at two levels: accounting  for heterogeneous acquisition and processing and heterogeneous brain structures across healthy and disease subjects. The hypothesis is that non-homogeneous methodology (inclusion criteria, acquisition, and processing) implies more significant uncertainty due to a systematic effect. Thereby, we suggested a harmonization procedure that allows the comparison of data from multisite reducing this kind of experimental error. This result is innovative in its simplicity and focuses on expanding our independent morphological variables knowledge. The main caveat of our procedure is the impossibility of accounting for all differences between the samples and using this to define a universal harmonization according to each acquisition protocol. Our method does not depend on one reference sample or a gold standard definition. The removal of the shift for each sample is carried out based on the available healthy control data, which can always be added to the analysis. In other words, our method does not account for a possible uncertainty due to a global systematic effect that could shift all samples equally. By looking at the uncertainty in cortical thickness estimation, our results are compatible with previous work. Our findings of 0.02 mm ( Table 2) of uncertainty in cortical thickness are in agreement with a variability less than 0.03 mm found in a test-retest analysis by Han et al. (2006). Moreover, they suggest that uncertainties of 0.15 mm across platforms of manufacturers and 0.17 mm across field strength, at this study, coupled with processing uncertainties and estimated in 0.14 mm. A most recent study by Frangou et al. (2020) suggests a Mean Inter-individual variation of 0.07 ± 0.06 mm for hemisphere's cortical thickness at all age ranges in the study compared to 0.10 mm from our findings. This compatibility suggests that the possible global systematic effect is close to 0; thus, the harmonization could be regarded as being close to the universal one.
The typical values for the K, S, and I independent variables through the whole life span and their related uncertainties were obtained by this study. These curves not only can be used to discriminate pathological from healthy aging in crosssectional harmonized data but also are a powerful tool to morphologically describe non-typical cortex. Estimating the uncertainties is crucial to constructing the baseline curves and was carefully developed by this study. We estimated three types of uncertainties with linear mixed models and repeated measures: the uncertainty due to a systematic effect (combination of multiple datasets), repeated measures, and the natural variation in the human species. The error from acquisition and processing were coupled at the individual and sample levels. We expect future studies to disassociate the uncertainties in detailed factors (test-retest analysis, across platforms, field strengths, and acquisition parameters with the same sample of subjects) to measure morphological variables of interest better.
We have demonstrated the usage of the baseline curves by analyzing AD data. In terms of K, a measure of axonal tension, the AD curve is virtually flat as if it had reached a plateau of reduced brain tension. Compared with Healthy Controls, this plateau seems to be the same experienced by the oldest subjects. Interpreting in terms of the Mota and Herculano-Houzel cortex model, AD can be regarded as a form of accelerated aging from the axonal tension perspective. On the other hand, the evolution of S, I, and even the self-similarity dimension α suggest that the AD cortex is not morphologically similar to the health one, with a very different shape, isometric volume, and time evolution. However, it is interesting to note that both S and I converge to a value compatible with the oldest health subjects. In the future, correlating these results with the elasto-mechanic properties of the AD cortex will help to understand better the physical limits of the Mota and Herculano-Houzel model and expand it theoretically.
A significant result that could only be achieved by combining multiple data is the value of the self-similarity dimension α through the human lifespan and its evolution through atypical conditions. We extended the analysis in de Moraes et al. (2021) and verified how the slope α behaves with healthy aging. We found that the slope is compatible with the theoretical prediction for subjects between 20 and 60 years old. Defining this range of applicability of the theory is essential to understanding the period of life in which the model's basic assumptions are still valid. Thereafter, we confirmed the findings that the slope escapes the model after 60 years old, diverging from the theoretical value of 1.25. Before 15 years old, we have found a hint that the slope is also inferior to 1.25, in disagreement with the model. In addition, a non-linear pattern is suggested by the K baseline function before 20 years. One of the leading hypotheses to explain the deviation of the measured α from the expected value in both cases (before 15 and after 60 years old) is the breakdown of the homogeneity on the cortical surface. The suggestion of a homogeneity breakdown makes it particularly interesting to further study AD in the age range compatible with 1.25. Understanding the regimes outside the model is a fundamental step toward the theoretical understanding of cortical folding. Coupled with analysis of neuropathologies such as AD, these empirical results improve the cortical folding theory.
More than discriminating between disease and healthy subjects, understanding the pathology trajectory is vital to defining landmarks, indicating future injury, and collaborating to develop efficient treatment. Notably, in Alzheimer's Disease and its prodromal form, Mild Cognitive Impairment, subsequent    (B) For S, the AD and CTL patterns have similar intercepts, but statistically different slopes (AD slope p = 0.004 and CTL slope p < 0.0001; pairwise comparison estimate 0.0013, p = 0.004). (C) For I, that reflects brain volume, CTL has a decreasing volume with aging, while AD has a smaller slope (AD slope p < 0.0001 and CTL slope p < 0.0001; pairwise comparison estimate −0.0017, p < 0.001). The curves for the AD group were extended to early ages to ease the comparison.
phases of damage may start at least 20 years before the diagnosis (Frozza et al., 2018). The structural damage is an initial state than cognitive alterations, such as memory and clinical functions (Jack et al., 2010b). As pointed out by Fjell et al. (2014), healthy aging also contributes to reductions in cortical thickness and understanding the overlap between healthy aging and pathological aging triggers is crucial to segregating both events. Nevertheless, the possible similarity of structural damage in Alzheimer's Disease and aging is explored by multiple previous studies (Pacheco et al., 2015;Gutierrez Becker et al., 2018;Huizinga et al., 2018;Wang et al., 2021). We hypothesize that the structural changes due to AD are reflected in measures of the brain's global structure and shape. Therefore, cortical folding derived variables K, S, and I are sensitive to this perturbation. We verified this hypothesis by analyzing the difference in aging trajectories between healthy and pathological brains. At a more refined grain scale, the lobes, the results are compelling to confirm (Wang et al., 2021) suggestion that K, S, and I together would be a powerful tool to discriminate similarly but distinct events, such as AD and aging. In this case, the difference between events is in the cortical structure's changing velocity or rate. Then with the rates for K, S, and I at the lobes, we could explore whether a lobe would unfold faster. For AD, the temporal lobe has the biggest unfolding (K) and shrinking rate (I), while healthy aging is more aggressive in the Parietal lobe, followed by the temporal and frontal lobes.
The natural expansion of this study is to include Mild Cognitive Impairment (MCI) subjects. As seen before, the MCI is an intermediary step between healthy aging and AD in the disease progression, cognitive impairment, and structural changes (Jack et al., 2010b;Frozza et al., 2018;de Moraes et al., 2022). Also, we intend to continue this study by applying this approach to a longitudinal dataset composed of MRI structural images with MCI subjects that eventually convert to AD and compare to the cortical folding components of non-converters. Including this data would allow having a more reasonable aging trajectory by having the transition period from healthy aging to pathological aging. Furthermore, we expect to explore K, S, and I function with age for the age range comprising newborns to 25 years old in future studies. This study will help better understand the value of the slope α at this age, revealing if there is an inflection point that can be connected to the development/aging transition in human brains (Davis, 2021), such as suggested by the data.
This study successfully achieved its goal of combining multiple samples with a simplistic harmonization procedure that adjusts K, S, and I values for the 3,650 subjects. Hereafter we estimated K, S, and I typical values and aging rates to discriminate against Alzheimer's Disease subjects based on their changing rates. To validate the proposed methodology, we analyzed results for Cortical Thickness, which are in concordance with previously published findings. At a glance, the significant differences found in aging rates for the cross-sectional data are suggestive that a brain with Alzheimer's Disease has premature and accelerated aging in terms of morphology.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. Data from previous publications are available from their repositories: Wang et al. (2016Wang et al. ( , 2019a. Data from AHEAD (Alkemade et al., 2020) and AOMIC (Snoek et al., 2021) datasets with the cortical folding variables extracted by these authors are available at https://doi.org/10.5281/ zenodo.5750619 (de Moraes et al., 2021). Data acquired at IDOR does not have clearance for public sharing patients' information that could lead to identification due to local Ethics Committee approval restrictions but will be shared upon request.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Copa D'Or Research Ethics Committee under the protocol number CAAE 47163715.0.0000.5249. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
FM, VM, FT-M, and BM contributed to the conception and design of the study. FM organized the database. FM and VM performed the statistical analysis and wrote the first draft of the manuscript. FT-M and BM supervised and were responsible for the project administration and funding acquisition. All authors contributed to the article and approved the submitted version.

FUNDING
This study was funded by the Research Support Foundation of the State of Rio de Janeiro (FAPERJ), National Council for Scientific and Technological Development (CNPq), and intramural grants from Instituto D'Or de Pesquisa e Ensino (IDOR). BM was supported by Instituto Serrapilheira (Grant Serra-1709-16981) and CNPq (PQ 2017(PQ 312837/2017. ADNI data included in this work collection and sharing was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering,