Multi-stage Biomarker Models for Progression Estimation in Alzheimer’s Disease

Multi-stage Biomarker Models for Progression Estimation in Alzheimer’s Disease Alexander Schmidt-Richberg1, Ricardo Guerrero1, Christian Ledig1, Helena Molina-Abril2, Alejandro F. Frangi2, Daniel Rueckert1, and on behalf of the Alzheimers Disease Neuroimaging Initiative 1Biomedical Image Analysis Group, Imperial College London, London, UK [email protected] 2Center for Computational Imaging and Simulation Technologies in Biomedicine (CISTIB), University of Sheffield, Sheffield, UK Abstract. The estimation of disease progression in Alzheimer’s disease (AD) based on a vector of quantitative biomarkers is of high interest to clinicians, patients, and biomedical researchers alike. In this work, quantile regression is employed to learn statistical models describing the evolution of such biomarkers. Two separate models are constructed using (1) subjects that progress from a cognitively normal (CN) stage to mild cognitive impairment (MCI) and (2) subjects that progress from MCI to AD during the observation window of a longitudinal study. These models are then automatically combined to develop a multi-stage disease progression model for the whole disease course. A probabilistic approach is derived to estimate the current disease progress (DP) and the disease progression rate (DPR) of a given individual by fitting any acquired biomarkers to these models. A particular strength of this method is that it is applicable even if individual biomarker measurements are missing for the subject. Employing cognitive scores and image-based biomarkers, the presented method is used to estimate DP and DPR for subjects from the Alzheimer’s Disease Neuroimaging Initiative (ADNI). Further, the potential use of these values as features for different classification tasks is demonstrated. For example, accuracy of 64 % is reached for CN vs. MCI vs. AD classification. 1 Introduction Alzheimer’s disease (AD) is a progressive neurodegenerative condition and the most common form of dementia. Patients that show first symptoms like general memory loss are usually diagnosed as suffering from Mild Cognitive Impairment (MCI), an early stage of dementia. Later, throughout the disease, these symptoms are followed by behavioural changes and further cognitive and functional decline. Patients become less able to perform simple tasks and increasingly depend on carers’ support. However, an objective staging of the disease is nontrivial due to the considerable variability in the patient’s age at disease onset and the individual rate of progression [3]. Estimating the current disease severity and the future rate of progression is of interest for patients and caregivers but also has potential to improve clinical trials as more homogeneous study groups can be recruited. Most existing approaches with this aim employ the results of cognitive tests like the Clinical Dementia Rating – Sum of Boxes (CDR-SB) as biomarkers to characterise disease progression. Many of those methods (see e.g. [5] for a detailed overview) use Markov transition models and Cox proportional hazard models to estimate disease progression. However, as detailed in [5], the main limitations of these methods are “the use of a limited number of health states to capture events related to disease progression over time”, and the fact that “a single symptom, such as cognition, is not able to characterise AD progression”. Therefore, additional meaningful biomarkers that better describe anatomical changes can be quantified from imaging data. For example, [4,12] use an event-based model to determine the order in which CSF-, image- and cognition-based biomarkers become abnormal and then employ this information to assign a subject to one of several discrete disease stages. As modelling disease progression by a number of discrete stages is a strong simplification, some approaches have been developed that acknowledge the course of disease as a continuous process. For example, Yang et al. [10] assume an exponential-shaped trajectory of the ADAS score. The authors then estimate atime shift γindicating the disease progress of a subject by fitting its ADAS scores to this curve. Similarly, Delor et al. [2] compute a disease onset time by adjusting subjects according to their CDR-SB score. The approach presented in this work builds upon these methods. Here, quantile regression is used to estimate typical trajectories of clinical biomarkers (see Sect. 2). In detail, two models are trained, one for the transition from CN to MCI and one for the MCI-to-AD conversion. These models are then combined to a multi-stage model for the whole course of the disease. Thereafter, a probabilistic model is derived that allows the estimation of a subjects current disease progress and rate of progression based on measured biomarker values. The approach is flexible with regard to the considered biomarkers, which can be based, for example, on cognitive scores, neuroimaging, or both. Moreover, missing measurements are handled in a natural way, this means, the approach can be employed even if the set of observed biomarkers is incomplete. The proposed disease progress estimation is evaluated in Sect. 3using clinical data. Its applicability for different classification tasks is demonstrated at the end of this section. 2 Methods To model disease progression, the existence of a set of biomarker values yb sv acquired from multiple subjects s∈S={1,...,n S}during multiple visits v∈Vs is assumed. Here, b∈Bsv denotes the index of the biomarker. Each biomarker vector is associated with the time tsv ∈Tof acquisition, measured in days after the first (baseline) visit, as well as the diagnosis dsv that was given during each visit. The number of visits can vary for each subject, Vs⊆V={1,...,n V}. Also, the biomarkers acquired at each visit might differ, such that Bsv ⊆B={1,...,n B}. In the training phase of the presented method, characteristic trajectories of biomarkers in the course of disease progression are learned based on a number of training subjects (Sect. 2.1). These models are then employed in the test phase to estimate how far and how fast test subjects have progressed along the disease trajectory (Sect. 2.2). 2.1 Model Learning Aim of the model training phase is to learn the temporal trajectory of biomarker evolution throughout the disease by determining the probability that a certain biomarker bhas a value ybat a specified time point. More technically, each measured biomarker value yb sv is understood as an observation of a response variable Ybat a disease progress (DP) psv ∈R(the explanatory variable or covariate). The conditional distribution of Ybgiven pis then denoted by fYb(y|p). Adisease progression model M(p) comprises the distributions of all biomarkers in Bon a domain P⊂R, such that M(p)={M1(p),...,MnB(p)}with Mb(p):=fYb(y|p)forp∈P. Another way of representing the model is by its q-quantile functions yb q(p), which can be derived directly from fYb(y|p)(for example, the median trajectory is denoted by yb 0.5(p)). The learning of a model consists of three main steps. First, the training subjects have to be temporally aligned to establish correspondences between the time points of observation. Progression models are then estimated using quantile regression to learn the probability distributions fYb(y|p). Since temporal alignment is based on the point of conversion from either CN to MCI or MCI to AD, two separate models are learned. These two models are then combined to a multi-stage progression model. Temporal Alignment of the Training Data. Temporal alignment aims at associating the time points tsv of biomarker acquisition to the corresponding DP psv. In detail, the goal is to find a strictly monotonically increasing time warp function τ(t) that maps the subject-specific acquisition time tsv ∈Tto the population-based disease progress psv ∈P, such that psv =τ(tsv). During model training, the time point t0 sat which the clinical diagnosis changes and thus indicates transition to a more severe disease state is set to p= 0, that means psv =τ(t0 s;tsv):=tsv −t0 s.(1) For this reason, a specific model is trained for each transition phase (here, CN-to-MCI and MCI-to-AD). To identify t0 s, the visits v∗and v∗∗ with the last CN (MCI) and the first MCI (AD) diagnosis are determined. The time point of conversion is then assumed to be the average of these two visits, i.e. t0 s:= 0.5·(tsv∗+tsv∗∗ ). CN/MCI model MCI/AD model Fig. 1. Approach for automatically determining the optimal model offset δp. Fig. 2. Example of the model composition. First, separate models are trained for CN-to-MCI and MCI-to-AD converters. An optimal offset between these models is then automatically determined and model training is repeated with the whole set of samples (Colour figure online). Learning Disease Progression Model. The conditional distributions fYb are learned independently for each biomarker using quantile regression via vector generalised additive models (VGAMs) [11]. In contrast to logistic or exponential regression [10], VGAMs do not depend on prior assumptions on the functional form for each predictor variable other than their smoothness. However, the domain Pof Mb(p) is limited to the progress interval contained in the sample set. This means Pis given by P=[p−,p +], with p−:= mins,v(psv) and p+:= maxs,v(psv) being the earliest and latest observed DP, respectively. Therefore, the models are extrapolated by a linear extension of the underlying predictor functions (see [11] for details) and P=Ris assumed in the following. Model Composition. To combine the CN/MCI and MCI/AD models, an optimal offset δpis determined by optimising the similarity of the models in the overlapping region. Given δp, the end point p[1] +of the CN/MCI model M[1] corresponds to p[1] +−δpin the MCI/AD model M[2]. Similarly, M[2](p[2] −) corresponds to M[1](δp+p[2] −)(cf.Fig.1). The quality of the fit is quantified by ˆ δp:= argmax δp 1 2M[1](p[1] +)−M [2](p[1] +−δp)+M[1](δp+p[2] −)−M [2](p[2] −) with M[1](p)−M [2](q)= 1 |B|b∈Bf[1] Yb(y|p)−f[2] Yb(y|q)dy being the area between the corresponding density functions (averaged over all biomarkers). After determining ˆ δp, the multi-stage model is retrained using the full set of samples with p= 0 defined as the point of conversion from CN to MCI (see Fig. 2). 2.2 Progress Estimation Once the disease progression model is built, the aim is to estimate the progress of any given subject. However, the point of conversion t0 sis usually unknown and thus Eq. (1) cannot be employed. Progress estimation is accomplished by finding the most likely time warp τ(t) that optimally fits the evolution of the biomarkers, as measured from the patient, into the progression model M. Let ts=(ts1,...,t snV)Tbe the vector containing the time points of all visits of subject sand τ(ts)=(τ(ts1),...,τ(tsnV))T. Let further ys=(ysv)v∈Vsbe the biomarker vector measured for s, with ysv =(yb sv)b∈Bsv denoting the values acquired at visit v. Based on ts, the most probable time warp ˆτsgiven ysis determined by maximising the logarithm of the likelihood function L(τ(ts)|ys). This means ˆτs:= argmax τ log L(τ(ts)|ys) = argmax τ log fY(ys|τ(ts)) (2) with Y=(Y1,...,YnB). The joint probability of all observations yb sv is then fY(ys|τ(ts)) =  v∈Vs fY(ysv |τ(tsv)) =  v∈Vs  b∈Bsv fYb(yb sv |τ(tsv)) . whereat all biomarker observations are assumed to be independent of each other. A simple translational time warp parameterisation is given by τ(p0;t):=p0+t. (3) Here, the disease progress (DP) p0∈Ris an offset that indicates how far the subject has progressed in the course of disease at the time point of the first visit. However, this simple model cannot accommodate for different rates of progression, which are known to exist between subjects [3]. If |Vs|>1, the extended affine time warp definition τ(p0,r;t):=p0+rt (4) can be employed, where r∈R+is a scaling factor indicating the disease progression rate (DPR). The optimal values ˆp0 sand ˆrsfor DP and DPR are determined by maximising Eq. (2) over all p0and r, i.e. ˆτs(·)ˆ=τ(ˆp0 s,ˆr;·). In general, more complex time warps definitions are possible. 3 Experiments and Results This section evaluates the presented approach using data from the Alzheimer’s Disease Neuroimaging Initiative (ADNI, Sect. 3.1). First, in Sect. 3.2, progression models are trained for all available biomarkers. These models are then employed to estimate DP and DPR. The applicability of these values for different classification tasks is shown in Sect. 3.3. 3.1 Materials For this study, all subjects enrolled in either ADNI1, ADNIGO or ADNI2 (up to an acquisition date of 30/01/2014) were considered. For model training, all subjects were selected that converted either from CN to MCI (60 subjects) or from MCI to AD (248 subjects). The number of available biomarkers can vary between the training subjects (see below). To obtain consistent results in all experiments, the test set consists of the subjects with a stable diagnosis for which all biomarkers are present at baseline (bl), month 12 (m12) and month 24 (m24) visit. In total, these are 158 (88 male, 70 female) subjects classified as cognitive normal (CN), 90 (50 male, 40 female) patients with early MCI (EMCI), 94 (63 male, 31 female) with late MCI (LMCI), and 88 (45 male, 43 female) patients diagnosed with AD. The set Bof biomarkers considered in this work consists of cognitive scores and image-based features as detailed in the following. Biomarker values are corrected for age using a linear regression on the baseline samples of all control subjects. Cognitive Scores Bcog:A number of cognitive tests are performed by subjects participating in the ADNI study at each visit. The achieved scores are used as biomarkers. These tests include the Minimental State Examination (MMSE), the Alzheimer’s Disease Assessment Scale (ADAS 11 and ADAS 13), the Functional Activities Questionnaire (FAQ), the Clinical Dementia Rating – Sum of Boxes (CDR-SB) and the Rey Auditory Verbal Learning Test (RAVLT). Not every test score is available for each subject and visit, such that the absolute number of available training samples varies between 393 to 399 samples from 60 subjects for M[1], and 1452 to 1480 samples from 248 subjects for M[2]. Volumes of Brain Structures Bvol:Further, the volumes of 35 distinct brain structures are used as biomarkers. For this, MR scans are first automatically segmented into 134 regions using the whole brain segmentation proposed by [8], which is based on multi-atlas label propagation with expectation-maximisation based refinement (MALPEM). Here, brain atlases from the MICCAI 2012 Grand Challenge on Multi-Atlas Labeling1are employed. Corresponding manual expert segmentations are provided by Neuromorphometrics, Inc.2under academic subscription. The 30 atlas segmentations are transformed to an unsegmented scan and fused into a consensus probabilistic segmentation estimate using a local weighting approach. Subsequently, all 134 probabilistic label estimates are further refined using image intensity information. To reduce the total number of models, left and right cortex are fused to single structures, resulting in 35 distinct anatomical regions. For procedural reasons, only segmentations for images acquired before 20/11/2013 were available, such that the total number of training samples for each structure is 219 (59 subjects) and 955 (247 subjects) for the two models. 1https://masi.vuse.vanderbilt.edu/workshop2012. 2http://Neuromorphometrics.com. Fig. 3. Examples for the learned multi-stage disease progression models (Colour figure online). Biomarkers Derived from Manifold Learning Bml:Features obtained from MR images using manifold learning (ML) have been shown to contain valuable information about disease severity and progression [6]. The main idea of ML is to find a meaningful, low-dimensional representation of a high-dimensional feature space, such that similar scans also have similar coordinates in the lowdimensional manifold. This is achieved in three steps. First, the image regions that are most relevant with regard to information about disease state are automatically learned using sparse regression as in [6]. To compensate for varying intensity values in the images caused by different scanners and acquisition protocols, local binary patterns (LBP) are computed in a 26-connected neighbourhood for voxels within these regions and used as features in the high-dimensional space. The manifold is then learned using Laplacian eigenmaps. The local geometry is determined via a sparse similarity graph, built using the sum of squared differences (SSD) as similarity measure. Connections in the graph are made between the knearest neighbours, with the additional constraint that an instance can only be connected to one instance per subject. The manifold coordinates are computed for all training subjects with at least 5 visits, which results in 185 samples from 41 subjects available to train the CN/MCI model and 859 samples (155 subjects) for the MCI/AD model. The manifold dimension is chosen to be d= 20, that means 20 features are obtained per subject per visit and denoted as D1 to D20. 3.2 Model Learning and Progress Estimation Progression models are trained for all 60 biomarkers. CDRSB and FAQ models could not be built for the CN-to-MCI transition because the majority of samples is clustered at y= 0. These biomarkers were therefore omitted for determining Fig. 4. Visualisation of the disease progress (DP) estimated with different biomarker settings (columns: number of visits; rows: biomarker sets). The distribution of DPs is indicated by the grey values. The x-axis gives the DP with solid and dashed lines at the points of conversion at p= 0 and p=δp, and the mean model range p−and p+. The red bars show the median and 25th/75th percentile of the estimated DPs (Color figure online). the model offset. Based on the remaining biomarkers, an ˆ δp= 1860 days is determined. Composed models are then retrained based on all samples and example results visualised in Fig. 3. Using the composed models, disease progression of all test subjects is estimated as proposed in Sect. 2.2. Plausibility of the estimated DPs is evaluated with regard to their ability do differentiate between the different diagnoses, that means to what extent an ordering of CN <EMCI < LMCI <AD is achieved. To this end, the disease progress ˆp0 sis estimated on the search space [−2500,4500] for all CN, EMCI, LMCI and AD subjects in the test set using several biomarker configurations. On the one hand, different sets of biomarkers Best are considered for estimation: Bcog,Bvol,Bml, as well as all imaging-based biomarkers Bimg := Bvol ∪Bml. On the other hand, biomarkers from one (baseline), two (baseline and m12) and three (baseline, m12 and m24) visits Vest are employed. The distribution of the estimated DPs depending on the diagnosis is visualised in Fig. 4. 3.3 Application: Classification Image-Based Classification of Disease Stage. One of the main research topics of image-based analysis of Alzheimer’s disease is the classification of subjects according to their diagnosis based on structural MR images. The high interest is highlighted by a classification challenge held in the course of MICCAI 20143. The estimated DPs are therefore employed as single features to distinguish between CN and AD, CN and MCI, MCI and AD, and all three classes simultaneously. Support Vector Machines (SVMs) are used as classifiers. Corresponding to Sect. 3.2, DPs are estimated using the translational time warp 3http://caddementia.grand-challenge.org. Table 1. Results for the image-based classification of subjects according to their diagnosis. Different imaging biomarkers acquired at a different number of visits are compared. For each test, the accuracies (ACC) of a 10-fold cross validation are given. Biomarkers Visits CN/AD CN/MCI MCI/AD CN/MCI/AD Best Vest ACC ACC ACC ACC Bml {bl}0.89 0.71 0.82 0.63 {bl,m12}0.91 0.71 0.82 0.63 {bl,m12,m24}0.91 0.69 0.82 0.61 Bvol {bl}0.83 0.68 0.80 0.59 {bl,m12}0.86 0.68 0.78 0.58 {bl,m12,m24}0.85 0.69 0.78 0.58 Bimg {bl}0.88 0.72 0.81 0.63 {bl,m12}0.90 0.71 0.83 0.64 {bl,m12,m24}0.90 0.71 0.83 0.64 definition (3) based on different biomarker settings and compared to each other. The cognitive scores are excluded because the diagnosis is made largely based on the CDR, such that including them biases the classification (e.g., classification accuracy reaches 1.0 for CN vs. AD and 0.96 for MCI vs. AD using Bcog at baseline). The results of a 10-fold cross validation are shown in Table 1. Classification of Stable and Progressive MCI. Of high clinical interest is the task of identifying subjects with memory complaints that will develop Alzheimer’s disease within a given period of time [9]. For example, this information is valuable to select subjects for clinical trials. Therefore, classification between MCI/AD converters and non-converters is performed using different features: A) DPs estimated with the time warp definition Eq. (3) using m12 measurements, B) DPs using baseline and m12 measurements (translational time warp Eq. (3)) and C) DPs and DPRs computed with the affine time warp Eq. (4) based on baseline and m12 visits. The set of converters (progressive MCI, pMCI) consists of all subjects that convert from MCI to AD between m12 and m36 visit. It is to be noted that these subjects are part of the training set used for learning the progression models, which introduces a bias in the test. However, adding one subject to a set of more than 1500 samples for the quantile regression barely changes the models, such that this effect can be neglected. Non-converters are test subjects with a stable MCI diagnosis (sMCI). In total, the data consists of 231 sMCI and 106 pMCI subjects. Weighted SVMs are used for classification to compensate for the unbalanced sets. Classification results are given in Table 2. Classification of Subjects with Rapid Cognitive Decline (RCD). The rate of cognitive decline is known to vary considerably between subjects. It is therefore of interest for family and researchers alike to predict if a subject