Fully bayesian longitudinal unsupervised learning for the assessment and visualization of AD heterogeneity and progression

Clustering evaluation

The reported results are based on 750 000 iterations with 500 iterations thinning where 250 000 iterations were the burn-in period, which therefore saved 1000 Markov Chain Monte Carlo (MCMC) samples. The distributions of the estimated parameters started converging after the burn-in samples and they remained stable thereafter for the remainder of the simulations. The deviance for the different models decreased with the increasing number of clusters (Supplementary Table 1). The different initializations produced various outputs from which the one with the packages’ default settings was the worst in terms of deviance. The model with initialization in the means of the clusters from our previous study [13] and the addition of uniform noise for eight clusters was optimal in terms of quality.

Figure 1 shows the multidimensional scaling coordinates of the component-subject probability matrix. Subjects are coloured dependent on the cluster to which they belong. Data from six subjects were excluded (Supplementary Figure 3 and Table 2): four subjects from the outlier clusters 7 and 8 based on the maximum probability rule (Figure 1A), and three subjects with uncertain classification with highest posterior density (HPD) intervals, analogous to confidence intervals in frequentist statistics (one subject from cluster 7 (already excluded) and 2 subjects from cluster 2) (Figure 1B). The remaining 66 subjects were used for further analysis. The separation between the six clusters in terms of probability for their subjects to belong to the same cluster is seen in Figure 1C, where clusters 1, 2 and 3 are clearly separated from each other. Visualization of the 1^st, 2^nd and 5^th multidimensional scaled (MDS) components shows the separation between clusters 4, 5 and 6 (Figure 1D).

Cluster characterization

Three main patterns of atrophy were found in the dataset: i) typical AD pattern (clusters diffuse 1, 2 and 3) (Figure 2B), ii) a minimal atrophy pattern (Figure 2A) and iii) a hippocampal sparing pattern (hippocampal sparing early and late onset) (Figure 2C).

The minimal atrophy cluster is characterized by initial atrophy in the entorhinal cortex and longitudinal thinning in adjacent inferior temporal gyrus (Figure 2A). The atrophy patterns in the three diffuse clusters (reported as typical AD) more closely follow the pattern of neurofibrillary tangles (NFT) spread as suggested by Braak and Braak [24]. However, differences between these three patterns do exist and may be attributed to age (even after correcting for this effect). The atrophy in the diffuse 3 cluster is more advanced (Figure 2B) and these subjects have lower cognitive performance (Supplementary Table 1). Two clusters are observed within the hippocampal sparing AD subtype. The degree of atrophy as well as the age at onset of dementia differentiate these two clusters (Figure 2C, Supplementary Table 1). The early onset hippocampal sparing cluster has a greater level of atrophy at baseline and accumulates atrophy faster over time, in contrast to the late onset hippocampal sparing cluster. In both hippocampal sparing clusters, the precuneus and the inferior parietal gyri are atrophied (Figure 2C). The 1^st and 3^rd quartile images show the dispersion around the mean cortical atrophy of each cluster (Supplementary Figure 1).

The six clusters did not differ in terms of sex distribution, but they did differ in years of education (Supplementary Table 1). The lowest and highest median years of education are observed in the diffuse 2 cluster (12 years) and the hippocampal sparing early onset (18 years). The two clusters with hippocampal sparing patterns of atrophy differ in several aspects, such as the age of onset of dementia. The minimal atrophy cluster has the slowest decline over time in the clinical dementia rating scale (CDR), while the hippocampal sparing early onset cluster has the steepest decline (Supplementary Table 1). The hippocampal sparing early onset cluster also has a steep decline in constructional praxis, and the greatest deficits in ideational praxis at baseline. Although the minimal atrophy cluster has the best scores in all the Alzheimer's Disease Assessment Scale (ADAS) subscales at baseline, the hippocampal sparing late onset group has a better score in the word recognition task at baseline but declines very fast during the next two years. Total Mini mental state examination (MMSE) scores also differed between clusters (Figure 3). The minimal atrophy cluster had the highest baseline MMSE scores and maintained them during the available assessments. The diffuse 3 group had the lowest baseline MMSE scores and a steep decline. However, the steepest decline in MMSE was observed in the hippocampal sparing early onset cluster that started at levels comparable to the minimal atrophy cluster.

Comparison to previous results

The subjects in the cross-sectional study [13] that were assigned to the diffuse 1 subtype are now distributed in more than one cluster with the highest prevalence in the diffuse 1 and 2 clusters (Table 1). Two subjects from the cross-sectional diffuse 2 cluster are now in the diffuse 3 and one in the outlier cluster 8. All the seven subjects from the cross-sectional hippocampal sparing subtype are still in the hippocampal sparing clusters. Three subjects assigned to the limbic predominant atrophy pattern in the cross-sectional study are now in the outlier cluster 7, diffuse 1 and the HPD uncertain group. The subjects in the minimal atrophy group are still mainly in minimal atrophy (17 subjects out of 20) while two subjects are assigned to the diffuse 1 cluster and one subject to the hippocampal sparing late onset cluster. Out of the four cerebrospinal fluid (CSF) Aβ1-42 negative AD subjects that are included in the current study, one subject is assigned to the longitudinal diffuse 2 cluster (was in the cross-sectional diffuse 1 cluster), one to the longitudinal outlier cluster 7 (was in the cross-sectional limbic predominant cluster) and two are assigned to minimal atrophy (both subjects were in the cross-sectional minimal atrophy cluster) (Table 2).

Participants

We used data from the Alzheimer’s disease neuroimaging initiative (ADNI), a project launched in 2004 in the US and Canada from Michael W. Weiner, MD. The initial goal of the ADNI-1 cohort was to gather neuroimaging data to better detect and track AD in its early stages. The inclusion criteria for AD patients were the following: 1) to fulfil the NINCDS/ADRDA probable AD criteria, 2) a CDR global score between 0.5 and 1, and 3) an MMSE total score between 20 and 26. The exclusion criteria for AD included: the use of psychotropic medication that could affect memory, history of significant head trauma, evidence of significant focal lesions at the screening MRI, and the existence of a significant neurological disease other than AD. For the CU subjects, inclusion criteria were an MMSE total score between 24 and 30 and a CDR global score equal to 0. Exclusion criteria for CU subjects comprised presence of depression, mild cognitive impairment (MCI) or dementia. For more information on the ADNI study, please see http://adni.loni.usc.edu/about/.

We included all subjects with longitudinal sMRI data and available CSF data (101 AD and 113 CU) from our previously published cross-sectional study [13]. In total, 75 subjects were excluded due to bad longitudinal image quality and processing results (see below). At baseline, 94% of the AD subjects were amyloid-beta (Aβ)1-42 positive, while only 31 CUs were included since we wanted them all to be negative for Aβ1-42 and phosphorylated tau (pTau). The cut-offs for Aβ1-42 and pTau used are discussed by Shaw et al. [33]. The CU sample was further limited by additional inclusion criteria: 1) remain as CU across all the available follow-ups and not only the follow-ups used in this study (0-36 months of continuous follow-up for the 31 CU subjects), 2) have longitudinal MRI for all the time points of the analysis.

Altogether, 104 individuals were included in the final analysis (Table 1), 72 AD patients (72 subjects had baseline and 12-month MRI scans, and 57 subjects had a 24-month MRI scan) and 31 CU (baseline, 12- and 24-month MRI scans).

MRI acquisition and preprocessing

The MRI dataset consists of high-resolution sagittal 3D 1.5T T1-weighted Magnetization Prepared RApid Gradient Echo (MPRAGE) volumes (voxel size 1.1×1.1×1.2 mm3). Full brain and skull coverage were required and detailed quality control (QC) was applied to all the images [34].

Images underwent pre-processing with the longitudinal stream of FreeSurfer 6.0, where a subject-specific template is used [35]. For this study we utilized cortical thickness values for 34 cortical regions (Desikan-Killiany [36] atlas) and 7 subcortical volumes (hippocampus, amygdala, putamen, caudate, thalamus, accumbens, pallidum) from each hemisphere (Supplementary Table 2). The estimated total intracranial volume (eTIV) was also extracted for the statistical modelling of the volumetric data [37]. All data was processed through theHiveDB system [38]. The FreeSurfer output underwent visual QC and 29 AD and 48 CU subjects were excluded due to low output quality or since less than two continuous time points existed per subject after the QC.

Statistical analysis

Data standardization

The cortical thickness and subcortical volume ROI data of AD patients were standardized based on the sample of CU subjects, including mean centring and unit variance scaling. In this study, we adapted this longitudinal data procedure in order to account for the atrophy that is caused by the normal ageing process in the CU group over time. This ensures that the ageing time interval is accounted for. The z-values were calculated using the following formula $z_{j,t}^i=x_{j,t}^i-{\widehat{\mu }}_{j,t}^{CU}\text{/}{\widehat{\sigma }}_{j,t}^{CU},$ where x is the original measurement of subject i, in the time point t for the region j, while $\widehat{\mu }$ and $\widehat{\sigma }$ are the mean and standard deviation of the CU group at time t and region j. After this calculation, each value will resemble an atrophy level corrected for normal ageing levels and also normal decline over time, which was not previously done, and is crucial for biological and clinical interpretation of brain atrophy.

Statistical longitudinal clustering

We used a generalized linear approach, which allows us to incorporate fixed and random effects that can serve in different ways in sMRI and other modalities. The algorithm clusters the random intercepts and slopes of each individual’s outcomes of interest (ROI measures in this study), with repeated measurements instead of repeated measurements data for each individual subject. A pair of subjects with similar estimated trajectories of atrophy (similar starting value/intercept and slope over time) are grouped together, while subjects with different trajectories are assigned to different groups. As previously discussed by Fraley et al. [39], if we know the number of clusters K, we can formulate the unobservable cluster allocation of subjects i as P(U_i = k; w) = w_k, k=1,..,K and i=1,..N. Here w is the vector of unknown cluster proportions that are positive and sum to 1. The meaning of U_i is that U_i = k when an observation Y_i is produced by the model density f_i,k(y_i, p_k, p), where p_k are cluster specific parameters and p are population parameters. The marginal density of Y_i is $\hspace{0.17em}{f}_i\left(y_i;\theta \right)={\displaystyle {\sum }_{k=1}^K{w}_k{f}_{i,k}(y_i,p_k,p)}$, where θ = {w^T, p_k, p} is the vector of unknown model parameters. Finally, clustering is based on the estimated values $\hspace{0.17em}{\widehat{p}}_{i,k}\hspace{0.17em}$ that resemble the chances of an individual i belonging in cluster k [39].

The conditional mean response for each region j in a regression can be expressed as

$E(Y_{i,j,l}\text{|}{\beta }_j,b_{i,j})=\hspace{0.17em}{x}_{i,j,l}^T\hspace{0.17em}{\beta }_j+z_{i,j,l}^T{b}_{i,l},\hspace{0.17em}l=1,\dots ,n_{i,j}$(1)

Here, x_(i,j,l) are vectors of known covariates (fixed effects), z_(i,j,l) is a vector that includes time values when the observations were made, β_r is the vector of unknown fixed effect for region j,b_j are i.i.d. random variables that express the j_th response of subject i. A Gaussian distribution is used to model the ROI data (general linear model), but data of ordinal or nominal nature can be analyzed by changing the link function on the left part of the equation (1). For each individual i the conditional distribution of the joint random effects vector B_i=(B^T_i,1,..,B^T_i,J) over j regions is multivariate normal,

$B_i\text{|}\hspace{0.17em}{U}_i=k\hspace{0.17em}~N({\mu }_k,D_k)$(2)

In (2), the i_th subjects belong to the k_th cluster (that is U_i=k, for a K cluster solution), μ_kis the unknown mean vector over j regions and D_k is the cluster-specific positive definite covariance matrix. For each individual i, Y_i,j,l (j=1,…,J and l=1,…,n_i,j) are conditionally independent given the random effects B_i. The random vectors Y₁,…,Y_N, as well as the random effect vectors B₁,…,B_N that were mentioned above, are independent. In summary, the μ_k and D_k, comprise the cluster-specific parameters that we will estimate given the data to explore the various grey matter atrophy patterns. The dependence among the Y_i,js (different markers j) is included in the non-diagonal components of the matrix D_k.

Accounting for external effects that might drive the resulting clusters within the model is convenient in this kind of analysis. Therefore, fixed effects β_j in (1), common to all clusters (population level effects) are estimated for each of the external variables and brain region that we want to assess during the clustering analysis. Adding the regression dispersion parameters ϕ_j to the β_j vector, we summarise the model parameters common to all clusters. Those parameters are also assigned distributions. For more information specific to the distributions of those effects and their hyperparameters, please see the supplementary material (model specification).

This fixed effect approach allows for the fitting of the resulting cluster profiles (atrophy maps) for different combinations of fixed effects to investigate their regional contribution. Since the longitudinal data from almost all cohorts with MR acquisitions typically have different numbers of visits per subject (irregularly sampled), we chose a model that can utilize all available measurements of each individual subject to calculate regression intercepts and slopes (vectors b_i,l). In the model, hierarchically centred generalized LMM are assumed with a non-zero unknown mean, b_l for l = 1, … L regions (see formula 1). The model that combines all the aforementioned features, i.e., Multivariate Mixture of Generalized Mixed effect Models (MMGLMM) [25], is applied to longitudinal trajectories of atrophy (see packages “mixAK” and “coda”, R version 3.0.0 or higher). We chose a Bayesian approach for the clustering optimization because it includes prior distributions in the parameter estimates. This enables the algorithm to investigate the parameter space even in cases of small samples where likelihood information is limited.

The clustering algorithm estimates different outcomes. One outcome is the different cluster components. Each estimated multivariate Gaussian component resembles a pattern of atrophy that is observed in the dataset. Each individual subject is assigned a probability to belong to any of the components (soft clustering), rather than being assigned to a single component. The assignment of subjects into clusters is based on the maximum posterior probability rule (an individual is assigned to the component with the highest individual component probability). This is a much more realistic approach in comparison to hard clustering approaches used in most previous data-driven studies [7, 8, 10, 17, 30, 40], since the heterogeneity in AD is modelled here as a continuum and allows for mixed patterns instead of single patterns. Hence, the data-driven algorithm provides explicit information on whether a subject has a distinct atrophy pattern or a mixture of patterns through the estimation of subject component probabilities. The proposed framework clusters subjects of a cohort into groups (provides probability of subjects to belong in any of the clusters) and not patterns of atrophy into groups for a cohort (clusters of regions/vertices), as in the study by Marinescu et al. [14].

A schematic representation of the proposed analytical framework is portrayed in Figure 4. The time from the first visit (baseline) was defined as a random effect for the sake of comparability with previous cross-sectional studies on AD subtypes [13, 17]. The intercept of the model will correspond to the atrophy levels on the first visit and the slope will show how these atrophy levels change over time. The fixed effects of the model are age, sex, education, years from the onset of dementia, total intracranial volume, and baseline CSF Aβ1-42 and pTau181 levels. The resulting clusters are visualized in terms of their fitted values on the median intercept (i.e. baseline), 12 months and 24 months after the baseline observation for a specific set of fixed effects. Only fitted values below two standard deviations of the CU mean are presented [41]. Measures of dispersion (1^st and 3^rd estimate distribution quartile) are also visualized in order to assess within-cluster variance. With those measures we can interpret how different the subjects within each cluster are. We also present the cortical maps of each individual and time point that was used in the analysis in the supplementary material to show how well the estimated components represent the individuals to which it is assigned.

The statistical model that we chose has all the features that were described above and its original specification can be found in [25]. The optimization was performed using the R language, version 3.4.1 [42]. The model is fully Bayesian and thus the output of the Markov chain Monte Carlo (MCMC) simulation is exploited to make inference on the population and cluster-specific parameters. To adequately explore the distributions of the estimated parameters and speed up convergence of the algorithm, we optimized the model from different initial values based on i) the packages’ default values [25], ii) previous study results [13] and iii) cross-sectional clustering on the baseline data including k-means clustering and hierarchical agglomerative clustering as well as the addition of uniform noise to increase randomness in the initialization [43, 44]. To identify the optimal solution for our dataset, we initially optimized models for 2-8 clusters for all the different initializations, summing to 49 MCMC chains. Then we assessed i) the model deviances (-2*logLikelihood) [25], ii) the quality of parameter convergence with respect to MCMC with high autocorrelation (visual inspection of the MCMC trace plots and auto-correlation values) [43] and iii) the quality of clustering with respect to observations with low classification certainty (See Supplementary Table 1). In our hybrid model evaluation approach, all three quality criteria were considered as important in the selection process (scaled to the same interval, 0-1) [45].

Statistical comparisons

Annual changes in cognitive measures for the different clusters were estimated with linear regression. For the MMSE total score, trajectories of decline were also calculated for the clusters with a mixed effect linear regression. Fixed effects were used to assess the differences at baseline and in decline between the clusters and random effects accounted for the repeated measurements. Significance of fixed effects was not assessed due to sample size and the interpretation was based on the coefficient standard errors. The above-mentioned hierarchical model statistical analyses (comparison between clusters) were carried out using R 3.6.3 software.