DNA methylation profiles were used in this work, more specifically, 22905 cpg sites from 6 GSE profiles, the details of which are shown in Supplementary Text 1, 2. The aging predictor and disease predictor were modeled by the classification ensemble algorithm of the classification tree (shown in detail in the section “Materials and Methods” below), and the best predictor was selected based on 10-fold cross-validation (Table 1 and Supplementary Text 1). The learning curves of the training data set are shown in Figure 2.

According to the cross-validation results, the top 46, 35, 14, 23, and 34 dimensions of cpg sites were identified as risk biomarkers related to aging, AD, PD, PSP, and FTD, respectively. Next, these five prediction results were verified in the independent test data set. This ensured that the classification results in the test data had adequate accuracy (Table 1). The receiver operating characteristic (ROC) curves of the predictors are also shown in Figure 2. The area under the ROC curve (AUC) 0.818, 0.693, 0.642, 0.903, and 0.857 in the predictor of aging for AD, PD, PSP, and FTD, respectively. Overall these results showed that the predictor had high accuracy and efficiency.

It has been reported that the identified biomarkers are closely related to aging and ND (Table 2; Salmerón et al., 2001; Hu et al., 2004; Katoh, 2008; Qiu and Ghosh, 2008; Fujimoto et al., 2009; Su et al., 2009; Kozarova et al., 2011; Euskirchen et al., 2012; Berwick and Harvey, 2014; Lopez-Pelaez et al., 2014; Alquézar et al., 2016; Chailangkarn et al., 2016; Densham et al., 2016; Husain et al., 2016; Pîrşcoveanu et al., 2017; Kurihara et al., 2019; Monteiro and Forti, 2019). Interestingly, these related functions (i.e., DNA damage repair response, mitochondrial dysfunction, or Ca ion homeostasis disorder), further indicated the imbalance of cellular homeostasis, and abnormal cell death was further induced across various ND types. Therefore, cellular homeostasis disorder may not only disrupt normal brain functioning but also lead to pathological changes related to ND.

To study the accelerated aging pattern in ND, the aging score was calculated based on 46 aging risk biomarkers (shown in the section “Materials and Methods” below). Both the median and median scores for different age groups are shown in Table 3. With the increase in age, the aging scores in both disease and normal aged individuals showed an upward trend. The results also showed that the accelerated aging pattern in ND compared with the normal aged people, which was consistent with previous studies (Liu et al., 2020). Compared with the chronological age, the aging score was more informative. The Kolmogorov–Smirnov (K–S) test was used to test the aging scores of each ND type as well as the healthy aged sample coming from the normal distribution. The results showed that the p-value was close to 0 and rejects the original hypothesis (shown in Supplementary Table 1). Consequently, the Kruskal–Wallis test was then executed to verify whether the aging score reflected a significant difference between ND and normal individuals, where different age groups, as well as different ND types, were compared. The results are shown in Figure 3, Table 3, Supplementary Figure 1, and Supplementary Table 2. These results indicate that the aging scores of ND sample individuals exhibit a significantly accelerated aging pattern, compared with those of normal aged individuals (p < 0.05).

To better study the potential mechanisms between aging and ND, the aging acceleration differential network was constructed, where FDR < 0.1 was used to ensure reliable correlation. The (partial) coefficients were also summarized and used to compare the relationship of each pair of cpg cites in the context of aging. Aimed at verifying the scale-free characteristics, the probability corresponding to logarithmic transformation and its degree were used to test the power-law distribution (shown in the section “Materials and Methods”). The node degree distribution curve and Pearson correlation coefficient of different ND types are shown in Figure 4.

The results revealed that the aging acceleration network was with the scale-free pattern. The degree and frequency were inversely correlated, and only a small proportion of genes had a high degree. Furthermore, the Fisher’s exact test was carried out to calculate the similarity between the different networks based on the training data and the test data, respectively, and the result showed a p-value very close to 0 (1e-23238, 1e-131120, 1e-3103, 1e-44540, and 1e-139070 for AD, PD, PSP, FTD, and all NDs, respectively). In addition, the node with the highest degree in the network was also significantly correlated with maintaining cellular homeostasis (Table 4; Crowe et al., 1994; Dong et al., 1997; Keats et al., 2007; Tang et al., 2011; Kaneko et al., 2016; Kula et al., 2019; Anerillas et al., 2020).

To further investigate the potential mechanisms between aging and ND, each shortest path of aging-AD, aging-PD, aging-PSP, and aging-FTD were identified based on the aging acceleration differential network using the Dijkstra algorithm. Enrichment analysis was then performed based on each of the shortest paths, the special functions (i.e., the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and the Biological Process (BP) term in Gene Ontology (GO), shown in Figure 5 and Table 5). These indicated multiple dysfunctions of different ND types (McGeer et al., 1989; Behl, 2000; Nagatsu et al., 2000; Everse et al., 2011; Obulesu and Lakshmi, 2014; Caputi and Giron, 2018; Cui and Xu, 2018; Yan et al., 2018; Liu et al., 2019; Porro et al., 2019; Starr, 2019; Ashrafizadeh et al., 2020; Burgaletto et al., 2020; Luo et al., 2020; Paul et al., 2021), respectively.

In addition, common enriched functions were investigated to study the crucial mechanisms across different ND types, where the top ten KEGG pathways are shown in Figure 6. For example, the Toll-like receiver (TLR) signaling pathway was the top KEGG pathway related to all four ND types. The TLRs activated their downstream pathways and then induced NF-κB and Pro-Il-1β, both of which are related to neuroinflammation and the pathogenesis of a variety of neurological diseases (Azam et al., 2019). Furthermore, the activation of TLRs can also induce the inflammatory response of macrophages by activating the transcription cascade (Lauterbach et al., 2019). The inflammatory response induced by activation of TLRs is considered to be closely related to AD, PD, PSP, and FTD (López González et al., 2016). Therefore, the importance of the Toll-like receiver signaling pathway is reflected in a variety of ND types. Furthermore, cytokine receptor interaction is also critical in ND progression. For example, microglia can release cytokines to trigger a cellular inflammatory response and participate in the occurrence and development of NDs (Colonna and Butovsky, 2017). Cellular apoptosis is also a fundamental process in the progression of nervous system diseases (Radi et al., 2014).

In order to further investigate the complex cellular homeostasis imbalances across different ND types, a global sensitivity analysis was performed using the MCMC method (shown in section “Materials and Methods”). To study the relationship between aging and ND markers, the differential MCMC K–S statistics were calculated based on each “aging-ND” pair; then the top ranked “aging-ND” pairs were identified for each disease (as described in the section “Materials and Methods”). The aging markers shared by different ND types as well as the disease markers within each specific ND type were then explored, based on both aging and ND MCMC differences.

TSC2 was the risk marker with the highest absolute MCMC difference across the four types of ND. It has been reported that the phosphorylation of TSC2 could restore the expression of NF-κB protein, which is vital to reduce the production of inflammatory mediators and ND markers (Kumar et al., 2017). In addition, the aging markers identified in all four diseases were identified according to the obtained “aging-ND” pairs. Among them, there were two aging markers with high frequencies. CHGN-1 (CSGALNACT1, frequency = 40) is a key enzyme for the production of CSPGs, participating in the demyelination, remyelination, axonal degeneration, and regeneration of the central nervous system (Saigoh et al., 2016). SLC15A2 (PEPT2, frequency = 24) can affect the production of pro-inflammatory cytokines by macrophages (Gupta et al., 2013). The excessive release of pro-inflammatory cytokines causes neuronal damage (Shen et al., 2004).

In AD, CDCA7L was with the highest absolute MCMC difference. It was associated with the intelligence quotient (Pan et al., 2011), which could even inhibit the Monoamine oxidase A(MAOA) promoter as well as neuronal cell death (Ou et al., 2006).

EFNB2 had the highest absolute MCMC difference in PD. EFNB2 was involved in adjusting the development of the nervous system and the neuronal migration (Lévy et al., 2018), and could also activate the Eph/efn forward signaling pathway as well as cell apoptosis (Zhong et al., 2019).

In PSP, the marker with the highest absolute MCMC difference was BCL2, which could affect the activity of the mitochondrial complex (Chong et al., 2020). BCL2 could also enhance the anti-apoptotic effect of nerve growth factor (NGF) and promote the survival and differentiation of nerve cells (Troullinaki et al., 2019).

In FTD, MICAL1 had the highest absolute MCMC difference. Through interacting with STK38 and STK38L, it acted as a negative regulator of cell apoptosis (Zhou et al., 2011). MICAL1 is also involved in the regulation of lamina specific connections in the nervous system (Schmidt et al., 2008).

The potential network markers were also found out based on each “aging-ND” pair from the results of sensitivity analysis by summarizing the betweenness in the aging acceleration differential network, as shown in Table 6, Figure 7, and Supplementary Figure 2. As a result, DUSP12 had the largest betweenness across the four ND types. DUSP12 could regulate the c-Jun N-terminal kinase (JNK) signaling pathway by dephosphorylating its substrate, which was critical to cell differentiation, apoptosis, and other neural functions in ND progression (Ha et al., 2019). It also has been reported that the overexpression of DUSP12 inhibited the production of pro-inflammatory cytokines and chemokines (Cho et al., 2017).

Furthermore, the genes (cpg site) within each shortest “aging-ND” path were also identified. For example, the top network (betweenness) marker in AD was RALGPS2. It has hitherto been revealed that RALGPS2 silencing induces cell apoptosis, by improving cell cycle inhibitors p27 and p21 (Santos et al., 2016). The top network marker in PD was MEOX2 (GAX). The deletion of MEOX2 can lead to decreased capillary density as well as resting cerebral blood flow, and then promote the loss of hypoxic angiogenesis in the brain. It causes hypoxia reaction and cell death (Wu et al., 2005). QRSL1, whose mutation may lead to dysfunction of mitochondrial energy production and mitochondrial disorder (Albers and Beal, 2002), was the top network marker in PSP. The top network marker in FTD was C3. C3 uptake by cells reduces stress-related cell death (e.g., oxidative stress or starvation). In the process of inflammation, storing C3 in cells can prevent some substances that induce cell death (Kulkarni et al., 2019). In brief, these network markers were also informative and indicate the relationship between cellular homeostasis and NDs.

The DNA methylation profiles were obtained from the Gene Expression Omnibus (GEO) database (Supplementary Table 1), including GSE15745, GSE51923, GSE53740, GSE57361, GSE66351, and GSE138597, along with the age index. These datasets were from seven different platforms: GPL6104, GPL8178, GPL8490, GPL5175, GPL13534, GPL 11154, and GPL 21145.

The steps of obtaining DNA methylation profiles were as follows:

The final sample included 366 healthy youth samples (age ≤ 50, 250, in the training data and 116 in the test data), 442 normal old individuals (age > 50, 300+142), 128 AD individuals (85+43), 36 PD individuals (25+11), 123 FTD individuals (85+38), and 42 PSP individuals (30+12). The DNA methylation data set included 22905 cpg sites (Supplementary Text 2).

After randomization as well as a random disorder, the healthy population samples were divided into training data set and test data set. The ratio of training data set samples to test data set samples were close to 2:1. The ReliefF algorithm was used to select key features, then the first 50 models were studied to train predictors. The optimal model was selected by 10-cross validation. To verify the accuracy of the aging predictor, the selected model was verified in the test data set.

(1) The normal aged group (age > 50) were labeled as 1 and the young healthy group (age ≤ 50) were labeled as 0.

(2) The 22,905 cpg sites were sorted by the ReliefF algorithm.

(3) The predictor was generated using the ensemble learning algorithm.

The ensemble learning algorithm established 100 decision tree models, and then put the data in the 100 classifiers for decision-making. The classification result was optimized to get the most answers from the 100 classifiers. For 100 classifiers, weak weight was given to the classifiers with more wrong classification results using the formula of the weight coefficient. The calculation method was as follows:

and the classifier weight was:

where m was the number of predicted variables; D_n was the weight of each sample;

ε_n was the classifier error; x_n was the member of predicted variables; y_n represented the corresponding attribute value of x_n and h_n was each attribute.

One hundred weak classifiers were superimposed to generate the strongest classifier. The optimal model was selected by 10-fold cross validation. Ultimately, the model with the highest accuracy rate was chosen. The identified features were considered as aging and each ND marker, respectively.

To construct disease predictors, 442 healthy aged individuals, 128 AD individuals, 36 PD individuals, 123 FTD individuals, and 42 PSP individuals were randomly chosen. In the classification process of each disease predictor, the ND sample was labeled as 1 and the normal aged sample was labeled as 0. To avoid the unbalance of samples in machine learning, the normal aged group (training data) were divided into 3, 9, 7, and 3 subgroups to compare with AD, PD, PSP, and FTD (Supplementary Table 2), respectively.

For each sample, the aging score was calculated as follows:

The transformed age was predicted based on the aging scores using linear regression:

where b is the regression coefficient for aging score.

To further reveal the relationship between aging and ND, the aging acceleration network was constructed based on the training data set and test data set.

As a result, five types of aging differential networks were constructed: AD, PD, PSP, FTD, and all four NDs together.

Global sensitivity analysis was used to investigate different cellular homeostasis imbalances in the background of aging acceleration, based on the MCMC method. Both the common and specific characteristics across different ND types were analyzed. To explore the common characteristics among different ND, we made a presumption that if the marker indicated the high risk scores in more than one ND type, then it might reveal the common mechanisms across multiple ND types.

The MCMC method was used for sampling from certain posterior distributions following a given probabilistic background in a high-dimensional space. The key step in MCMC was to construct a Markov chain whose equilibrium distribution equals the target probability distribution. It proceeded as follows:

where P(θⁿ| X) and P(θ^∗| X) are the posterior probabilities of the n-th accepted sample and the new sample, q(θⁿ→θ^∗) represents the transition probability from the n th accepted sample to the new sample, and q(θ^∗→θⁿ) is the transition probability from the new sample to the n th accepted sample.

In this article, the mean value of the aging score as well as four types of ND scores were used to evaluate the common characteristics across multiple types of ND. The mean value of the aging score and each ND risk score was used to evaluate the special characteristics of each type of ND.

where F1 was the cumulative distribution of samples with a minus value after normalization, whereas F2 was the cumulative distribution of samples with the plus value. The interval for the K–S statistic was set to two (based on the sign of normalized DNA methylation profiles). The MCMC pseudocode is provided in Supplementary Text 4, 5.

The biological functions of the existing genes were found by enrichment analysis. Gene Ontology (GO) terminology and KEGG pathway were taken from the gene set enrichment analysis (GSEA) platform (version 7.2)¹. A hypergeometric test was used to estimate the enrichment degree of the KEGG pathway or go BP term. The hypergeometric test formula was:

where N was the gene set of the whole gene, M was the known genes (e.g., KEGG pathway, or BP terms), N was the number of identified genes in each shortest pathway, and k was the number of common genes between the known genes and the identified candidate genes (in each aging-ND shortest path). The p-values of each pathway were controlled by the Benjamin-Hochberg (BH) method. To ensure the reliability of the results, FDR < 0.05 was selected.