Detection of Alzheimer’s disease by displacement field and machine learning

Introduction

Alzheimer’s disease (AD) is an abnormal process of aging. AD is a special category of senile dementia (SD) which leads to short-term and long-term memory, thinking, and behavior (Dong et al., 2015a). Research on AD has attracted scholars from all over the world because of its importance and effect on the society. Symptoms of AD may become severe enough to interfere with daily life, and to death (Goh et al., 2014; Hahn et al., 2013). There is neither cure nor treatment for AD. In 2006, more than twenty million people in the world suffered from this disease (Murphy et al., 2011). AD is predicted to influence 1 in every 85 people worldwide about thirty years later, and more than forty percent of prevalent cases need high level of care (Brookmeyer et al., 2007).

Now that the earth is growing into an aging society, AD has caused heavier burdens on families and society than before (Esposito et al., 2013; Kantanen et al., 2014). In China, AD accounts for more than half of SD, which extracts a total economic loss of more than eighty billion yuan every year, and is responsible for nearly sixty billion yuan in healthcare costs every year (Song & Wang, 2010). In United States, healthcare on people with Alzheimer’s disease currently costs roughly $100 billion per year and is predicted to cost $1 trillion per year by 2050 (Miller, Erlien & Piersol, 2012; Yuan, Machado & Arias-Carrion, 2014).

Nowadays, it is beneficial to develop early and accurate detection methods for AD, which is also necessary for the treatment and management for controlling the deterioration of AD (Zhang et al., 2015b). A 3D scan of the whole brain becomes acceptable and affordable with recent advances in neuroimaging technology (Zhang et al., 2015e; Zhang et al., 2014b), especially by the help of the most popular imaging technique: magnetic resonance imaging (MRI). With its high-resolution magnetic resonance (MR) images, the diagnostic accuracy of AD are greatly enhanced. MR images already play an essential role in detecting AD from normal elder controls (NC).

In this study, we proposed to employ displacement field (DF) to track the morphometry from normal brains to AD brains. The main advantage is that this method does not need to segment region-of-interest (ROI) beforehand, namely, it obtains the DF for the whole image. The classical level-set method was employed to estimate the DF.

This work is structured as follows: ‘State-of-the-Art’ provides the state-of-the-art. ‘Materials and Methods’ presents the materials and methodology. ‘Experiments and Results’ gives the experimental results. ‘Discussion’ discusses the results and the proposed method. ‘Conclusions’ concludes this work. The nomenclature in this work is listed in Table 7.

State-of-the-Art

In the past, most diagnosis work was carried out by measuring region of interest (ROI) of brain MR images, since researchers already know several typical AD-related regions and corresponding shape deformation (Anagnostopoulos et al., 2013; Kubota, Ushijima & Nishimura, 2006), such as the enlarged ventricles, the shrinkage of hippocampus (HC), and the shrinkage of cortex (Pennanen et al., 2004). Somehow, the ROI-based methods suffer from some shortcomings: (i) The ROI methods need a priori information and expert knowledge. (ii) The detection accuracy relies on the experiences of the interpreters (Yang et al., 2015). (iii) The mutual information among the voxels is difficult to implement (Lee, Park & Han, 2013; Xinyun, Wenlu & Xudong, 2011). (iv) We need to explore other potential regions that may be connected to AD (El-Dahshan et al., 2014). (v) Automatic segmentation of ROI is not feasible in practice, and examiners tend to segment the brain manually (Zhang et al., 2013).

Recently, a new type of methods, the “whole-brain analysis”, gets popularity since it considers all voxels in the brain as a whole. It does not need to segment the brain beforehand, and it does not need any biomarker for classification task. The main disadvantage is the curse of dimensionality, which can be solved as the high-speed computers are very cheap nowadays (Álvarez et al., 2009). The whole-brain analysis heavily relies on pure computation, and it can be accomplished by only computer scientists, after physicians help to label the data as either AD or healthy. Generally, the whole-brain analysis treats the whole brain as a ROI, and it consists of two stages: feature extraction and classification. We reviewed more than 10 literatures, and analyzed them in detail.

Scholars have presented various methods to extract efficient features for AD and other types of pathological brain detection.¹ In addition, various classification models and methods exist, nevertheless, only a few of them are suitable for MR images. El-Dahshan, Hosny & Salem (2010) extracted the approximation and detail coefficients of 3-level discrete wavelet transform (DWT). They used artificial neural network (ANN) and K-nearest neighbors (KNN) classifiers. Plant et al. (2010) used brain region cluster (BRC). They suggested to use information gain (IG) to evaluate the interestingness of a voxel, and applied clustering algorithm to identify groups of adjacent voxels with a high discriminatory power. They used support vector machine, Bayes statistics, and voting feature intervals (VFI), with the aim of pattern classification. Park (2012) employed manifold learning for classification. Their use of the first and second distance measure resulted in an 18% and 46% error rate for classifying between patients with AD and normal patients, respectively. Chaves et al. (2012) utilized large margin-based methodology for AD detection in SPECT and PET images. Their system yielded accuracy, sensitivity and specificity values of 90.67%, 88% and 93.33% (for PET) and 92.78%, 91.07% and 95.12% (for SPECT), respectively. Saritha, Joseph & Mathew (2013) was the first to use wavelet-entropy in pathological brain detection. They employed spider-web plots to reduce the size of feature space. They also used the probabilistic NN (PNN) for classification. Zhang et al. (2015a) suggested that removing spider-web-plot yielded the same classification performance. Savio & Grana (2013) presented a novel method that employed deformation-based morphometry (DBM) method. They tested five features, and found three features performed excellent as trace of Jacobian matrix (TJM), modulated GM (MGM), and Geodesic Anisotropy (GEODAN). Furthermore, they utilized Pearson’s correlation (PEC), Welch’s t-test (WTT), and Bhattacharyya distance (BD), to measure the significance of voxel site. Kalbkhani, Shayesteh & Zali-Vargahan (2013) modeled the detail coefficients of 2-level DWT by generalizing autoregressive conditional heteroscedasticity (GARCH) statistical model, and the parameters of GARCH model were considered as the primary feature vector. They tested the KNN and SVM models. Wang et al. (2014) proposed to utilize the undersampling (US) technique on a three-dimensional image. They used singular value decomposition (SVD) to select features. Finally, they combined kernel SVM (KSVM) with decision tree (KSVM-DT). Zhou et al. (2015) followed Saritha’s method, and again employed wavelet-entropy for feature extraction. Naive Bayes classifier (NBC) was employed for abnormal brain detection. Harikumar & Kumar (2015) analyzed the performance of ANN, in terms of classification of medical images, using wavelets as feature extractor. Their classification accuracy achieved 96%. Zhang et al. (2015c) employed discrete wavelet packet transform (DWPT). They utilized Tsallis entropy to get wavelet packet entropy features. They introduced a generalized eigenvalue proximal support vector machine (GEPSVM). Nazir, Wahid & Khan (2015) suggested to use filters for the removal of noises, and extracted color moments as mean features. Finally, they achieved an overall accuracy of 91.8%. Dong et al. (2015b) employed the eigen-brain (EB) to extract features, and then they employed Welch’s t-Test (WTT) for the aim of reducing features. They proposed to use SVM with radial basis function (RBF). Damodharan & Raghavan (2015) combined tissue segmentation and neural network for brain tumor detection. Zhang et al. (2015f) proposed a novel classification system that implemented 3D discrete wavelet transform (3D-DWT) to extract wavelet coefficients the volumetric image.

After reviewing the latest literature, two findings can be obtained: (1) The DWT based features are efficient. However, we presented a novel feature of displacement field, which is the first to be used in MR images. (2) SVMs were commonly used, compared with conventional decision tree, artificial neural network (Wang et al., 2015a; Zhang et al., 2015d), and other classifiers (Collins & Pape, 2011). Hence, we continue to use SVM. Besides, two variants of SVM are introduced in this study: generalized eigenvalue proximal SVM (GEPSVM) and twin SVM (TSVM), with the aim of augmenting the classification performance further.

Materials and Methods

Co-registration and brain-masking

All three-dimensional MR brain images of each subject were motion-corrected, and coregistered to generate an averaged three-dimensional image, with the aim of increasing signal-to-noise ratio. Afterwards, those three-dimensional images were spatially co-registered to the Talairach space and then were brain-masked.

Figure 1 offers an example of the preprocessing on 3D images with resolution of 1 mm × 1 mm × 1.25 mm. The motion-correction procedure registered the images of three scans, and then generated an image in the original acquisition space with resampling to 1 mm × 1 mm × 1 mm. Finally, the averaged image was normalized to the Talairach coordinate space, and brain-extracted (Fig. 1).

Key slice selection

Calculating the displacement field on the whole brain was time-consuming. Therefore, we introduced the “key-slice (KS) selection” method that picks up key slices containing structures indicative of AD from NC. The procedure was as follows: we established a criterion called “inter-class variance (ICV)” v as (1)${\displaystyle v\left(k\right)=\Vert {\mu }_{\text{A}}\left(k\right)-{\mu }_{\text{N}}\left(k\right){\Vert }^2}$ where k is the index of key slice, ‖ ⋅ ‖² represents the l₂-norm, μ_A and μ_N represents the mean of gray-level values of the kth slice of AD subjects and NC subjects, respectively, and v denotes the ICV. We select the KSs whose ICVs are larger than half of maximum ICV, with 10× undersampling factor (i.e., every ten slices).

Besides, the direction of key slice can be either sagittal, or axial, or coronal. Scholars found images along coronal direction give a well-defined view compared to those along sagittal and axial directions. Figure 2 shows an image along coronal direction has an advantage because it covers 3 most important AD-related regions, which are treated as indicative of AD. These AD-related regions consist of the cerebral cortex (CC), the HC, and the ventricle. When sagittal or axial directions are used, two or even more slices are needed to cover those tissues. Therefore, we chose the coronal direction for key slice selection, in order to use less slices.

Shape registration

The shape registration consists of both global and local registrations. The global registration estimates the rigid parameters, which was accomplished in ‘Co-registration and brain-masking’. Then, the local registration finds the “displacement field (DF)” between the moving and reference images (see Fig. 3).

The rigid registration is an essential preprocessing procedure before non-rigid registration. It can fix the deformation stemming from the position, move, and pose of patients. Then, the following non-rigid registration can reflect the shape deformation of diseases.

The problem of finding the non-rigid estimation can be thought of a motion estimation task between a normal brain (moving) and an AD brain (reference). Several types of solutions are available to solve this task, such as spline function based methods, phase-correlation methods, fluid methods, optical-flow methods, elastic methods, etc. Among the above methods, the first type is parametric, hence, this kind of solution needs to solve optimal spline-based function parameters (Petibon et al., 2013). The second type needs a mass of computation resources; it is difficult to determine the local search range, and it is impossible to guarantee to find global optimal points (Lu, Yang & Zhang, 2014). The other three methods are non-parametric. They find the DF by solving directly a predefined physical model using partial differential equation (PDE) (Yang et al., 2007).

The level-set motion estimation method (Huang et al., 2014; Lee, Sandhu & Tannenbaum, 2013; Vandemeulebroucke et al., 2012) is a rather novel method, which is formed on the basis of the level set evolution theory. The moving image I₁ morphs iteratively along its gradient direction, till it deforms close to the given reference image I₂. The displacement field is in the form of (2)${\displaystyle \frac{\text{d}V}{\text{d}t}=\left(I_2-I_1\left(V\right)\right)\frac{\nabla I_1\left(V\right)}{\left|\nabla I_1\left(V\right)\right|}}$ where V represents the displacement field, and I₁(V) the deformed image of I₁ by V. The above equation can be solved by iterative algorithms described in reference (Kodipaka et al., 2007).

Figure 4 gives an example of the displacement field between a glioma brain and a normal one. Figure 4A presents a moving image while Fig. 4B is a reference image. Their overlap is shown in Fig. 4C. Apparently, the reference image is larger than the moving image, and they are not quite a match. Figure 4D shows the result of rigid-registration, and Fig. 4E shows the two images now match well except the cortex, the temporal and occipital lobe. Figure 4F shows the result of non-rigid registration. Compare to Fig. 4B, they closely match each other, and Fig. 4G also suggests a nearly-perfect match. The displacement field is shown in Fig. 4H, with two enlarged areas in Figs. 4I and 4J.

The values of displacement field are complex, i.e., the real part represents the displacement field along the horizontal direction, and the imaginary part along the vertical direction. Those values of Cartesian system are transformed into polar coordinate system, viz., the magnitude-direction spaces. Hence, the inputs to the classifiers are two fold in this study: (1) the directions of the displacement vector of the whole brain and (2) the magnitude of the displacement field of the whole brain. We did not use the sign of the displacement vector, which were commonly used in literature (Das Gupta & Tamadapu, 2013; Kodipaka et al., 2007). The reason is the direction field already contains the information of sign of the displacement.

Statistical analysis

Generalization error was obtained by K-fold CV, and K was set to 10 because of two reasons: (1) to balance between computational cost and reliable estimates, and (2) for fair comparison since the common convention is to set K with the value of 10 (Zhang et al., 2014a).

For a 10-fold CV, the dataset is randomly divided into 10 mutually exclusively folds of nearly equal size. In each run, 9 subsets are used for training, and the rest for validation (see Fig. 5). Above procedure repeats 10 runs, such that each subset is used once for validation. The 10 results over validation set are combined together along the diagonal blocks in Fig. 5, with the aim of producing an individual out-of-sample evaluation. The 10-fold cross validation repeated 50 times, viz., a 50×10-fold cross validation was implemented.

Experiments and Results

The programs were developed in-house using MATLAB (Version 2015a; Natick, MA), and run on an IBM laptop with 3 GHz Intel i3 dual-processor and 8 GB random access memory (RAM).

Key-slice selection

The curve of “inter-class variance” versus coronal slice index is shown in Fig. 6A. Ten coronal slices, from 60 to 150 with increasing steps of 10, were selected. The reason is that their corresponding ICVs were all higher than half of the maximum value. Figures 6B–6C shows respectively the axial and sagittal view of the selected key-slices, where the red lines represent the key-slices. The results are coherent with our formal work (Dong et al., 2015b).

Displacement field

Figure 7 shows the displacement field of an AD over the key-slices. Each column corresponds to a key slice. Here CI means coronal index. We show the moving image, the reference image, the overlap of both, the rigid registration, the overlap of rigid registration and the reference image, the non-rigid registration, the overlap of non-rigid registration and the reference image, and the displacement field for each key slice. From the row of displacement field, we can see how that a healthy brain is warped into an AD brain through the displacement fields (green arrows in the bottom row).

Region detection

We implemented the AD-related region detection procedure to different AD subjects as ‘Method of region detection’ described. Figure 8 shows the related regions, which are labelled by green points. Some areas are slightly outside of the brain, such as Fig. 8A, this is because the algorithm may also consider the distortion of background. This does not influence following labelling procedure.

Discussion

The results in Table 5 compare the proposed three classifiers (SVM, GEPSVM, and TSVM) with ten state-of-the-art methods. Plant’s results (Task 1 in Table 3 in Plant et al., 2010) presented the means together with 95% confidence intervals. Wang’s results (Table 7 in Wang et al., 2014) were obtained through a single K-fold cross validation analysis. Savio’s (Table 5 in Savio & Grana, 2013) and Dong’s results (Table 9 in Dong et al., 2015b) gave the means with SD.

Among the proposed methods, the proposed “DF + PCA + TSVM” yields the accuracy of 92.75 ± 1.77, sensitivity of 90.56 ± 1.15, specificity of 93.37 ± 2.05, and precision of 79.61 ± 2.21. Additional to it, the proposed “DF + PCA + GEPSVM” offers the accuracy of 91.52 ± 1.63, sensitivity of 88.93 ± 1.80, specificity of 92.27 ± 1.79, and precision of 76.66 ± 2.33. The proposed “DF + PCA + SVM” obtains the accuracy of 88.27 ± 1.89, sensitivity of 84.93 ± 1.21, specificity of 89.21 ± 1.63, and precision of 69.30 ± 1.91.

In terms of average accuracy, the proposed “DF + PCA + TSVM” result is as large as 92.75%, better than nine approaches of AD prediction, e.g., BRC + IG + SVM of 90.00% (Plant et al., 2010), BRC + IG + Bayes of 92.00% (Plant et al., 2010), BRC + IG + VFI of 78.00% (Plant et al., 2010), MGM + PEC + SVM of 92.07% (Savio & Grana, 2013), GEODAN + BD + SVM of 92.09% (Savio & Grana, 2013), US + SVD-PCA + SVM-DT of 90% (Wang et al., 2014), EB + WTT + SVM of 91.47% (Dong et al., 2015b), EB + WTT + RBF-KSVM of 86.71% (Dong et al., 2015b), and EB + WTT + POL-KSVM of 92.36% (Dong et al., 2015b). Nevertheless, the average accuracy of the method “DF + PCA + TSVM” is only less than one approach of “TJM + WTT + SVM” of 92.83% (Savio & Grana, 2013).

Why TSVM is better? There are two reasons. First, the non-parallel support vector machines provide more flexible and complicated hyperplanes than standard support vector machine. Second, the twin support vector machines formulate each of the two QP problems as a standard support vector machine, which makes it superior to generalized eigenvalue proximal support vector machine.

There were many other methods (Arbizu et al., 2013; Chaves et al., 2013; Cohen & Klunk, 2014; Dukart et al., 2013; Gray et al., 2012) proposed for detecting AD from NC, however, they dealt with images produced by other types of modalities: PET, SPECT, DTI, etc. Hence, it is inappropriate to compare the proposed methods with them. We will test our methods on SPECT and PET images in the future.

Table 6 shows that the displacement field finds the discriminant associated with the following regions reported in latest references: Angular Gyrus (Carbonell et al., 2014), Anterior Cingulate (Schultz et al., 2014), Cingulate Gyrus (Yu et al., 2014), Culmen (Almeida et al., 2011), Cuneus (Yokoyama et al., 2015), Fusiform Gyrus (Rieck et al., 2015), Inferior Frontal Gyrus (Eliasova et al., 2014), Inferior Occipital Gyrus (Liu et al., 2012), Inferior Parietal Lobule (Tramutola et al., 2015), Inferior Semi-Lunar Lobule (Almeida et al., 2011), Inferior Temporal Gyrus (Narayan et al., 2015), Insula (Fletcher et al., 2015), Lateral Ventricle (Gonzalez-Marrero et al., 2015), Lingual Gyrus (Cooley et al., 2015), Medial Frontal Gyrus (Jensen et al., 2015), Middle Frontal Gyrus (Yoo et al., 2015), Middle Occipital Gyrus (Cai et al., 2015), Middle Temporal Gyrus (Han et al., 2015), Paracentral Lobule (Kang et al., 2013), Parahippocampal Gyrus (Shimoda et al., 2015), Postcentral Gyrus (Willette et al., 2015), Posterior Cingulate (Pekarya, Sattin & Lloyd, 2015), Precentral Gyrus (Wang et al., 2015c), Precuneus (Villeneuve et al., 2015), Sub-Gyral (Ouyang et al., 2014), Superior Parietal Lobule (Yamashita et al., 2014), Superior Temporal Gyrus (Ramos et al., 2015), Supramarginal Gyrus (Redolfi et al., 2015), Uncus (Bangen et al., 2014).

Notwithstanding, some regions reported to be associated with AD are not interpreted by displacement field. Those areas include caudate nucleus (Montagne et al., 2015), claustrum (Pirone et al., 2015), lentiform nucleus (Dong et al., 2015b), subcallosal gyrus (Dong et al., 2015b), and subthalamic nucleus (De Reuck et al., 2014). Why are those areas not detected by our DF method? The reasons are tri-fold: (1) We only preserved the displacement field with magnitude longer than 5. Reducing the value of 5 may include more potential regions and noises. A feasible solution is to reduce from 5 to 3 and meanwhile develop robust anti-noise method. (2) Some literature used advanced imaging modalities, such as MRSI and fMRI for metabolism detection and function analysis. Therefore, we may also include these imaging-modality techniques. (3) The key-slice selection procedure may miss important regions. Hence, we will try to reduce the slice separation, although it will increase the computation burden.

The advantages of displacement-field are two-fold. On one hand, it reaches excellent classification performance, which was comparable to latest approaches (Dong et al., 2015b; Plant et al., 2010; Savio & Grana, 2013; Wang et al., 2014). On the other hand, it can directly locate AD-related regions. The disadvantages of displacement field are: (i) The accuracy of displacement field estimation relies on the accuracy of rigid registration, which is a necessary pre-processing (see Fig. 3). (ii) The estimated DF may exist in the background (see Fig. 8A), which needs to be removed by masking technique.

Conclusions

The contributions of this study consist of four points: (i) We proposed to use displacement field in the application of Alzheimer’s disease detection, and we proved its effectiveness. (ii) SVM and its two variants (GEPSVM and TSVM) were tested, and we proved TSVM performed better than SVM and GEPSVM. (iii) The proposed system achieved comparable sensitivity, specificity, and accuracy with ten state-of-the-art algorithms. (iv) The proposed computer-aided diagnosis (CAD) system can locate AD-related regions in the brain, which complies with 28 recent peer-reviewed articles.

Future research direction shall center in following aspects. First, the displacement field will be generalized to 3D, so key-slice selection can be discarded. Second, we shall try to embed kernels to classifiers. Third, spectroscopy method, like extracting the quantity of Docosahexaenoic Acid (DHA) (Tan et al.), may help to increase the classification accuracy. Forth, the proposed approach might be helpful in identifying brain mechanisms underlying endogenously defensive mechanism to neuroinjury and neurodegeneration (Yuan & Hou, 2015). Fifth, the least square technique (Hamedi, Salleh & Noor, 2015) may be applied to SVM and its variants, in order to further reduce training time. Sixth, the swarm-intelligence algorithm will be introduced to help to enhance the algorithm performance (Wang et al., 2015b). Finally, the stress field calculates the stress at every point within a particular region (Petit et al., 2015). It will be tested in future research.

Supplemental Information