Predicting Four-Year’s Alzheimer’s Disease Onset Using Longitudinal Neurocognitive Tests and MRI Data Using Explainable Deep Convolutional Neural Networks

Source of data

Data from this study was obtained from the publicly available Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (http://adni.loni.usc.edu). The aggregated dataset derived from subjects recruited in ADNI1, ADNIGO, ADNI2, and ADNI3.

Participants and cohort definition

This study was composed of two separate training steps, each of which used a different cohort of data: 1) the first step used a cohort of single-timepoint cross-sectional data for classification task (classification step) using all patients that were diagnosed with either NC or dementia on their initial visit, and 2) the second step used a cohort of longitudinal data for prediction task (prediction step) using of all the MCI patients. For the prediction step, only patients with three visits, all six months apart, were considered. If a patient was diagnosed with MCI in the first visit and remained as MCI four years after the first visit, then they were classified as stable MCI (sMCI). Patients who were diagnosed with MCI on the first visit but were diagnosed with dementia four years later were classified as progressive MCI (pMCI). With this criteria, 797 patients with NC and 582 patients with dementia were selected for the classification cohort. To prevent data imbalance and bias during training, we randomly selected 500 patients from each group to ensure a balanced dataset. For the prediction cohort, 129 sMCI and 86 pMCI were selected to form a dataset of 215 patients. Table 1 summarizes these two cohorts along with their demographic information.

Outcomes

The primary outcome was conversion of MCI to clinical diagnosis of AD within 4 years from the time of the current visit.

Input variables

For both cohorts, structural T1-weighted 3D MRI scans were used as inputs to 3D convolutional neural networks (CNN). In the prediction cohort, we randomly selected one MRI scan from each of the three visits, and for the classification cohort we selected a single MRI scan from the first visit. To reduce variability in the dataset, all MRI scans were oriented and then linearly registered to the standard 1 mm MNI152 image template using FSL FLIRT. MRI scans were then skull-stripped using FSL BET to remove extraneous non-brain tissue. Raw volume intensities were rescaled into the range [0, 1] through min-max normalization.

In addition to the MRI scans, we also extracted the following clinical and cognitive tests for each patient: 1) Mini-Mental State Examination (MMSE), 2) Clinical Dementia Rating –Sum of Boxes (CDRSB), 3) Modified Preclinical Alzheimer Cognitive Composite with Trails test (mPACCtrailsB), 4) Modified Preclinical Alzheimer Cognitive Composite with Digit test (mPACCdigit), 5) Logical Memory Test –Delayed Total Recall (LDELTOTAL), 6) Alzheimer’s Disease Assessment Scale-Cognitive Subscale 11 (ADAS11), 7) Alzheimer’s Disease Assessment Scale-Cognitive Subscale 13 (ADAS13), 8) Rey’s Auditory Verbal Learning Test –forgetting (RAVLT_forgetting), 9) Rey’s Auditory Verbal Learning Test –immediate (RAVLT_immediate), 10) Rey’s Auditory Verbal Learning Test –learning (RAVLT_learning), 11) Functional Assessment Questionnaire (FAQ), 12) Trail Making Test-B (TRABSCOR), 13) ADAS-Cog Delayed Word Recall (Q4), and 14) APOE4. To minimize the effect of outliers during training, each variable was standardized to have a mean of 0 and standard deviation of 1 derived from the training set. We denoted these non-imaging variables as NI.

Model development and specification

Single-timepoint Alzheimer’s disease classification

A modified DenseNet-BC model with 3D MRI scans as an input was used classify NC or dementia. The densely connected convolutional blocks labeled “Dense Block” consist of a series of convolution operations stacked on top of each other. An overview of this architecture can be seen in Fig. 1. The DenseNet used for this study had four densely connected convolutional blocks of lengths 6, 12, 32, and 24 with the normal pooling transition layer in between the dense blocks. We set the compression factor between dense blocks to 0.5, the growth factor to 32, and removed any layer-wise dropout. There was a single dropout layer before the final linear layer with a drop probability of 0.5.

OPEN IN VIEWER

Fig. 1

DenseNet architecture. A) The DenseNet used in this study consists of 4 dense blocks (B) connected by convolution and pooling operations. The final prediction is done by pooling the output of the last DenseBlock and feeding that through a linear layer. B) The DenseBlocks used in the DenseNet architecture. It consists of several sequential convolutions, where the input to one convolution is the concatenation of all the previous outputs. The “+” operation in this figure denotes concatenation.

Multi-timepoints longitudinal disease progression prediction

We trained a smaller model that uses the DenseNet CNN described above to extract salient information from each MRI scan for predictions. The smaller model predicts if a patient will convert to dementia from MCI based on a feature embedding derived from the pretrained model and concatenated from multiple time points of data. To create the feature embedding, the MRI scans from each timepoint were fed into the pretrained DenseNet model. We then discarded the final classification layer, and only extracted the latent feature maps from the last convolutional layer after global average pooling and right before the classification layer, which is a one-dimensional tensor of 288 values. We also concatenate the neurocognitive scores of each timepoint as additional clinical variables to the respective tensor and then stack all three timepoints tensors together. The shape of the final feature embedding was (3, 302), since there were 14 clinical variables being concatenated for each of the three timepoints.

We then constructed a 1D CNN (Fig. 2) that uses the above mentioned multi-timepoints embedded and concatenated feature as an input and computes whether the given patient will convert from MCI to dementia within four years. This 1D CNN consisted of two 1D convolutions followed by 1D max pooling, and finally a linear layer to predict class logits.

OPEN IN VIEWER

Fig. 2

High-level overview of how the prediction model works. MRI scans for each timepoint are run through a pretrained CNN, and the output is combined with clinical variables for each timepoint. All three timepoints are then stacked on top of each other and then run through a smaller CNN to predict if a patient will convert from MCI to Dementia.

Training pipeline

For classification, data were split into an 80% training and 20% testing dataset. The model was optimized through a Stochastic Gradient Descent optimizer with a batch size of 5, learning rate of 0.01, weight decay of 0.001, and Nesterov momentum factor of 0.9. The model was implemented using the PyTorch Python library and all training was done on Google Compute Platform virtual instances with Tesla T4 GPU acceleration.

For prediction, a smaller model was trained to predict if a patient will convert from MCI to Dementia within four years using three timepoints of brain MRI scans and cognitive tests. After training the classification model, the weights were frozen and used to train a smaller model that predicts whether a patient will convert from MCI to Dementia within 4 years using multiple timepoints of input MRI scans and neurocognitive scores. To prevent data leakage between cohorts, the prediction cohort only consisted of patients that were diagnosed with MCI on their first visit. Nested cross-validation was used to train and optimize the model while making sure it is robust enough to perform well on new data. Figure 3 shows how the cross-validation folds were acquired. The prediction cohort was split into five outer folds. In each training fold, four of these folds are combined to form the training set, and the fifth fold is used as an independent testing set. Each training set is split into 3 inner folds. Two of these folds were used to train the model, and the third fold is used to evaluate model performance and tune hyperparameters. After optimal hyperparameters were selected in the inner fold, the model was then retrained on the entire training set and evaluated on the independent testing set. To prevent any data imbalances, all training data splits were stratified, so there were an equal number of converters and non-converters within each subset. To comprehensively evaluate the model performance, the balanced accuracies, sensitivities, and specificities of the trained model on each validation and test set were recorded. Because the feature embeddings were much smaller than the MRI scans, the batch size was increased to 8 and the learning rate was decreased to 0.0001, and the other optimization parameters remained the same.

OPEN IN VIEWER

Fig. 3

Training, testing, and validation splits acquired through nested cross-validation. The entire dataset is split into five folds; four of them are used as a training set and the fifth is used as an independent testing set to evaluate model performance. The training set is split into three inner folds, and two of them are used to train the model while the last set is used for hyperparameter tuning.

Note that, the phase one training was done using only the CN and AD subject, while the phase two training was done using the MCI patient only, and there was no “data leakage” from the validation data from the phase 1 training to the phase two training. To evaluate the effect of combining multiple data modalities, ablation studies were also designed. One alternative model was also trained separately using only MRI images, and another model was trained using only NI variables. To evaluate the effect of including multiple timepoints of data, all the experiments were run using a single timepoint of patient data, and then again with three timepoints of patient data.

Saliency-map-based model explainability

After training, the best-performing model was re-evaluated on the testing set and heatmaps were generated based on where the model focused on the most when making predictions. To visualize the most salient brain regions, the Grad-CAM technique was modified to work in 3 dimensions. Since the model reduced the resolution of the MRI volume to 2x3x2 voxels, the 3D Grad-CAM technique was applied to the output of the first convolutional block, which had an output resolution of 32x32x32 to obtain a better visualization of the heatmap. Finally, the heatmaps were averaged over all 3 timepoints. This approach enabled visual highlighting of the sections of the images that were most significant to the network.

Model performance

Figure 4 shows the balanced accuracy of the classification model for each epoch as it trained. The CNN achieved an optimal accuracy of 0.88 around epoch 100 and had a loss of 0.38. After this point, the validation loss and accuracy began to level off and gains in performance decreased with increasing epochs.

OPEN IN VIEWER

Fig. 4

Training accuracy and loss curves for classification model. Data is split 80/20 into a training and testing set, and then used to train the DenseNet model to classify 3D MRI scans as having NC or Dementia. The model overfits as it approaches epoch 100 and achieves an optimal testing accuracy of 88%.

Figure 5 shows the average balanced accuracy of the prediction model for each epoch as it trained using 4 timepoints of input data. This model achieved its optimal average accuracy of 0.83 and optimal loss of 0.41 around epoch 500, at which point the model performance plateaued. Multiple experiments were run to compare how the model performance changes based on the input modality, and the validation and testing performance from each of these experiments for one time point and three time points are recorded in Table 2. For single timepoint experiments, the NI+MRI model outperformed both the single data modality models with a testing balanced accuracy of 0.8075. Similarly, the three timepoint NI+MRI model outperformed both the single data modality models with a testing balanced accuracy of 0.8347. In all experiments, the three timepoint models outperformed the one and two timepoint models and the NI+MRI models outperformed the individual NI and MRI models. This accuracy is an average of the optimal accuracy in each testing fold, rather than the average accuracy at a single epoch. In the experiments that achieved a balanced accuracy of 0.8347, the optimal accuracy in each testing fold occurred at epochs 127, 148, 429, 456, and 510.

OPEN IN VIEWER

Fig. 5

Average training and validation accuracy of prediction model over all cross-validation folds for each epoch. The darker line represents the average accuracy over all folds during training, and the lighter shaded region is the lower and upper accuracy bounds encountered during training. The green curve corresponds to the performance on the training subset, and the red curve corresponds to the performance on the validation subset.

To visualize the important regions that the model focused on when making the prediction, heatmaps were generated using the 3D Grad-CAM algorithm. The average heatmap over all patients within the testing set is shown in Fig. 6. The regions on MRI that were important for prediction were putamen, thalamus, amygdala, frontal pole, frontal gyrus, and planum polare.

OPEN IN VIEWER

Fig. 6

Heatmap visualization for a single patient in the prediction task. The model predicted if the patient would convert from MCI to Dementia and used the Grad-CAM algorithm to compute a heatmap of the salient features it looked at from the MRI scan. This heatmap was superimposed onto the original MRI scan and then averaged over all three timepoints. Areas highlighted in red were most important in the prediction, and each row represents a different slice in the different anatomical planes.

Transfer learning

The neural network used for MRI feature extraction in this study was first trained to classify between normal cognition and dementia using MRI scans. This transfer learning approach has been found to be more effective than training a single large model on the final prediction task because it requires much less data for each cohort and shrinks the parameter space for the final task, resulting in faster training times [19]. However, this also means the performance of the prediction model is dependent on the quality of features extracted by the classification model. Earlier studies were able to train 3D CNNs with an optimal classification accuracy of 82% [24] and 86% [13] and use them as feature extractors, and thus our classification model was a good candidate for transfer learning.

Multiple time points

To make predictions that factor in the temporal relationship between multiple timepoints of feature embeddings, we trained a 1D CNN, which achieved an improved performance in predicting risk of dementia onset in four years. This novel approach decoupled the computationally intensive 3D spatial MRI analysis from the smaller model that focuses on the temporal relationship within the data. The 1D convolutions also allowed us to combine the MRI feature embeddings and NI values and use them for prediction without any bias towards one type of data or the other. To our knowledge, there are no other studies that use a 1D CNN for longitudinal analysis of dementia. In the current study, we choose to use 1D CNN that can take the concatenated feature embeddings as multi-channel input, thereby effectively captureing the temporal information from multi-timepoint. Alternative approaches can also be considered to capture the temporal information from longitudinal data across multiple timepoint, such as RNN, specifically LSTM [30–32]. These RNN-based approach would be more effective to learn longitudinal patterns for longer sequence data with more number of timepoints or with varing length of time sequence. On the other hand, the research question that the current study aims to tackle focuses on utilizing a fixed longitudinal time window of four years, for which a 1D CNN is suitable to capture the inter-timepoint through the convolution kernal.

Heatmaps

Deep learning models are often seen as “black boxes” that can be difficult to interpret. This is especially problematic in medical image diagnosis, where it is important to understand why a particular diagnosis was made. While efforts have been made to develop methods to explain the decision-making process of deep learning models, the explainability of these models is still limited.

Heatmaps enabled visualization the distribution of the attentions of the model across different functional and anatomical regions in the brain. Such visualizations reflect brain regions that were most relevant to ML algorithms to predict MCI and AD conversion, and could potentially improve the explainability of deep-learning models. In our study, the most salient brain structures shown on the heatmaps were the putamen, thalamus, amygdala, frontal pole, frontal gyrus, and planum polare. Although the brain structures implicated in the progression from MCI and AD in the literature varied depending on dataset, variables included in the models, prediction models used among others, most of these regions are known to be frequently affected by neurodegenerative diseases like dementia. A meta-analysis by Zakzanis et al. found that the temporal lobes, amygdala, thalamus, and hippocampi exhibit the most significant differences when comparing normal brains to those with dementia [33]. Other studies using the ADNI dataset have also shown that the gyri, anterior hippocampus, and amygdala tend to have more pronounced deterioration within the pMCI cohort, which supports the results from our heatmaps [34]. Some studies that applied CNNs for classification in the ADNI dataset also found the same structures to be important in heatmaps generated through gradient and occlusion-based algorithms [35]. Some regions like the frontal gyrus and planum polare were consistently highlighted by our neural network but were not found to be subject to significant deterioration in other studies concerning dementia. Other brain regions that have been shown to be associated with development of AD, such as the default node network and ventricles, were not consistently highlighted in the heatmaps.Differences in findings could be due to different dementia subtypes, heterogeneity of various cohorts, sample sizes, data types included, and methodologies, among others.

Prior studies

Previous studies have shown that models that used multiple timepoints or multiple modality data modalities improve prediction of MCI to dementia conversion compared to those using a single timepoint or a single data modality. Grassi et al. trained several models to predict the onset of dementia from patient demographic information and clinical test scores and found that these models reached an optimal balanced accuracy of 79% [22]. Similarly, Zhang et al. also obtained an accuracy of 78.79% for this task, but they trained a 3D CNN that used structural MRI scans as an input rather than using non-imaging data [21]. Both studies used only a single data modality as inputs to their models. Other studies by Spasov et al. and Huang et al. have explored combining MRI data with non-imaging data and PET data to train more comprehensive models and achieved accuracies of 86% and 76% respectively [23, 25]. These studies were successful in combining multiple types of imaging and non-imaging data to make predictions but did not explore how the change in this data over time can aid prediction. Our study built on prior studies to include longitudinal MRI and non-imaging data to predict the onset of dementia.

Previous studies that have utilized longitudinal data to train models for this task have achieved performance levels like those of models based on single timepoints and have also produced useful visualizations of important brain regions that are involved in prediction. Ocasio et al. evaluated different CNN approaches and found the sequential convolutional approach yielded slightly better performance than the residual-based architecture, the zero-shot transfer learning approach yielded better performance than fine tuning, and CNN using longitudinal data performed better than CNN using a single timepoint MRI in predicting MCI conversion to AD [18]. The best CNN model for predicting MCI conversion to AD at three years yielded a balanced accuracy of 0.793. Heatmaps of the prediction model showed regions most relevant to the network including the lateral ventricles, periventricular white matter and cortical gray matter. Our study built on Ocasio et al.’s study by replacing the Siamese network architecture with a two-phase training methodology, which allowed us to incorporate more than two timesteps of patient data. Additionally, we adopted a much larger residual CNN architecture and incorporated non-imaging data into the final model.

Limitations and future directions

The sample size of MRI data was small compared to feature space that could result in overfitting, although the two-stage framework helped train the embedding model on the limited data set. A drawback of the two-stage framework was that the prediction model learned the temporal relationships between lossy feature embeddings rather than full MRI scans. Future studies may investigate an end-to-end approach to predicting conversion from MCI to AD where the model can learn temporal relationships using the entire MRI scan. Large and diverse datasets are needed to improve generalizability.

This study only explored a single type of CNN architecture and evaluated one way of creating a feature extractor model. This model does not need to be created using a classification cohort consisting of normal/dementia patients, so future studies should explore other ways of pretraining a feature extracting network such as convolutional autoencoders. Additionally, we only used a subset of the available NI data in ADNI. With more advanced data filling and analysis, it may be possible to include more descriptive clinical variables or cognitive tests.

The study used only anatomical MRI data. Multiparametric MRI (such as diffusion-tensor imaging, task functional MRI and resting-state MRI) will be incorporated into these models in the future. Similarly, other modalities such at PET and additional non-imaging clinical data can also be included in the model [20]. Other DL methods, such as deep survival analysis [36, 37] should also be explored to the prediction of MCI conversion to AD.

Conclusions

We developed a deep-learning model to predict MCI to AD conversion within 4 years using longitudinal whole-brain 3D MRIs and along with neurocognitive test scores. This framework sets the stage for further studies of additional data time points, different image types, and non-image data to further improve prediction accuracy of MCI to AD conversion. Accurate prognosis could lead to better management of the diseases, thereby improving the quality of life.