Abstract
Age-related changes in the BOLD response could reflect neurovascular coupling modifications rather than simply impairments in neural functioning. In this study, we propose the use of a sparse dynamic causal model (sDCM) to decouple neuronal and vascular factors in the BOLD signal, with the aim of characterizing the whole-brain spatial pattern of hemodynamic sensitivity to healthy aging, as well as to test the role of hemodynamic features as independent predictors in an age-classification model. sDCM was applied to the resting-state functional magnetic resonance imaging data of a cohort of 126 healthy individuals in a wide age range (31 females), providing reliable estimates of the hemodynamic response function (HRF) for each subject and each region of interest. Then, some features characterizing each HRF curve were extracted and used to fit a multivariate logistic regression model predicting the age class of each individual. Ultimately, we tested the final predictive model on an independent dataset of 338 healthy subjects (173 females) selected from the Human Connectome Project in Aging and Development cohorts. Our results entail the spatial heterogeneity of the age effects on the hemodynamic component, since its impact resulted to be strongly region- and population-specific, discouraging any space-invariant–corrective procedures that attempt to correct for vascular factors when carrying out functional studies involving groups with different ages. Moreover, we demonstrated that a strong interaction exists between certain right-hemisphere hemodynamic features and age, further supporting the essential role of the hemodynamic factor as independent predictor of biological aging, rather than a simple confounding variable.
Significance Statement
By inferring region-wise hemodynamic profiles at the individual level, this is the first study providing an exhaustive whole-brain characterization of the hemodynamic sensitivity to healthy aging, reporting further evidence of the vascular changes across the adult lifespan. Using a predictive framework, we analyzed the statistical influence of advancing age on individual regional hemodynamic attributes, offering a quantitative evaluation of the diverse hemodynamic bias across different brain regions. We then unveiled a specific set of hemodynamic predictors to discriminate young from elderly people, mainly describing vascular properties of right-hemispheric areas. This suggests the asymmetric nature of vascular degeneration processes affecting the human brain at the latest stage of life, other than a potential biomarker that could be relevant for brain-age prediction.
Introduction
BOLD functional magnetic resonance imaging (fMRI) is widely used to study neural functioning, often overlooking its indirect nature that arises from the coupling between neuronal activity and the hemodynamic response function (HRF). Indeed, many studies have characterized BOLD-derived measures such as functional connectivity (FC) by assuming that HRF is fixed across subjects and brain regions. This hypothesis of hemodynamic equivalence can be accepted when dealing with young healthy adults, but may not hold in older individuals (Geerligs et al., 2017), in which vascular differences can severely alter the neurovascular coupling. Given its critical role in defining cerebrovascular integrity (Tomoto et al., 2020), the impact of age on the vascular component of BOLD is one of the greatest sources of signal variability. As a result, a failure to consider age-related changes in the neurovascular coupling can induce a bias in the study of the aging brain functioning, other than a misunderstanding of its cognitive relevance (Tsvetanov et al., 2021b). A full quantification of the impact of age on the HRF is then needed since a spatial- and age-invariant canonical HRF would not accurately represent the variability induced by aging (G. Chen et al., 2023).
Among the proposed approaches to estimate the HRF, those involving the computational modeling of the fMRI signal have been widely investigated (Tsvetanov et al., 2021b). While the unpredictable resting-state brain activity cannot be inferred with the traditional general linear model (Raichle, 2001), decomposition approaches may be implemented to detect non-neuronal sources of BOLD activity (Tong et al., 2015; Tsvetanov et al., 2018). However, the selection of the spurious signals may vary depending on the chosen statistical criterion (Tsvetanov et al., 2021b). Scaling methods such as the resting-state fluctuation amplitude (RSFA) have been proposed as an index of vascular reactivity, but age-related RSFA differences cannot be fully explained by cerebral blood flow (CBF) modifications, making it dangerous to employ it as a normalization signal (Tsvetanov et al., 2021a). Additionally, many works have employed different deconvolution methods (Karahanoǧlu et al., 2013; Wu et al., 2013; Caballero-Gaudes et al., 2019; G. Chen et al., 2023). While some techniques evaluate the BOLD signal against a given threshold to identify resting pseudoevents, possibly introducing detection errors which tend to bias the recovery of HRF shape (Wu et al., 2013), others employ a fixed HRF shape throughout the brain (Karahanoǧlu et al., 2013), limiting the description of the high variability of HRF within and across subjects.
Differently from the previous approaches, dynamic causal modelling (DCM) implements a Bayesian inference procedure to effectively disentangle the neuronal and vascular parameters from the fMRI signal associated with each brain region (Friston et al., 2003; Daunizeau et al., 2011). Very few studies have investigated the relevance of DCM-based HRF in neurocognitive aging, focusing on specific regions and/or vascular features (Tsvetanov et al., 2016; Anderson et al., 2020).
Here we propose the use of a sparse version of DCM (Prando et al., 2020) to fully characterize the HRF at the subject level and separately for each region of interest. We then extracted a set of vascular features from each HRF curve (West et al., 2019) to explore their spatial variability in two groups of young and old participants, thus providing a whole-brain quantitative evaluation of the diverse hemodynamic sensitivity to healthy aging. We finally examined the hypothesis that certain hemodynamic features hold significance in predicting age by training an age-classification model and testing it on an independent second dataset of healthy subjects.
As far as we know, this is the first study to provide a quantitative evaluation of the spatial pattern of the hemodynamic sensitivity to advancing age, as well as to give a whole-brain overview of the role of the hemodynamic component as an independent predictor of biological aging, rather than a simple normalization or confounding variable.
Materials and Methods
We analyzed resting-state fMRI (rs-fMRI) data in 464 normal adults between the ages of 18−100 years from two independent cohorts. The first dataset was retrieved from the publicly available MPI-Leipzig Mind–Brain–Body (MPI-LMBB) dataset (Babayan et al., 2019), while the second data collection includes subjects recruited from the Human Connectome Project in Aging (HCP-A) and Development (HCP-D; Harms et al., 2018). In the present study, the MPI-LMBB was employed as the primary dataset for model derivation and internal validation (i.e., “Dataset 1”), while the testing dataset (i.e., “Dataset 2”), crucial for evaluating the generalizability of the derived predictive model, was formed through a combination of data from HCP-A and HCP-D. “Dataset 2” was exclusively employed in the final step of statistically evaluating the classification performance of the model, as detailed below, Statistical evaluation of classification performance.
Participants
Dataset 1
The first rs-fMRI dataset employed in this study consists of resting-state scans of 295 healthy subjects from the publicly available MPI-LMBB dataset (https://legacy.openfmri.org/dataset/ds000221/). The data selection was performed on the original dataset (consisting of 318 individuals) by excluding 23 participants with high motion [<400 volumes with a mean framewise displacement (FD) of <0.4 mm; Power et al., 2012] or affected by preprocessing failures and/or unavailability of rs-fMRI data (Moretto et al., 2022). Then, while the first half of the dataset (147 individuals) was employed for clustering purposes (see detailed description below, Data acquisition and preprocessing), a fraction of the second half was used for the setup of the model after discarding 22 middle-aged participants (31−59 years of age). As a result, a final sample of 95 younger (20−30 years old) participants and 31 older (60–80 years old) participants was included in the study.
Dataset 2
HCP-A and HCP-D are multisite projects designed to generate an in-depth imaging, behavioral, and biosample repository of typical brain changes across a broad portion of the human lifespan (Harms et al., 2018), with HCP-D recruiting participants from 5 to 21 years of age and HCP-A enrolling individuals from 36 to >100 years of age. Both datasets used in this study were drawn from the second release (Lifespan HCP 2.0 Data Release) and focused on typical brains only since HCP studies excluded participants who have been diagnosed and treated for major neuropsychiatric, neurological disorders or depression, as well as those affected by cognitive impairment, learning disabilities, or severe medical conditions (Harms et al., 2018). We did not consider young subjects collected from the former HCP-Young Adult study (Van Essen et al., 2012) since this dataset differs in some important scanning protocol choices from HCP-A/HCP-D that can entail differences in the downstream measures derived from MRI (Harms et al., 2018). From both datasets, we extracted data from 641 subjects aged 18–21 and 36–100 years from HCP-D and HCP-A, respectively, in order to approximately match the age range associated with the training dataset (NDA study DOI, http://dx.doi.org/10.15154/hk8k-pm84). As for the MPI-LMBB dataset, subjects affected by high head motion (<400 volumes with mean FD of <0.4 mm) were discarded (249 participants), as well as middle-aged individuals with age ranging from 41 to 59 years (54 participants). As a result, data from 338 individuals (180 young 18–40 years and 158 old 60–100 years participants) were available for the testing stage.
For both datasets, demographic details and repetition time (TR) are described in Table 1, while Figure 1 provides a schematic representation of the whole pipeline.
Overview of the key processing steps for the implementation of the age-classification model to assess the predictive power of sDCM-derived hemodynamic features in healthy aging. D1, Dataset 1; D2, Dataset 2; HPF, high-pass filter; BPF, band-pass filter; yo, years old.
Participant demographics and TR in the considered datasets
Data acquisition and preprocessing
Dataset 1
Data acquisition was performed with a 3 T Siemens Magnetom Verio scanner, equipped with a 32-channel head coil. The protocol included a T1-weighted 3D magnetization-prepared 2 rapid acquisition gradient echoes (T1w 3D-MP2RAGE) and rs-fMRI images acquired using a 2D multiband gradient-recalled echo echo-planar imaging (EPI) sequence [TR, 1,400 ms; flip angle, 69°; voxel size, 2.3 × 2.3 × 2.3 mm; volumes, 657; multiband factor, 4] and two spin echo acquisitions (Fig. 1, Box 1a). During rs-fMRI scans, the subjects were asked to keep their eyes open and to lie down as still as possible. Further details about the acquisition protocol can be found in Babayan et al. (2019).
The structural preprocessing pipeline included bias field correction (N4BiasFieldCorrection; Tustison et al., 2010), skull-stripping (MASS; Doshi et al., 2013), and nonlinear diffeomorphic registration (Avants et al., 2011) to the symmetric MNI152 2009c atlas (Fonov et al., 2011). Preprocessing of rs-fMRI data consisted of slice timing (Smith et al., 2004), distortion topup (Andersson et al., 2003), motion correction mcflirt (Jenkinson et al., 2002) and nonlinear normalization to the symmetric MNI atlas (Fonov et al., 2011) passing through the single subject pseudo-T1w image, which was obtained by multiplying the T1w 3D-MP2RAGE image with its second inversion time image (boundary-based registration; Greve and Fischl, 2009). As a second step, the GIFT toolbox (http://trendscenter.org/software/gift/) was used to decompose the functional preprocessed data into independent components (ICs) by performing a group spatial-ICA. The ICs were manually classified by visual inspection of both the spatial maps and the source power spectra, in accordance with Damaraju et al. (2014) and Griffanti et al. (2014). As a result, ICs that were related to banding artifacts and vascular or noise components were discarded. Further details on the IC selection can be gathered in Silvestri et al. (2022). Then, 10 principal components (PCs) related to CSF and white matter signal (five from WM, five from CSF) were regressed out from rs-fMRI time series as well as the six standard head motion parameters and their temporal derivatives (Jo et al., 2013; Fig. 1, Box 2).
Dataset 2
Neuroimaging data of HCP-A/HCP-D were collected with a Siemens 3 T Prisma system, equipped with a 32-channel head coil. Anatomical T1w images were acquired via a multiecho magnetization-prepared gradient-echo sequence, while rs-fMRI images were acquired using a 2D multiband gradient-recalled echo EPI sequence (TR, 800 ms; flip angle, 52°; voxel size, 2.0 × 2.0 × 2.0 mm; volumes, 478; multiband factor, 8; Fig. 1, Box 1b). Two sessions of eyes-open rs-fMRI (REST1 and REST2) were performed with opposite phase-encoding direction (four different runs). Additionally, a pair of opposite phase-encoding spin echo were acquired and used to correct functional images for signal distortion. Further details about the acquisition protocol can be found in Harms et al. (2018).
Concerning data preprocessing, the data package employed in the study is the preprocessed volumetric version provided by HCP (Glasser et al., 2013), which contains cleaned rs-fMRI files aligned across subjects using nonlinear volume registration. The dataset includes outputs of “HCP functional preprocessing” for resting-state scans, which is the result of applying Generic fMRI Volume Processing Pipeline, linear detrending and multirun FIX (MR-ICA-FIX) implemented in the HCP pipelines (Glasser et al., 2018). For each subject, only the first run (encoded as rfMRI_REST1_AP_hp0_clean) was employed in the following analysis to ensure phase-encoding consistency with “Dataset 1”.
Common preprocessing steps
Application of the functional atlas
For both datasets, we extracted time series data from a 100-area parcellation scheme of the cortex provided by the Schaefer atlas (Schaefer et al., 2018), which maps to seven resting-state functional networks (RSNs) referring to both hemispheres: visual network (Vis), somatomotor network (SomMot), dorsal attention network (DorsAttn), saliency/ventral attention network (SalVenAttn), limbic network (Limbic), control network (Cont), and default mode network (DMN/Default).
We also defined a set of 10 subcortical regions based on the AAL2 segmentation (Rolls et al., 2015). For each hemisphere, we selected five regions consisting of thalamus proper, caudate, putamen, pallidum, and hippocampus (Fig. 1, Box 3).
We then assigned to each subject a binary temporal mask accounting for brain volumes corrupted by head motion (FD > 0.4 mm), and we applied volume despiking to the time series by means of the icatb_despike_tc function of the GIFT toolbox. Moreover, the temporal traces were bandpass filtered (0.0078–0.2 Hz) so as to retain the common range of frequencies that is usually supported by the canonical HRF in the DCM framework (Henson and Friston, 2007) and spanned by the hemodynamic component at rest (J. E. Chen and Glover, 2015; Fig. 1, Box 4). Furthermore, as part of the normalization process, if the range of the fMRI data exceeded four units, we performed a per-subject rescaling to ensure a consistent range of four, aligning with the recommendations outlined in Zeidman et al. (2019).
Consensus clustering
Given the need to keep the computational load of sparse DCM (sDCM) at a reasonable level (Prando et al., 2020), a consensus clustering procedure was applied to reduce the number of cortical parcels derived from the Schaefer atlas and at the same time to keep the hemodynamic analysis in a whole-brain setting (Fig. 1, Box 3). As a result, we applied to the first half of the MPI-LMBB dataset (147 subjects) a Consensus Clustering Evidence Accumulation procedure, similarly to what described in Ryali et al. (2015), to determine the number of optimal clusters in rs-fMRI time series data. Specifically, this framework employs base and consensus clustering methods to get robust and stable clusters. First, various groupings are obtained for each individual dataset using base clustering method; then, consensus clustering is applied on these groupings to get stable clusters across individuals. In order to account for hemodynamic differences across spatially distant parcels, this procedure was performed selectively for subsets of adjacent cortical regions referring to the same functional network. This additional constraint ensured that only functionally homogeneous and spatially contiguous parcels could be grouped together. For sets composed of only two contiguous regions of interest (ROIs), the objective criteria for determining the optimal number of clusters could not be applied. In these cases, the grouping was performed if the pair Pearson's correlation was on average greater than the mean correlation within the other clusters. This approach ensured that the hemodynamic consistency of each cluster was preserved, contributing to reflect the whole-brain vascular complexity (Taylor et al., 2018; G. Chen et al., 2023) underlying variations of regulatory neurochemical processes (Rangaprakash et al., 2018). The resulting clustering procedure provided 62 cortical clusters, from which demeaned fMRI time courses (i.e., within-cluster mean BOLD signal) were extracted and supplied as inputs to sDCM together with the BOLD signals from subcortical sources (Fig. 1, Box 5). A comprehensive overview of the ROIs derived from the consensus clustering process (presented in the order used in the subsequent analysis) was provided in Extended Data Table 1-1 in the GithHub folder (see below, Data and code accessibility), together with the clustered functional atlas.
Hemodynamic inference
sDCM
For each subject, the resting-state hemodynamic responses and the effective connectivity (EC) matrix were inferred using the sDCM framework (Prando et al., 2020; Fig. 1, Box 6). More specifically, thanks to a discretization and statistical linearization procedure, the inversion of such a linear generative model allows to estimate subject-level parameters accounting for the hemodynamic variability across brain areas, differently from other DCM approaches holding a fixed HRF (Handwerker et al., 2004). In this view, intrasubject hemodynamic variability is captured by the formal modeling of the neurovascular coupling process through differential equations that mathematically describe the link of the hemodynamic and neuronal connectivity parameters with the measured fMRI signal. More in detail, the hemodynamic linearization and discretization procedure consists in translating empirical priors on the physiological parameters describing the nonlinear model of the hemodynamic response (Buxton et al., 1998) into a hemodynamic prior which is exploited in the estimation process of the HRF. This is motivated by the low temporal resolution of fMRI data, which usually ranges from 0.5 to 3 s, and by the idea that the hemodynamic response can be modeled as a finite impulse response model with the neuronal state as input and the BOLD signal as output. In order to ensure that the length of the input response was large enough to model relevant temporal dependencies, we set the hemodynamic length equal to 16.8 s for both datasets, verifying that the number of samples was sufficient to capture the region-wise HRF dynamics. Concerning connectivity parameters, this DCM variant implements a sparsity inducing mechanism that automatically prunes irrelevant connections, thereby avoiding the need to perform a selection between competing network structures (Prando et al., 2020). The algorithm has been further adjusted to account for the signal reliability of the temporal frames by introducing the binary temporal mask as a weighting measure during the estimation procedure. For further insights into model implementation and the inference procedure, readers are encouraged to refer to Prando et al. (2020). It follows that the application of sDCM resulted in the estimation of the EC matrix and 72 ROI-specific HRF profiles (referring to 62 cortical and 10 subcortical nodes) at the individual level, encoding the spatial variations of the neurovascular coupling.
To assess model fit, we employed two summary metrics introduced by Prando et al. (2020) to gauge the reliability of the inferred state-space DCM models and their corresponding parameters for both “Dataset 1” and “Dataset 2.” First, for each subject and ROI, we calculated the normalized root mean squared error (NRMSE) between the actual BOLD signal and that predicted by the individual-level model. Normalizing the RMSE facilitates comparison across datasets or models with varying scales. Specifically, we opted for interquartile range normalization (i.e., the difference between the 25th and 75th percentile of the observations) due to its reduced sensitivity to outliers. Second, we assessed the across-subject distribution of Pearson's correlation between empirical FC (eFC) and that simulated by the inferred model (sFC).
Extraction of the hemodynamic features
First, we applied a PC analysis (PCA) on all the estimated HRF curves to provide a coarse overview of the hemodynamic variability captured by sDCM across all subjects and brain regions.
Then, as carried out in West et al. (2019), nine features defining the HRF shape and timing properties were computed from each individual ROI-level HRF (Fig. 1), including peak amplitude (“A peak”), time-to-peak (“t peak”), rise slope, trough amplitude (“A trough”), time-to-trough (“t trough”), fall slope, time from peak to trough (“Peak trough”), full-width at half-maximum (“FWHM”) of the peak, and area under the curve (“AUC”; Fig. 1, Box 7). It follows that 648 variables (nine features for 72 ROIs) were considered in the subsequent fitting procedure. For “Dataset 2,” the extraction of hemodynamic parameters was conducted after temporal interpolation of the hemodynamic responses to align with the TR of the training dataset. Further details about the technical implementation of the feature extraction step can be found in the script “extract_HRF_features.m” (see below, Data and code accessibility).
Age-classification model
In order to assess the prediction power of the hemodynamic features, we trained a logistic regression model classifying young-age versus elderly subjects (Fig. 1, Box 8) as follows:
Initially, a univariate method was employed to systematically examine the spatial sensitivity of brain hemodynamic characterization features with respect to healthy aging. Then, performing a classification task between young and elderly individuals stems from the objective of evaluating whether such changes involving hemodynamic features and their interacting effects can effectively predict the age class rather than merely serving as confounding variables.
In technical parlance, our methodological pipeline for identifying the optimal predictive model incorporated an iterative hybrid feature selection strategy (Jović et al., 2015; Remeseiro and Bolon-Canedo, 2019; Mochammad et al., 2021), wherein we initially applied a less computationally intensive filter to eliminate certain features (e.g., univariate selection). Following this, a more computationally intensive wrapper method, such as stepwise backward multivariate selection, is employed for fine-tuning. This process occurs subsequent to the removal of multicollinear features, which serves as an additional filter step aimed at discarding highly redundant features. By leveraging such a hybrid approach, computational time is minimized through the elimination of redundant features via the univariate method to diminish the dimensionality of the feature space. Then, the classification model accuracy is enhanced by exhaustively exploring all potential candidate features among the top-ranked ones using the wrapper method, which is deployed to pinpoint the most promising candidate subset maximally contributing to the accuracy of the classification model (Jović et al., 2015). By combining the strengths of both methodologies, hybrid methods typically yield high accuracy akin to wrappers and offer high efficiency akin to filters (Jović et al., 2015; Mochammad et al., 2021).
Moreover, the choice of training the model on the dataset with lower temporal resolution comes from the purpose of investigating hemodynamic factors that could be detectable from most commonly available MRI scanners. Then, testing the model on HCP data collected with a higher sampling rate allowed to assess its potential applicability to a wider range of fMRI acquisition and preprocessing pipelines.
Partitioning of the dataset
To overcome the risk of overfitting the predictive model on a predefined partition of the training dataset and to reduce the sensitivity of the model to outliers, the fitting procedure was iterated over 1,000 random partitions. Specifically, in each repetition we adopted a 20% hold-out validation approach to split the “Dataset 1” into two parts, training on one (i.e., training subset) while validating on the other with the trained model (i.e., validation subset). Given the unbalanced size of the two age classes, we performed a stratified split to maintain the class proportion in the two sets. The Pearson's correlation coefficient of each training covariate with age was also assessed to later consider it as a quality metric of the given data partition in the fitting of the final logistic regression model (see below, Minimal model identification). Then, the outputs of the single models were combined and statistically evaluated to define a robust final “minimal” logistic regression model.
Univariate feature selection
Given the high number of identified variables compared with the sample size of “Dataset 1,” the first step regarded the application of a univariate significance test for selecting features that significantly vary across the two age groups derived from the training subset. This step aimed at discarding predictors with negligible distribution differences between young and old populations, thus streamlining the set of considered features that exhibit significant changes and reducing the computational load of the following selection stages. For each of the 648 hemodynamic features, we conducted a two-sample t test by considering the two samples as extracted from two normal populations with unequal variances (e.g., for the parcel LHVis1, we statistically compared the distributions of each hemodynamic feature referring to the two populations). A Holm–Bonferroni’s correction was applied to control for multiple comparisons (p < 0.05). As a result, only features that show significant age differences were considered in the following fitting procedure.
Multicollinearity removal
While univariate filter methods can operate independently of machine learning models by evaluating and typically ranking individual features (Jović et al., 2015), they may overlook the correlation among selected features, potentially compromising model reliability due to multicollinearity (Mochammad et al., 2021). Strong correlations, indicative of collinear data, can result in counterintuitive and uninterpretable coefficients, destabilizing the model fit to the data (Belsley et al., 1980). For this reason, we proceeded with a correlation analysis to mitigate the collinearity effect among independent variables. Specifically, we assessed the Pearson's correlation between each pair of hemodynamic features retained in the prior selection stage. To alleviate very high redundancies (absolute Pearson's ρ > 0.85), we retained only the variable most correlated with the output. Specifically, for each pair of redundant features, we removed the feature that exhibited the least correlation with the “continuous version” of the output, i.e., the age values of subjects. This process was reiterated until no highly correlated pairs of features remained. The choice of a high correlation threshold at 0.85 was intentional, serving as a relatively high value to avoid overly aggressive variable removal, considering the subsequent multivariate selection step.
To gauge the effectiveness of this selection process, we employed a systemic measure for collinearity known as the condition number (CN), as proposed by Belsley et al. (1980). The CN, computed as the ratio of the largest to the smallest singular value of the design matrix X (i.e., the training matrix of independent predictors), serves as a valuable indicator of data collinearity (Tomaschek et al., 2018). A large CN signals significant collinearity effects, where even minimal differences in response or predictor variables lead to substantial consequences for estimated regression coefficients.
Multivariate feature selection
Subsequently, a multivariate selection of features was performed on the training subset. To further increase the statistical robustness of the final selection, we executed a bootstrap resampling 100 times to create multiple internal training sets. At each iteration, we generated the internal dataset by sampling with replacement from the current training subset a collection of data with the same cardinality and by maintaining the class proportion in the final selection of subjects. Then a stepwise backward feature selection was applied to each sampled normalized internal set by using the Bayesian information criterion (BIC) as the stopping rule. More specifically, a generalized linear model was fitted to each internal dataset using stepwise regression to add or remove predictors (i.e., the hemodynamic features) starting from the linear model containing an intercept and linear term for each predictor. The stepwise backward selection method initiates with backward selection and switches to forward selection whenever a variable is omitted, assessing whether any previously omitted variables should be included in the model. At each step, the function scans for terms to incorporate or exclude from the model, guided by fluctuations in the BIC as the primary selection criterion. The adoption of a stepwise backward selection ensures that all variables are initially considered in the model, capturing potential initial interaction effects differently from the forward approach.
Fit of the logit model
The robustness of each iteration-level final selection was achieved by finalizing the selection of covariates retained in a substantial portion of the internal training subsets. More specifically, the hemodynamic covariates that were selected in a significant portion (>50%) of the internal training sets in the previous stage were considered to fit a logistic linear model predicting the probability of a subject to be classified as young (pyoung). The model predictive performance was evaluated on the validation subset by computing the corresponding concordance index (c-index), which is numerically equivalent to the area under the receiver operating characteristic curve (AUCROC) as follows:
Minimal model identification
Having estimated a distinct classifier for each random partition, we combined single results to determine the hemodynamic features that were most frequently selected in the 1,000 iterations (Fig. 1, Box 9). First, we filtered out atypical model realizations by assessing the data representativeness of each training subset, which was quantified by computing the Euclidean distance between the correlation vectors of the training and the whole dataset with the output variable. We further excluded any model that performed worse than a random classifier on the validation subset (AUCROC or c-index < 0.5). We finally reduced the number of plausible variables by discarding selected features whose corresponding coefficient estimate was not statistically significant, according to a p > 0.05. This resulted in considering only features selected by the best-fitted models, which were ordered in descending fashion based on their selection frequency. On the basis of the identification of the knee point in the frequency curve (i.e., significant slope changes of the frequency curve), the most frequently selected variables were then employed to define the final minimal logistic regression classifier. The final model predicting the age group was then retrained on the 1,000 random sets to get the median regression coefficient for each hemodynamic feature, after selecting only good-fitted model realizations (AUCROC > 0.5; p < 0.05).
Statistical evaluation of classification performance
The final step was then to assess the performance of the minimal model on the testing dataset (“Dataset 2”). This evaluation was performed by computing the final ROC curve and the associated AUC (AUCROC) to provide an aggregate measure of performance across multiple classification thresholds. The PR curve and the associated AUC (AUCPRC) were also considered since it can be more informative than the ROC plot when evaluating classifiers on imbalanced datasets. Since the raw prediction of the minimal logistic model returns a continuous output variable within the range [0, 1], a cutoff value is essential to perform the final classification. As a result, ROC and PR curves were generated by assessing model performance across various cutoff values (i.e., spanning the array of young age probabilities predicted by the minimal model). Then, we selected the cutoff maximizing balanced accuracy (BA) to illustrate the classifier's output in terms of the confusion matrix for the best-performing scenario. Specifically, we assessed the BA as it is known to be more effective than accuracy in situations involving imbalanced classes (Haixiang et al., 2017). Additionally, we employed sensitivity, specificity, and precision to delineate the ROC and PR curves as follows:
Similarly, the same metrics were assessed on each validation subset and on the whole “Dataset 1” to provide sample confidence intervals for the performance in the random train-test split and an upper bound for the overall accuracy, respectively (Fig. 1, Box 10).
Data and code accessibility
This study is a secondary analysis of publicly available datasets. The download links for the datasets are provided above. No personal identifiable information was used in this study. The analyses were conducted utilizing custom-built code developed in-house and executed in MATLAB 2021b (MathWorks). The complete set of code required to execute the analysis pipeline and produce figures is openly accessible through the following link: https://github.com/FairUnipd/HRF_aging.git.
Results
Region-wise HRF
In both datasets, the model fit was generally satisfactory (“Dataset 1”: whole-brain median NRMSE, 33.91%; median eFC–sFC corr, 0.95; “Dataset 2”: whole-brain median NRMSE, 22.49%; median eFC–sFC corr, 0.91), indicating that sDCM effectively captures the physiological variability and system dynamics derived from the input BOLD data. Once the reliability of the HRF estimates is assessed, using PCA decomposition, we reported evidence of the variability of the HRFs inferred through sDCM across all subjects and brain regions. Figure 2 illustrates the distribution of HRFs as described by the first three PCs in each of the two age samples (“Dataset 1”). Additionally, the figure provides insights into the percentage of total variance explained and outlines the trends observed in the corresponding PC loadings. Although the spatial localization of clusters in terms of functional networks is quite reproducible in both groups, it is worth noticing that ROI scores differ significantly across the brain space and subjects on the basis on their hemodynamic features, providing clear evidence of the high variability characterizing HRF. As depicted in the panels illustrating HRF curves, hemodynamic profiles not only exhibit differences across aging groups but also show a strict dependence on the localization of the brain area under investigation. For example, an inversion of the amplitude difference is observed when moving from a cortical to a subcortical parcel. This observation motivates the subsequent ROI-wise analysis of hemodynamic curves.
The three-dimensional space identified by the first three PCs allows to describe the HRF variability across subjects in the two age samples and across ROIs. Each point condenses the hemodynamic description of an individual parcel. The color denotes network affiliation or the corresponding subcortical node. Points marked with a black outline refer to an example ROI, whose mean (solid line) and standard deviation (shading) of the HRF curves are plotted in the central-left and bottom-left panels, corresponding to young (blue) and old (red) healthy adults in “Dataset 1.” Furthermore, the central-right and bottom-right panels display the trend of the corresponding PC loadings, accompanied by the percentage of total variance explained by each PC (PC1, 58.8%; PC2, 32.2%; PC3, 6.8% for cortical ROIs; PC1, 70.5%; PC2, 20%; PC3, 7.1% for subcortical ROIs). Three-dimensional hemispheric and subcortical renderings in this figure and in the following ones were obtained by using the ENIGMA toolbox (Larivière et al., 2021).
Spatial and group variability of hemodynamic features
Figures 3 and 4 display the resulting median maps of hemodynamic features characterizing the two age groups. Both cortical and subcortical maps exhibit spatial heterogeneity. Specifically, the median peak amplitude varies substantially across ROIs, with higher values in visual and frontoparietal areas, as well as in caudate and pallidum subcortical nodes. However, this feature shows a general decreasing trend with aging which involves all brain regions. Similarly, trough amplitude (i.e., the poststimulus undershoot) is quite variable across ROIs with a spatial pattern that is inversely proportional to that associated with peak amplitude in both cortical and subcortical regions, meaning that a higher hemodynamic peak is often related to a lower (i.e., deeper) poststimulus undershoot. Concerning AUC, we found higher values in cortical visual and some SomMot areas as well, whose associated pattern highly resembles the one characterizing the peak amplitude. In addition, AUC shows high values in the caudate nucleus and putamen, while it remains low in the hippocampus and thalamus. These spatial profiles are mostly preserved in the older subjects, even if characterized by a slight decrease in strength. Together with peak amplitude, this feature entails a pronounced hemispheric asymmetry (slightly declining with age) described by a higher spatial heterogeneity regarding the right side of the brain. Regarding rate features (fall and rise slope), we found similar cortical maps (whose pattern recalls the one revealed by AUC) but opposite in sign, with an increasing (decreasing) trend for fall (rise) slope in the group of elders. Interestingly, this trend is reversed for subcortical nodes. Concerning response width (FWHM), the young sample shows a broader FWHM in frontoparietal areas and in the putamen and caudate nuclei, while the lowest values are encountered in the hippocampus and limbic areas in both groups.
Median spatial distribution cortical maps for the HRF features computed as described in Figure 1. Blue boxes refer to the young sample, while red boxes are related to the old sample (“Dataset 1”). The same scale is adopted for each feature to ease the comparison between the age groups.
Median spatial distribution subcortical maps for the HRF features computed as described in Figure 1. Blue boxes refer to the young sample, while red boxes are related to the old sample (“Dataset 1”). The same scale is adopted for each feature to ease the comparison between the age groups.
On the other hand, timing features (i.e., time-to-peak, time-to-trough, peak-to-trough) generally show a weaker variability than the previous hemodynamic features. Time-to-peak varies very little across the cortex and is stable across subcortical nodes, revealing the reduced sensitivity of this feature to age-related changes. In the older subjects, time-to-trough shows a pronounced increase in specific SomMot and DorsAttn areas in both hemispheres. Concerning subcortical areas, we observed the highest values in the right caudate nucleus and putamen in the younger group, while the lowest values characterize the left hippocampus and pallidum in the elders. Regarding peak-to-trough, its cortical median map exhibits a very similar cortical pattern. However, a slightly different spatial scheme is encountered in subcortical nuclei, where the right caudate presents the highest value in the younger sample, while the hippocampus and pallidum in older participants are interested by the lowest peak-to-trough.
Univariate feature selection
The univariate analysis involving HRF differences between younger and older groups across the 1,000 iterations yielded distinct results depending on the selected ROI and the feature of interest (Fig. 5a) Specifically, significant differences were mostly observed in areas of lower-order functional networks (i.e., visual and somatosensory motor), as well as in more cognitive RSNs such as DorsAttn, Cont, and DMN. Selection frequency is rather hemispheric balanced (left-hemispheric ROIs, 48.15%; right-hemispheric ROIs, 51.85%), although the most selected HRF feature may vary based on the considered ROI. However, it is worth noticing that amplitude (i.e., peak amplitude, trough amplitude, and AUC) and rate features (fall and rise slope) generally prevail over the others in terms of sensitivity to age-related hemodynamic changes. Concerning subcortical nodes, thalamic areas are the most frequently selected for significant differences in peak amplitude and fall slope. For each variable with a selection frequency > 0, Figure 5b reports the corresponding values of difference between the group medians. While peak amplitude, AUC, and rise slope are generally greater in the younger group, undershoot amplitude is usually damped (i.e., higher values) in elders, as well as the fall slope. Note that for these features, the general trend is reversed in subcortical areas. On the contrary, FWHM differences are consistent in sign across the brain space. As expected, temporal features show very weak differences (mainly values prevailing in the old group) across nodes, which are mainly localized in Vis, DorsAttn, and some subcortical parcels.
Results of the univariate feature selection step. a, Frequency plot illustrating the percentage of frequency selection for each feature across the 1,000 dataset partitions, reporting a balanced trend across hemispheres. Each marker refers to a different hemodynamic feature (as indicated in the legend), while the background color code discriminates between functional networks and/or subcortical ROIs. Within each cortical network, ROI numeric labels are organized in ascending order, while subcortical nodes are listed in a left-to-right sequence. b, For each feature–ROI pair characterized by a selection frequency >0 in panel a, the stem plot reports the difference between the corresponding median values in the two age samples as extracted from the maps in Figures 3 and 4. Positive values indicate a significantly higher estimate in the young group (blue area), while negative values refer to a significantly higher estimate in the old group (red area). As in panel a, each marker refers to a different hemodynamic feature and within each cortical network ROI numeric labels are organized in ascending order, while the subcortical nodes are listed in a left-to-right sequence.
Multicollinearity removal
The selection frequency of the preserved features following the multicollinearity reduction step is depicted in Figure 6 (panel a), alongside the distributions of across-subject CNs during the successive removal of highly redundant features (panel b). Remarkably, the frequency landscape remains largely stable compared with the previous stage, showing only marginal reductions in frequency values (e.g., observe the trough amplitude for LHVis2). This stage primarily focuses on eliminating significant redundancies among features within the same ROI. For instance, there is a notable decrease in selection frequency for the fall slope of LHCont2 compared with other features within that region, highlighting the efficacy of this step in refining feature sets. It follows that this retention pattern ensures the maintenance of hemispheric balance in the feature selection process (left-hemispheric ROIs, 48.89%; right-hemispheric ROIs, 51.11%).
Results of the multicollinearity reduction step. a, Frequency plot illustrating the percentage of frequency selection for each feature across the 1,000 dataset partitions, reporting a balanced trend across hemispheres. Each marker refers to a different hemodynamic feature (as indicated in the legend), while the background color code discriminates between functional networks and/or subcortical ROIs. Within each cortical network, ROI numeric labels are organized in ascending order, while the subcortical nodes are listed in a left-to-right sequence. b, Declining trend in across-subject CN distributions with the gradual elimination of highly correlated features. Each data point in the boxplots represents a CN value from 1 of the 1,000 iterations. To enhance visualization, each iterative CN curve was interpolated to achieve a maximum of 10 cycles.
Conversely, the CN trend demonstrates a consistent decline as highly correlated pairs of variables are systematically removed from the internal design matrix derived from the training subset. This progressive reduction signifies the diminishing impact of significant collinearity effects, leading to an overall decrease in the CN.
Multivariate feature selection
As illustrated in Figure 7a, a totally different scenario is provided after the multivariate selection process. First, the results suggest a prominent hemispheric asymmetry since retained hemodynamic features are mostly localized in the right-hand side of the brain for both cortical and subcortical domains (left-hemispheric ROIs, 29.11%; right-hemispheric ROIs, 70.89%). The emergence of a “right bias” becomes evident only during this multivariate selection phase, indicating that the information provided by rightward ROIs plays a major and pivotal role for effective classification through the logistic model.
Results of the multivariate feature selection step. a, Frequency plot illustrating the percentage of frequency selection for each feature across the 1,000 dataset partitions, reporting a lateralized trend in favor of right-hemispheric areas. Each marker refers to a different hemodynamic feature (as indicated in the legend), while the background color code discriminates between functional networks and/or subcortical ROIs. Within each cortical network, ROI numeric labels are organized in ascending order, while the subcortical nodes are listed in a left-to-right sequence. b, Stem plot reporting the multivariate (first stem) and univariate (second stem) selection frequency values for each of the top 10 (to maintain consistency with the results illustrated in panel c) hemodynamic features (in descending order). The blue area highlights the features that were included in the final minimal predictive model (Fig. 8). c, Boxplots depicting the distribution of data representativeness values versus model complexity trends for the 1,000 trained models (shown as blue points).
On the one hand, this result contrasts with task-based aging models such as hemispheric asymmetry reduction in older adults (Cabeza, 2002) or compensation-related utilization of neural circuits hypothesis (Reuter-Lorenz and Cappell, 2008), supporting a decreased hemispheric lateralization with advancing age, which can be viewed as a specific form of compensatory processes manifesting during task execution. On the other hand, the identified age-induced asymmetry is consistent with theories associated with the right-hemispheric aging model (Dolcos et al., 2002), indicating that age-related cognitive declines disproportionately affect the right hemisphere.
Specifically, the most predicting variables are found in Vis, SomMot, DorsAttn, Cont, and DMN networks. Concerning subcortical nuclei, the right thalamus and pallidum are the only nodes that play an important role in the model prediction. As regards feature class, amplitude features prevail after multivariate selection, corroborating their crucial contribution in predicting age. Figure 7b outlines the top 10 HRF features in descending order according to their multivariate selection frequency, together with the corresponding univariate selection frequency computed in the previous step. It is evident that the two selection stages follow two completely different trends, motivating the importance of accounting for multivariate interactions of hemodynamic features for model fitting. Despite the more abrupt slope change following the second top-selected feature, we also observed a noticeable slope change in the multivariate selection curve after the first five features (i.e., second knee point), which guided the identification of these predictors to be included in the final minimal logit model. The decision to choose the second knee point was informed by the observation that the vast majority of models with optimal data representativeness (i.e., low Euclidean distance between the correlation vectors of the training subset and the entire dataset with the output variable) usually retained between three and five features (i.e., model complexity ranging from four to six when including the intercept term), as evidenced in Figure 7c.
Minimal model identification
Figure 8a depicts the distributions of regression coefficients obtained by fitting the minimal model on the 1,000 random sets. Despite the presence of outliers, the width of the boxplots is quite limited, suggesting a good consistency of estimated β values in terms of variance across model realizations despite the physiological intergroup variability observed across the random training subsets (Fig. 7c). Regarding cortical parcels, the four selected ROIs are mostly localized in the right frontoparietal cortex [involving dorsal attention “PrCv,” “FEF” (“DorsAttn5”), “Post5” (“DorsAttn4”), and control “PFCl3” areas] and in a great part of the right visual central/temporal–occipital region. On the other side, right thalamic nucleus is the only subcortical node implicated in the model prediction. To what concerns HRF feature class, we could deduce that the age groups can be discriminated on the basis of the AUC in the first dorsal attention ROI (positive relation with logit of pyoung), peak amplitude of the right thalamus and the visual parcel (negative and positive relation with logit of pyoung, respectively), peak to trough for the second dorsal attention ROI (negative relation with logit of pyoung), and trough amplitude in the case of the control network region (negative relation with logit of pyoung).
Results of the minimal model identification. a, Distributions of the regression coefficients obtained by fitting the minimal model across the 1,000 set realizations. Each marker refers to a different hemodynamic feature (as indicated in the legend), while the color of the boxplots discriminates between functional networks and subcortical ROIs. The corresponding coefficients of variation (CV %) are as follows (from left to right): 30.47, 28.67, 27.51, 33.60 and 48.13%. Spatial localization of the corresponding ROIs (right hemisphere) is provided by the brain surface renderings on the top, while corresponding HRF curves are reported in the below panels showcasing the mean (solid line) and standard deviation (shade) for the young (blue) and old (red) groups. b, PDFs of the performance metrics (AUCROC, AUCPRC, balanced accuracy, and overall accuracy) obtained by evaluating the minimal model realizations across all 1,000 random sets.
On the other hand, Figure 8b displays the estimated probability density functions (PDFs) for performance metrics of the 1,000 model realizations. Overall, the results suggest that the minimal model performs well across all the set partitions in terms of (balanced) accuracy, as well as for precision, specificity, and sensitivity (i.e., AUCROC and AUCPRC) since all distributions mostly range from 0.7 to 1.
Classification performance of the minimal model
Having derived the median values of regression coefficients depicted in Figure 8a, we tested the “median” minimal model both on the whole “Dataset 1” and “Dataset 2.” Figure 9, a and b, displays the ROCs and PRCs in both datasets, respectively. The ROC and PR curves indicate good classification performance for both datasets (i.e., better than the random classifier), although PRC trend generally exhibits a higher variability in “Dataset 2.” Colored dots overlaid on the curves denote the cutoff scores maximizing the BA (specifically, 0.92 in “Dataset 1” and 0.79 in “Dataset 2”), whose corresponding value was employed to identify the confusion matrices in Figure 9, c and d. The asymmetry of the confusion chart in panel c suggests that the classification error is due to some young subjects being classified as old, while all elders are correctly classified by the minimal model. The predictive scenario is more balanced in “Dataset 2,” despite the fact that the percentage of right predictions is slightly higher in the old cohort.
Classification performance of the minimal model. a, ROCs of the minimal model for the prediction in “Dataset 1” (blue) and “Dataset 2” (green), with colored dots denoting the cutoff values for “Dataset 1” (blue) and “Dataset 2” (green); the dashed line corresponds to the ROC curve of the random classifier. b, PRCs of the minimal model for the prediction in “Dataset 1” (blue) and “Dataset 2” (green), with colored dots denoting the cutoff values for “Dataset 1” (blue) and “Dataset 2” (green); the dashed line corresponds to the PR curve [i.e., the positive rate P / (P + N)] of the random classifier. c, Confusion matrix of the minimal model for the prediction of the age class in “Dataset 1”; for each aging group, the class-wise precisions and class-wise recalls (computed as the percentage of the number of correctly and incorrectly classified observations for each predicted and true class, respectively) are indicated. d, Confusion matrix of the minimal model for the prediction of the age class in “Dataset 2”; for each aging group, the class-wise precisions and class-wise recalls (computed as percentage of the number of correctly and incorrectly classified observations for each predicted and true class, respectively) are indicated. D1, Dataset 1; D2, Dataset 2.
Together, these findings illustrate that despite the substantial computational demands of our chosen approach to balance the sample size and model stability, it ensures commendable performance across both datasets. Simultaneously, it ensures that the resulting model remains interpretable (i.e., it does not involve feature extraction methods) and facilitates meaningful insights into the underlying data.
Discussion
Heterogeneity of aging effects on HRF
The high intersubject and regional variability of the hemodynamic response evidenced by PCA analysis prompted us to model HRF locally in the brain and at the individual level (Handwerker et al., 2004; Rangaprakash et al., 2018; Kassinopoulos and Mitsis, 2019). Given the reported modifications of the hemodynamic system with healthy aging (Y. Li et al., 2018; West et al., 2019), we first investigated the hemodynamic variations by delivering a comprehensive depiction of the median spatial pattern for each feature and evaluating significant hemodynamic disparities between the two populations via the univariate selection step.
In general, amplitude features shared correlated trends across the lifespan, with significant age-related decreases for peak height (Hutchison et al., 2013; Fabiani et al., 2014; Ward et al., 2015; West et al., 2019) and AUC, which are reflected in a damping of poststimulus undershoot in the same networks. Analogously, West et al. (2019) reported the same spatial trends for peak and trough amplitude, while it finds an inverse (not significant) tendency for AUC. As expected, spatial variability of rate features varied with aging by a decreasing (increasing) trend for rise (fall) slope, accordingly to West et al. (2019). The emergence of such aging HRF properties possibly relates to a significantly slower neurovascular coupling in elderly people (Cherkaoui et al., 2021). In this view, also timing features generally presented an increasing age-related trend in cortical ROIs, which accords with West et al. (2019).
Interestingly, the patterns of significant age-related changes were generally reversed in subcortical areas, where the HRF is usually faster and narrower than in cortical ROIs (Lewis et al., 2018; Kim et al., 2022). Accordingly, Lee et al. (2009) represented prior findings of augmented subcortical perfusion in healthy older adults. Specifically, aging seems to have a greater impact on thalamic HRF, which presented significant differences in terms of a greater peak amplitude and faster rise and fall phases. This contrasts with Galiano et al. (2020), reporting that perfusion reduces with age in the medial thalamus. Also, the right pallidum, whose AUC appears to significantly increase with age, was the only subcortical nucleus to exhibit a significant decreasing impact on the association between age and CBF in Y. Chen et al. (2011). This lack of concordance between results could derive from the low signal-to-noise ratio in subcortical nodes (Goodro et al., 2012; Singh et al., 2018) and the impact of the age-related volume reduction (Wang et al., 2019).
The prominent role of right-hemispheric regions in the age-classification task
Since univariate age-related differences are largely confined to the local level, we then assessed the combined prediction power of the proposed HRF features that were retained in the univariate selection step, in order to classify young versus elderly subjects using a standard machine learning tool. The reason for performing such classification task proceeded from the aim to assess whether the hemodynamic component could reflect brain aging independently from neural information, motivating its potential role as biomarker of brain age rather than a mere confounder. Moreover, this supervised multivariate learning step allowed ranking the most discriminative features by also considering interacting effects between hemodynamic variables, differently from the previous steps. Previous studies have already reported evidence of the predictive power of HRF in healthy aging by examining a limited number of cortical nodes and networks (Tsvetanov et al., 2016; Anderson et al., 2020) and/or a reduced set of hemodynamic features (Cherkaoui et al., 2021). On the contrary, the proposed predictive framework provided a whole-brain spatial characterization of the hemodynamic sensitivity to healthy aging, by fully describing age-related vascular alterations of the HRF shape in terms of amplitude, rate, and timing features.
Our results showed that right-hemispheric ROIs generally have an improved prediction power to discriminate age compared with their left counterpart. The prominent role of interhemispheric asymmetry in neurovascular coupling emerged only when considering the interacting contribution of predictors (compare Figs. 5–7). In line with this asymmetry, Gauthier et al. (2013) documented a reduction in hypercapnic CBF within the right frontal areas, which is linked to the age-related stiffening of vessels, even in the presence of unchanged neuronal activation. These results may indicate the vascular origin of the right-hemispheric aging model (Dolcos et al., 2002) and other studies reporting hemispheric asymmetries in resting functional activity with age (Beason-Held et al., 2008; Z. Li et al., 2009; Zuo et al., 2010). Based on our observations, these functional outcomes could be biased by the pronounced hemodynamic alterations affecting the rightward ROIs. As a result, this may reflect asymmetric neurodegenerative processes of vascular nature rather than neural plasticity, which imply the possibility of distinct vascular aging processes parallel to neural ones evidenced by task-based studies (Cabeza, 2002; Reuter-Lorenz and Cappell, 2008).
Through the minimal model, we identified the brain regions that are key predictors of the age class. Concerning cortical networks, the observed vulnerability of Cont and DorsAttn to the hemodynamic decline might lead to reconsider previous works reporting decreased within-DorsAttn connectivity and cohesion with increasing age (Betzel et al., 2014; Zonneveld et al., 2019), further suggesting the potential bias induced by the hemodynamic component. Likewise, hemodynamic differences emerged in the Cont network, which was shown to have a general within-network decrease of local activity in older adults (Geerligs et al., 2015; Bethlehem et al., 2020). Interestingly, Avelar-Pereira et al. (2017) examined the interactions between the right Cont and DorsAttn, reporting that older adults show altered connectivity values compared with the young, while connectivity between the left Cont and DorsAttn was not significantly different between the age groups. Specifically, the significant age-sensitive HRF detected in the right “PFCl” recalls the work of Grady et al. (2016), which reported that this region shows stronger connections to the DorsAttn in older adults. Similarly, the peak amplitude detected in right visual central areas and the right thalamus seems to be a crucial predictor of brain age. Accordingly, Beason-Held et al. (2008) pointed out a decreased CBF in occipital and occipitotemporal cortices that may impact visual perception in older subjects. Moreover, the authors illustrated a longitudinal increase of CBF in the right-hemispheric subcortical thalamus, proposing that this trend may reflect preservation of function over time in view of age-related declines in global CBF.
Overall, the minimal model exhibited outstanding predictive performance in both “Dataset 1” and “Dataset 2,” reaching an overall BA of 92 and 79%, respectively, without losing in precision. As a result, the pattern of unveiled hemodynamic predictors may represent a reliable biomarker for brain age that could be replicable for a wide range of fMRI acquisition protocols.
Limitations
fMRI temporal resolution significantly matters for accurate timing and rate feature estimation (Rangaprakash et al., 2018), whose reliability in our study was limited by the high acquisition TR. Especially for subcortical nuclei, a refined HRF assessment would require high spatiotemporal resolution to resolve the small nuclear subdivisions (Singh et al., 2018). Moreover, as the proposed regional analysis was conditioned by the choice of the parcellation atlas, it would be advantageous to examine the reproducibility of the predictive hemodynamic pattern across various atlas variations (Cherkaoui et al., 2021). Third, we did not include EC parameters as candidate predictors of age (Tsvetanov et al., 2016) since our aim was to test the independent predictive ability of vascular factors. However, including connectivity factors could provide a more complete interpretation of the unique and shared contributions of neuronal and vascular components to brain-age prediction.
Regarding modeling limitations, DCM inversion scheme typically treats the temporal and spatial dynamics of the hemodynamic response independently, thus overlooking intrinsic spatiotemporal relations (Aquino et al., 2014). Moreover, the multivariate selection method carries a considerable computational burden, which may be enhanced by integrating embedded selection methods (Jović et al., 2015). Additionally, regression-based methods assume linear effects between response and predictors while lacking a principled approach to incorporate higher-order interactions. Future studies could benefit from exploring feature extraction and deep-learning techniques. However, in such cases the incorporation of a larger sample size becomes imperative. Moreover, the induced feature combinations may lead to a potential loss of physical meaning, thereby undermining the readability and interpretability that are paramount in biomedical domains (Remeseiro and Bolon-Canedo, 2019).
Conclusions
By inferring ROI-wise hemodynamic profile at the individual level through sDCM, we provided a whole-brain characterization of the spatial pattern of brain hemodynamic sensitivity to healthy aging. Through the univariate testing, we evaluated the statistical impact of increasing age on each regional hemodynamic feature, proposing a quantitative assessment of the spatially heterogeneous hemodynamic bias that could affect brain functional assessment. We then unveiled a specific set of hemodynamic predictors to discriminate young from elderly people, mainly describing vascular properties of right-hemispheric areas devoted to cognitive functions. The discerned asymmetry in vascular degeneration processes affecting the aging human brain suggested a spatial pattern that could emerge as a candidate biomarker for predicting the “age” of the brain (Cole et al., 2017; Green and Hillersdal, 2021). Such findings could pave the way for personalized predictions related to neurovascular diseases and mortality risks among older adults.
Footnotes
G.B. and D.B. were supported by the DEI Proactive grant "Personalized whole-brain models for neuroscience: inference and validation” from the Department of Information Engineering of the University of Padova (Italy). The Human Connectome Project (HCP)-Aging dataset reported in this study was supported by grants of the National Institute on Aging of the National Institutes of Health under Award Number U01AG052564 and by funds provided by the McDonnell Center for Systems Neuroscience at Washington University in St. Louis. The HCP-Aging 2.0 Release data used in this report came from http://dx.doi.org/10.15154/1520707. The HCP-Development dataset reported in this study was supported by grants of the National Institute of Mental Health of the National Institutes of Health under Award Number U01MH109589 and by funds provided by the McDonnell Center for Systems Neuroscience at Washington University in St. Louis. The HCP-Development 2.0 Release data used in this report came from http://dx.doi.org/10.15154/1520708.
The authors declare no competing financial interests.
- Correspondence should be addressed to Alessandra Bertoldo at alessandra.bertoldo{at}unipd.it.
