Study participantsi) ADNI dataset

Imaging and phenotypic data were obtained from the ADNI database (https://adni.loni.usc.edu/) [32]. ADNI was launched in 2003 as a public-private partnership led by principal investigator Michael W. Weiner. The primary goal of the ADNI is to test whether serial MRI, positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessments can be combined to measure the progression of MCI and early AD. In the ADNI database, we considered participants who were initially categorized as CN or as having MCI but progressed to AD and censored individuals who did not convert to AD. We applied the following criteria to select participants at risk of progression to AD: (1) availability of multimodal imaging dataset (i.e., T1-weighted [T1w], rs-fMRI) with sufficient image quality (i.e., 3 T), (2) head motion with mean framewise displacement < 0.3 mm, and (3) participants diagnosed with AD within five years of the scan date (CN or MCI → AD) and individuals diagnosed as CN or with MCI during the five-year follow-up period (CN → MCI or MCI → MCI, thus right-censored) and (4) excluded reversion cases (e.g., AD → MCI). The five-year limit of time to conversion was imposed so that the baseline MRI data might possess prognostic power for the conversion event that occurred later, as previously reported [34, 35]. This resulted in the selection of 68 unique participants from the ADNI-GO/2/3 datasets who were at risk of conversion to AD. Due to the limited number of unique participants, we included one session from the study participants that revisited the study (e.g., 6 or 12 months), constituting a total dataset of 115 samples. The sample consisted of participants who were CN or with MCI (mean age, 72.9 years; 51% male), and the median conversion time from baseline to AD was 1,277 days, with a 48% conversion rate. In summary, each sample had a baseline scan and an associated conversion date; however, some samples were obtained from the same participant. We considered demographic and clinical variables, including age, sex, Mini-Mental State Examination (MMSE) score, apolipoprotein E4 (APOE4) status, and education level, as these factors are known to contribute to the complex etiology of AD. The discovery dataset was defined as the ADNI risk dataset, and the demographic information of the 115 samples is detailed in Supplementary Table 1.

ii) OASIS dataset

To validate our model, we constructed two external validation datasets from the OASIS3 database: the OASIS risk dataset and the OASIS diagnosis dataset. To construct the OASIS risk dataset, we adopted the same participant selection process as that used in the ADNI risk dataset. The OASIS risk dataset consisted of 33 unique participants curated to reflect the characteristics of the ADNI risk dataset. In addition, the OASIS diagnosis dataset included 803 unique participants, each labeled with their respective diagnoses. More details on these datasets are provided in Supplementary Table 2 and Supplementary Table 3, respectively. Institutional Review Board (IRB) approvals were obtained from the original study depicting the ADNI and OASIS databases. In the ADNI dataset, consent forms were approved by each participating site’s IRB. The IRB of Washington University School of Medicine approved the OASIS study and written informed consent was obtained for all participants. All ADNI and OASIS datasets have been fully anonymized, with no protected health information included.

Data acquisitionI) MRI data

Imaging data were acquired using 3 T GE, Philips, and Siemens scanners. In ADNI-GO/2 phases, T1w images were obtained using a 3D magnetization-prepared rapid acquisition gradient echo (MPRAGE) sequence (repetition time [TR] = 6.8 ms; echo time [TE] = 3.2 ms, flip angle = 9°). The rs-fMRI data were acquired using a 2D echo planar imaging (EPI) sequence (TR = 3,000 ms; TE = 30 ms; flip angle = 80°; number of volumes = 140, and voxel size = 3.31 \(\times\) 3.31 \(\times\) 3.31 mm3). In ADNI-3 phase, T1w images were obtained using a 3D MPRAGE sequence (TR = 2,300 ms; TE = 2.98 ms; flip angle = 9°), and rs-fMRI data were acquired using an EPI sequence (TR = 3,000 ms; TE = 30 ms; flip angle = 90°; number of volumes = 197; and voxel size = 3.4 × 3.4 × 3.4 mm3). Imaging data in OASIS3 were acquired using Siemens TIM Trio 3 Tesla scanners. T1w images were obtained using a 3D MPRAGE sequence (TR = 2.4 ms; TE = 3.2 ms; flip angle = 8°), and the rs-fMRI data were acquired using a 2D EPI sequence (TR = 2,200 ms; TE = 27 ms; flip angle = 90°; number of volumes = 164, and voxel size = 4 \(\times\) 4 \(\times\) 4 mm3).

ii) PET data

For PET images from ADNI-GO/2/3, standardized dynamic protocols were applied. The amyloid-β PET consisted of AV45- and FBB-PET. For AV45-PET, a protocol lasting 50–70 min post-intravenous injection of 370 MBq of [18F] Florbetapir was used, with a scan duration of 20 min divided into four 5-min frames. For FBB-PET, the protocol lasted 90–110 min post-intravenous injection of 300 MBq of [18F] Florbetaben, with a 20-min scan divided into four 5-min frames. For FDG-PET, a 30–60 min protocol post-injection of 185 MBq of [18F] Fludeoxyglucose was used, with a 30-min scan divided into six 5-min frames. For Tau-PET, a 75–105 min protocol post-injection of 370 MBq of [18F] AV-1451 (Flortaucipir) was applied, with a 30-min scan divided into six 5-min frames. All PET data we employed were from baseline tracer scores matched to the baseline fMRI sessions. In the PET tracer comparison, only baseline session PET data were used. If PET data were unavailable due to missing modality, they were excluded from the PET tracer score analyses.

Data preprocessingI) MRI data

MRI data preprocessing of the ADNI and OASIS data was performed using fMRIprep [36]. T1w images were corrected for intensity non-uniformity. A T1w-reference map was computed after the registration of multiple T1w images, and the skull was removed. Brain tissue segmentation of the cerebrospinal fluid, white matter, and gray matter was performed on the brain-extracted T1w images using FSL. Volume-based spatial normalization to the standard MNI152 space was performed through nonlinear registration using ANTs [37]. Brain surfaces were reconstructed using FreeSurfer [38]. For each rs-fMRI session, slice timing was adjusted, and the head motion was corrected by registering the volume of each time point to the reference volume of the middle time point. The reference time-point data were coregistered to the T1w-reference image using boundary-based registration and then registered onto the MNI152 template. Spatial smoothing with an isotropic Gaussian kernel of 6 mm full-width half-maximum (FWHM) was applied, and removal of residual motion artifacts was performed using independent component analysis-based automatic removal of motion artifacts (ICA-AROMA) [39]. Finally, the non-steady first 10 s volumes were removed, and a bandpass filter was applied with a range of 0.008–0.1 Hz.

ii) PET data

The PET images were preprocessed as follows. Raw PET images were coregistered across different frames to reduce the motion effect, and then the frames over 5-min intervals were averaged. These images were reoriented to conform to a standard image grid of 160 × 160 × 96 matrix with 1.5 mm3 voxel size, and intensity normalization was performed. Spatial smoothing was applied with 8 mm FWHM. Further details are given on the ADNI website (https://adni.loni.usc.edu/data-samples/adni-data/neuroimaging/pet/). The PET images were coregistered to the corresponding T1w image with six degrees of freedom and subsequently nonlinearly registered onto the standard MNI space.

Imputation of time to conversion

The time taken to convert to AD can be considered the survival time commonly found in medical studies [34, 35, 40]. In clinical trials and longitudinal studies, a subject may experience a particular event, such as death or disease onset. The longitudinal data collected by the ADNI are suitable for survival analysis, as the study was conducted on the same subjects, and the time to conversion to AD can be utilized as a response variable. Identifying individuals at heightened risk of conversion to AD presents several research challenges, particularly the issue of precisely determining the timing of conversion to AD. Since the time to conversion to AD is interval-censored, we are certain that conversion occurred within the interval between the date of the last MCI diagnosis and that of the first AD diagnosis. The time to conversion to AD was estimated based on the imputed survival times. A Weibull accelerated failure time (AFT) model with clinical covariates was fitted using the icenReg R package [41]. After imputation, interval-censored data were converted into right-censored data so that conventional survival analyses could be applied.

Functional connectivity gradient estimation

We estimated cortex-wide FC gradients and subcortical-weighted manifold degrees from baseline fMRI data. We used the Schaefer atlas with 200 parcels for the cortex [42] and the Desikan-Killiany atlas for the subcortex [43]. First, FC matrices were generated by calculating Pearson’s correlations of the time-series signals between two distinct brain regions for each individual. Fisher’s r-to-z transformation was applied to the individual-level matrices. A group-level connectivity matrix for the ADNI samples was constructed by averaging the individual connectivity matrices. This group connectivity matrix was subsequently thresholded, retaining only the top 10% of elements per row. An affinity matrix was then computed using a normalized angle kernel, which helped capture the similarity in connectivity patterns between cortical regions. We estimated low-dimensional eigenvectors (i.e., gradients) from the affinity matrix using diffusion map embedding [44], which is a nonlinear manifold learning technique. The diffusion map embedding algorithm is robust to noise and provides better computational efficiency than other nonlinear manifold learning techniques [45, 46]. We followed the diffusion map parameter settings used in previous studies [30, 47]. After constructing the ADNI group-level gradients, we aligned them to the group-level gradients of an independent Human Connectome Project (HCP) database (https://github.com/CAMIN-neuro/caminopen/tree/master/gradient_align) [48], which was derived using the procedure for the ADNI dataset. This further alignment mitigated the bias of demographic attributes specific to the ADNI dataset and stabilized the manifold [49]. The computed group-level gradients were then used as templates for the individual-level gradient alignment. Alignment was accomplished using Procrustes rotation [50]. The entire gradient estimation process was performed using the BrainSpace toolbox (https://github.com/MICA-MNI/BrainSpace) [51]. For subcortical regions, we estimated subcortical-weighted manifolds through element-wise multiplication of subcortico-cortical functional connectivity and cortical gradients [29, 30]. The nodal degree was then computed, which effectively reflected the subcortico-cortical functional connectivity weighted by the cortical eigenvectors. We used the first two dominant cortical gradients, G1 and G2, along with the two corresponding subcortical-weighted manifold degrees, S1 and S2.

PET-SUVR estimation

PET scans enable the visualization of amyloid plaques and the assessment of glucose metabolism and tau spreading in the brain, offering information about the underlying pathological changes in AD. The standardized uptake value ratio (SUVR) for PET images was calculated using the cerebellum or brainstem as the reference region. Voxel values corresponding to the regions of interest (ROIs) defined by the atlases, specifically the Schaefer atlas [42] for cortical regions and the Desikan-Killiany atlas [43] for subcortical areas, were averaged and subsequently divided by the reference value. Following previous studies [52, 53], the reference regions were defined as the whole cerebellum for AV45, FBB, and tau-PET, and the brainstem for FDG-PET. These values within each ROI yielded regional SUVR values. Additionally, to facilitate the comparability of AV45-PET and FBB-PET measurements, we transformed their SUVR into centiloid (CL) units using the established transformation formula (https://adni.loni.usc.edu/wp-content/themes/freshnews-dev-v2/documents/pet/ADNI%20Centiloids%20Final.pdf). Finally, we calculated the global PET tracer score by averaging the regional SUVR values across all ROIs.

Prediction of conversion to AD

We applied a penalized Cox regression model with an elastic net penalty to analyze the time to conversion to AD, using FC manifolds as regressors. In the survival analysis framework, we have data in the form \((_, _,_)\), where \(_\) is the observed survival time and \(_\) is the indicator variable for the censoring of the \(i\) th sample. Further, we assume \(_<_<\cdots <_\) be the increasing unique survival times. The common Cox proportional hazards model [54] assumes a semiparametric form for the hazard:

$$h_i\left(t\right)=h_0\left(t\right)\exp\left(\boldsymbol\beta^T\boldsymbol x\right),$$

(1)

where \(_\left(t\right)\) is the baseline hazard at time \(t\), \(}=(_,_,\dots ,_)\) is the column vector of the regression predictors, and \(}=_,_,\dots ,_)}^\) is the vector consisting of the regression coefficients. The estimation of \(}\) is achieved by maximizing the partial likelihood,

$$L\left(\beta\right)=^m}\frac_}}}e^_}},$$

(2)

where \(_\) is the risk set at the time \(_\), consisting of indices \(j\) such that \(_\ge _\). The performance of the survival model was evaluated using the concordance index (C-index), which indicates how well the model predicts the ordering of times for a specific event. Because the estimated FC manifolds have high dimensionality, penalized models may be a more efficient approach. In penalized regression, the least absolute shrinkage and selection operator (LASSO) [55] method, characterized by an L1-penalty, along with the ridge [56] method using an L2-penalty, is typically employed. Ridge regression is particularly effective in addressing multicollinearity among features, whereas LASSO reduces dimensionality by reducing regression coefficients to zero. An elastic net [57] penalty combines the advantages of both LASSO and the ridge, providing an optimal solution. We estimated the \(}\) via the coordinate descent algorithm using the glmnet R package [58]. In our study, we denote \(}=(}}_,}}_)\) where \(}}_\) is the FC manifold vector, and \(}}_\) indicates the clinical variable vector. To control for the potential effects of the demographic and clinical variables \(}}_\) including age, sex, MMSE, APOE4 status, and education level, we considered these as additional covariates and imposed the penalty term only on the FC manifold \(}}_\). Accordingly, \(}}_=_,\dots ,_\right)}^\) and \(}}_=_,\dots ,_\right)}^\) are regression coefficients for the FC manifolds and the clinical variables, respectively. We defined a linear combination of \(}}^}\) as the risk score. In summary, we used FC manifolds (G1, S1, G2, and S2) and controlled for covariates as regression predictors. In our model, \(}\) can be estimated by maximizing the partially imposed penalized log partial likelihood \(l\left(}\right)=\text(L\left(}\right))\),

$$\widehat\beta=}_\left[l\left(\boldsymbol\beta\right)-\lambda P_\alpha\left(_\right)\right]=\left(_,_\right),$$

(3)

where

$$P_\alpha\left(\beta_\right)=\left(\alpha\Sigma\vert\beta_\vert+\frac12\left(1-\alpha\right)\Sigma\beta_^2\right),\;\left(0\leq\alpha\leq1\right).$$

(4)

The input feature scale is normalized using the mean and standard deviation. For model estimation, we implemented a five-fold cross-validation to select hyperparameters: L1-ratio \(\alpha\) and regularization parameter \(\lambda\). Hyperparameter selection was based on the highest C-index derived from the cross-validation procedure. This procedure was repeated 20 times, and the resulting outcomes were reported. For robustness, we averaged the estimated coefficients and risk scores from all optimally identified models. To assess the effect size related to the conversion to AD, we calculated the region-wise sum of the absolute coefficients of the FC manifolds (G1 and G2 or S1 and S2). This effect quantifies the extent to which an alteration in the functional connectivity of a region influences the likelihood of conversion to AD while controlling for other covariates.

Survival analyses of conversion to AD

We stratified participants at risk of conversion to AD into fast and slow conversion groups based on a zero cutoff of the predicted risk score and assessed between-group differences between the two groups using a log-rank test. To further examine whether PET scores were different between the two groups, we compared averaged global amyloid-β positivity (i.e., CL), glucose metabolism, and tau accumulation derived from PET scans between the groups using the Wilcoxon rank-sum test.

Considering the longitudinal nature of the data, the trajectory of the conversion to AD was traced for each participant. We identified four potential trajectories for participants at risk of conversion to AD: CN to MCI, MCI to censored MCI, CN (or via MCI) to AD, and MCI to AD. We denoted the progression paths from the diagnosis at the time of the baseline scan to either AD conversion or the last follow-up with an arrow (e.g., MCI → AD). We used a t-test to compare differences between the groups based on diagnosis at the time of the scan and an analysis of variance (ANOVA) with a post hoc analysis based on the Tukey honest significant difference test to compare differences between the conversion trajectories. Finally, the associations between the risk scores and both the clinical dementia rating sum of boxes (CDRSB) and Alzheimer’s Disease Assessment Scale (ADAS) were explored. The CDRSB is a widely used tool that quantifies the severity of dementia symptoms across multiple domains, while the ADAS provides a comprehensive assessment of the cognitive and non-cognitive functions of AD.

Robustness analyses and external validation

To demonstrate the robustness of our findings, we conducted sensitivity analyses using several methodologies. First, we considered the ridge penalty along with the elastic net penalty as an alternative to the penalized Cox regression model. Second, we separately fit the models using only the first (G1 and S1) and second (G2 and S2) manifolds. Third, while estimating the functional gradients, we applied different thresholds by retaining 5% and 15% of the elements per row in the FC matrix.

We validated our findings using an external validation dataset obtained from the OASIS3 database [33]. The replication involved two external validation sets: (1) the OASIS risk dataset (longitudinal) and (2) the OASIS diagnosis dataset (cross-sectional). (1) The OASIS risk dataset was created using the same participant selection criteria as that used in the ADNI risk dataset and was designed for external validation purposes. (2) The OASIS diagnosis dataset contained labels (CN, MCI, and AD) based on the corresponding diagnoses. We estimated the FC gradients from both datasets using the same procedure as that used for the ADNI risk dataset and aligned the individual gradients to the ADNI-to-HCP gradient template. Subsequently, the aligned FC gradients and subcortical-weighted manifold degrees were estimated. We then normalized the FC manifolds and demographic and clinical covariates using the mean and standard deviations derived from the ADNI risk dataset. We predicted the risk scores for the two external validation datasets from all identified models from 20 iterations and averaged the risk scores. For the OASIS risk dataset, we fit another Weibull AFT model using clinical covariates from the OASIS risk dataset and estimated the time to conversion to AD. We evaluated the C-index and stratified the participants into fast and slow conversion groups using a zero threshold. We then examined their conversion trajectories and risk scores. For the OASIS diagnosis dataset, we evaluated the risk score distribution for the CN, MCI, and AD groups and tested for significant differences among them. We also analyzed the association between the risk scores and CDRSB.