Detection procedures for Patients of Alzheimer’s Disease using Waveform Features of Pupil Light Reflex in response to Chromatic Stimuli

INTRODUCTION: Various studies conducted to predict Alzheimer’s disease (AD) indicate that some pupil light reflex (PLR) features may contain symptoms of the disease. An effective procedure that can predict the disease using PLRs is needed. OBJECTIVES: Two analytic approaches were examined in order to estimate the possibility of identifying Alzheimer’s patients using features of PLR waveforms from chromatic stimuli. In particular, an index of the probability of being an AD patient is introduced, and the features which contributed to PLRs the most were extracted. METHOD: PLRs for three colours of light pulses (red: 635nm, blue: 470nm, white: CIE x=0.28, y=0.31) at two levels of intensity (10 and 100 cd/m2) were observed at 60Hz for 10s. Pulses consisted of pre-stimulus (2s), light pulse (1s) and restoration phases (7s). 15 features were extracted from each PLR waveform, such as pupil constriction velocity, pupil response delay, etc. Seven AD patients (age:42-84, mean=68.1) and 12 similar-aged control subjects (age:62-89, mean=72.1). RESULTS: The first approach was based on factor scores of features of PLRs. Two factor scores were extracted from the 15 features across all measurement conditions, and logistic functions were introduced in order to calculate the probability of identifying AD patients. Function parameters were estimated using a Bayesian technique, such as the Markov chain Monte Carlo method (MCMC). In consideration of the number of participants and biased data distributions, the second approach was based on the sparse modelling technique. Least absolute shrinkage and selection operator (LASSO) was applied to sets of PLR features from each light stimulus, together with the ages of subjects, and optimised result sets were obtained. Prediction performance was higher than with the previous procedure. CONCLUSION: The use of PLRs features from chromatic stimuli for identifying AD was developed and evaluated.


Introduction
The pupil light reflex (PLR), which produces changes in pupil diameter in response to a light pulse of detect symptoms of Age-Related Macular Degeneration (AMD) [5]. Some critical studies have suggested that common sources may be the origin of both AD and AMD diseases [6][7][8][9]. Also, because most AD patients are elderly, the influence of aging on PLRs should be evaluated carefully.
The authors have been studying a diagnostic procedure for detecting AD symptoms using PLRs of various types of light pulses and observing the conditions the light pulses produce [10]. Though these results show the possibility of aiding the diagnosis of the disease, a more flexible procedure is required because the prediction performance is insufficient to detect patients with the disease [11]. Though the authors have been examining a procedure to detect AD patients using PLR waveforms, new techniques should be considered to resolve some issues with the existing procedure [11]. In particular, the number of AD patients used in the experimental survey was restricted, and the assessment of the prediction of accuracy was insufficient. Since most participants for the experimental survey are elderly, they may already possess some artefact factors that would affect PLR observations, such as cataracts, or other age-related issues. Of course, a larger number of patients may contribute to accurate estimation of the parameters of the detection procedure, though this increases the number of parameters which need to be estimated, however. The number of samples has been discussed in order to be establish a causal relationship.
When the information from the data is insufficient, or only sparse data is obtained, some statistical procedures can be applied to the modelling of the causal relationship [12,13]. These techniques are introduced to construct a model. Therefore, significant features of PLRs which can be used in the diagnostic procedure should be extracted as needed, even if the number of samples is insufficient. Other factors should be estimated in order to reduce the influence of the problem of insufficient data mentioned above. If the distribution of features can be estimated, the possibility of diagnosing AD patients may be predicted using a Bayesian process, such as the Markov chain Monte Carlo method (MCMC) [13]. During the analysis, restriction of the distribution of features often influences optimising procedures for the prediction models since the extracted features are distributed sparsely. In order to analyse sparse distributed data samples, some sparse estimation techniques have been developed, such as the least absolute shrinkage and selection operator (LASSO) technique [12,14,15], and examination of the effectiveness of this estimation procedure is required.
In this paper, features of PLRs are extracted and the differences between AD patients and control group subjects are discussed, based on the results of our previous study [11]. Some procedures for predicting the probability of diagnosing AD patients are compared [16]. Also, the additional sparse estimation technique is applied to this prediction in order to improve its accuracy. The following topics are addressed: 1. Features of PLRs are analysed and their factors are extracted. The differences in features and factor scores between AD patients and control group subjects are compared. Also, the influence of aging is evaluated.
2. Logistic regression is introduced to calculate the probability of diagnosing the disease in AD patients and control group subjects. The models of fitness of the groups are then compared.
3. The MCMC technique is introduced to estimate the parameters of the models, and the performance of the models is discussed.
4. The LASSO technique is introduced to improve prediction performance, and selected features are discussed.
The reminder of the paper is structured as follows. Section 2 reviews previous and related work. Section 3 introduces methodologies to develop prediction procedures, such as feature characteristics of PLRs. Prediction using logistic regression, parameter estimation, and prediction with LASSO technique are detailed in section 4. The performance and validity of the methods are discussed in the Discussion and Conclusion sections.

Diagnostic for Alzheimer's disease (AD)
In order to diagnose AD patients and persons presenting symptoms of dementia, such as those with mild cognitive impairment (MCI) [20], the MMSE (Mini-Mental State Examination) question inventory [21] is often used during clinical consultations. As MMSE is used by observers to mark conditions during subject observation, the validity of the results depends on the activity of the subject during the examination. Therefore, some revisions have been introduced to improve the accuracy of consultations [22]. Also, more objective indices of subject's responses have been required in order to help medical doctors during clinical consultations.
Eye pupils and PLRs have been referred as one of the bio-markers for detecting AD patients [23]. In particular, PLR responses have been analysed in order to detect symptoms of AD, since ipRGCs and ipRGC-based PLRs have been discovered [24,25]. However, the diagnostic procedure using PLRs for AD patients continues to examine patients using 2 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6 Table 1. Summary of the related works

Predictions of AD patients
Computational diagnostic procedures for detecting AD patients or patients with dementia have been studied using human behavioural features. Unfortunately, the typical features which present AD symptoms are not currently being examined. Various information from human actions can be applied to the observation and prediction of cognitive performance using machine learning techniques. A summary of related works is listed in Table 1.
Since speech is observed by medical doctors during consultations, acoustic and lexical features are employed as a diagnostic procedure for a small scale set of data (9 AD patients and 9 control participants) [17]. The sound of a patient's speech can be recorded during medical consultation, and several projects have begun to summarise large sets of speech corpus data. Research groups can conduct computational prediction and classification of the level of diseases in participants using a corpus consisting of over 100 AD patients [18,26]. Performance depends on the data sets and the prediction techniques, and it does not seem easy to compare the accuracy of the two proposed procedures with previous studies.
Another approach is to observe human gait behaviour [19]. Since we need to pay a lot of attention to our bodies while we are walking, some interactions may affect this behaviour, causing our attention to be distracted. Recently, image processing techniques which can easily extract gait features are being applied to diagnosing elderly people. The possibility of dementia in 103 elderly participants was predicted using a machine learning technique [19].
As data deviation grows when the number of participants increases, AD patient prediction performance or prediction of the degree of development of the disease, including prediction using features of PLRs, does not necessarily improve through the use of larger-sized sets of data. In particular, feature extraction or feature assessment for making prediction is still preferred, before adaptations that would increase the size of data set are made.
Also, features of PLRs have been applied to the prediction of AD patients using conventional machine learning techniques [11]. AD patient prediction performance should also be improved using other metrics from the previous studies, though Table 1 shows an acceptable level of performance for predicting control subjects. As AD patients possess various stages of the disease, the performance of binary classification using categories such as "patient" vs. "healthy" is sometimes not useful for the diagnosis, and an intermediate assessment such as the level of MCI may be required. For this approach, the probability of diagnosing AD patients would help participants to understand the progress of the disease.
Therefore, more appropriate feature selection and a more appropriate prediction procedure should be considered, using the data sets for PLRs with various chromatic stimuli [11]. The following sections introduce new techniques in order to improve performance.

Method
An experimental survey was conducted with the participants. The detailed procedure that was used has been described in a previous study [11]. In this paper, the data set was used to develop prediction procedures.

Participants
A conventional PLR experiment was performed on 19 participants (42∼89 years old, mean age:70.6), 12 of which were healthy individuals with normal vision (Control group: 62∼89 years old, mean age:72.1) and 7 of which were patients with Alzheimer's Disease (AD Patients: 42∼84 years old, mean age:68.1) who had already been diagnosed by medical doctors.
Informed consent was obtained from all participants prior to the experiment.

PLR measurement
The stimuli consisted of three chromatic lights, red (635nm), blue (470nm) and white (CIE x:0.28, y:0.31), 3 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6 . These stimuli were labelled as r10, r100, b10, b100, w10 and w100. The duration of observations was 10 seconds, with the first 2 seconds being a pre-stimulus phase used as a rest period, followed by a 1 second light pulse and a 7 seconds restoration phase. Pupil diameters were measured in mm at 60Hz using a system developed by some of the authors [27]. PLRs for each stimulus were observed in single trials using a repeated-measure design.
Examples of measurements for a healthy individual and for an AD patient are shown in Figure 1. In these figures, PLRs are illustrated in response to 6 stimuli, namely the 3 colours and two levels of brightness. All subjects in the experiment participated 6 times, and a short break was taken between each set of observations.

Feature definitions
A typical PLR waveform shape is illustrated in Figure 2. In the figure, the light pulse overlaps for a period of 2∼3 seconds. As the figure illustrates, there are pupillary response delays due to the shrinking of pupil and its restoration to normal size.
Some features are extracted to specify the PLR response, and these variable features are summarised in Table 2. They are pupil size, velocity of pupillary change, duration of change, and integration of the waveform. The features correspond with the information in the small circles in Figure 2. In total, 15 features were extracted from the recorded pupillary changes using a MATLAB program. These features were calculated for each PLR response.

Characteristics of features in response to stimuli
Comparison of features. The features extracted for each stimulus were compared between the two groups. The results are summarised in Table 6. When there is a significant difference between pairs of values, the values are displayed in bold face. As the table shows, there are many significant pairs for the b100 and r10 conditions, but few pairs for the white stimulus. In regards to the significant differences for the b100 condition, such as pupil size, velocity and acceleration of pupillary change, the pupil size for AD patients is relatively small, and responds slowly. All participants are elderly, and in addition to being AD patients, their ages may affect pupil responses. Therefore, the effect of the two factors (participant group and age level) on pupillary changes is examined using two-way ANOVA. These means change along with age levels. Most variables selected are related to velocity and time delay. Since there are few significant interactions between the two factors, they may be independent of each other in regards to PLR features.
Factor analysis. There are significant differences in some of the PLR features between the AD and control groups. Though their variables exhibit the qualitative tendencies of a pupillary change, the sources could not be determined, as the physical variables and measurement units are completely different; some are expressed using sizes and others are expressed as velocities or accelerations. 4 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6  Factor analysis was used to extract latent sources of the variables, which were measured repeatedly in order to reduce the number of dimensions of the features of PLRs. As the number of features is too great for all to be evaluated, some features are correlated with others. In general, the factors extracted represent the variations. Since overall integration (int) is a summation of two parts, such as int_con and int_rest, have been omitted during the analysis, thus 14 variables were measured. A two-factor structure was estimated using a principal component solution and a screw plot. A factor loading  Table 4. The fundamental variables of each participant, such as pupil size, commonly contribute to both factors. The second factor contains variables which are concerned with features of the progress of restoration of the pupil after a pulse of light, and the first factor contains the remaining variables of the features of PLR. As mentioned above, both factors contain two variables, and there is a significant correlation between these two factors (r = 0.38). Since even the contribution ratio of each factor when the other factors are eliminated is over 60%, the two factors can account for the deviation of features of PLRs. 5 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6 The factor scores (f actor 1, f actor 2) were calculated using the factor loading matrix, and the means for each stimulus are summarised in Figure 3, according to group. The horizontal axis represents the first factor, and the vertical axis represents the second factor. The error bars show standard errors. The two groups are indicated using suffixes, such as "p" for patients, and "c" for the control group.
When means for red (r10 and r100) or white (w10 and w100) light stimulus between the two groups are compared, they are located proximate to each other, but the means for blue (b10 and b100) light stimulus of the two groups are relatively far apart from each other. Pupil reactions in response to light colours are represented in Figure 3. In the results of t-tests applied to pairs of means for the two subject groups, there are significant differences in the first factor scores for b100 (t(15) = 2.64, p < 0.05), in the second factor scores for b100(t(15) = 2.88, p < 0.05) and for w100(t(17) = 2.15, p < 0.05). Also, differences with little significance were observed for the second factor scores for b10(t(17) = 2.07, p < 0.10) and for r10(t(17) = 1.77, p < 0.10).
The results show that between the two subject groups there are significant differences in both factor scores for the b100 condition. As Table 6 shows many significant differences, the factor scores also represent these differences. The second factor seems to reflect the difference in features of PLRs between the two groups.
Although the influence of age level on factor scores was analysed, it did not affect either factor.

Introducing logistic regression
Function optimisation. As mentioned in the Introduction, this paper introduces the probability of diagnosing AD patients in order to estimate the progress of the disease using logistic regression analysis to predict binary response variables (p) and features of PLRs. Here, p = 1 for the control subject and p = 0 for the AD patient, so that p is given using the following equation with logit function.
Suffix i represents the subject, and suffix j represents the light stimulus condition.
Logistic regression analysis was applied to the above factor scores for several conditions, such as the two factor scores for a specific light stimulus condition (light intensity: 10 or 100 cd/m 2 ), the four factor scores for one colour condition (blue, red or white, including light intensities 10 and 100 cd/m 2 ) and further combinations such as two colours or all conditions (6 conditions and 2 factor scores). Every model was evaluated for fitness of model using AIC (Akaike  Information Criteria), and for prediction accuracy using R 2 . During validation of parameters, the likelihood or fitness of the prediction is evaluated using AIC or other statistical indices in addition to the level of accuracy. Accuracy was measured using an appropriate threshold, and performance was summarised using two dimensional metrics such as true and false positives. The relationships are then illustrated as ROC (Receiver Operating Characteristics) curves. Figure 4 shows ROC curves for every stimulus condition. The surface area of the curve is also a measure of the AUC (the area under the ROC curve). Their indices are summarised in Table  5.
The results of analysis suggest that discriminant performance is higher for blue light stimuli, in particular for the b100 condition. The ROC curves show step-wise changes, since the number of participants influenced the results. However, the performance of AUC for b100 produced the highest reaction of any single stimulus condition. 6 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6 Model parameter estimation. Table 6 shows the possibility of discriminating between AD patients and control group subjects using a logistic regression function. However, the parameters of these functions can not be estimated sufficiently because the amount of data is too limited. In order to obtain accurate parameters for the estimation function, a Bayesian technique such as the Markov chain Monte Carlo (MCMC) method was introduced. This technique generates distributions of parameters, and thus the most appropriate parameters can be extracted. As a result, a more accurate estimation of the parameters can be conducted. In regards to the data generation procedure using the MCMC technique, the burn-in period was 2000 and the number of samples was 10000 [28]. The parameters were estimated from the distributions generated in the 8000 samples which remained after convergences was confirmed. Using this procedure for the b100 condition, the parameters are estimated as follows.
The magnitude of the coefficient for the second factor (f actor 2) is 15 times that of the coefficient for the first factor (f actor 1).
As an additional example, the parameters for blue light stimuli were estimated as follows: The magnitude of the parameters suggests that the coefficient for the second factor at a high level of brightness (b100) is relatively larger than for others (b10 and the first factor at a low level of brightness). As the patterns of coefficients depend on the light stimulus, the wavelengths of the stimuli may affect these reactions.
In regards to the discussion in the previous section, the probability of correct discrimination between control group subjects (1) or AD patients (0) may be illustrated using two dimensional information (f actor 1 and f actor 2), as shown in Figure 5. Figure 5 shows that the probability distribution against patients with AD depends mainly on the scores of the second factor. Also, the score of the first factor helps to more finely adjust the probability during the period where the curve is steep.
In regards to the experimental procedure, the features of PLR can be best measured during a one second high brightness level pulse of blue light (b100). Also, factor scores were calculated using the factor loading matrix shown in Table 4. Finally, the probability of diagnosing AD patients can be predicted using the above function.
However, the possibility of developing a more flexible procedure for use in future experiments involving additional new participants will be a subject of our further study. 7 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6  b10&b100  r10&r100  w10&w100  ps_base  -------/

LASSO technique application
Optimisation using the LASSO technique. In order to improve prediction performance while the sparse distribution of pupillary features for PLRs is considered, an optimisation procedure such as LASSO (Least Absolute Shrinkage and Selection Operator) [12,13,29] is introduced for discriminating AD patients and controlled subjects. In the previous section, factor analysis was introduced to reduce the dimensions of the features of PLRs. Two factor scores cover 63∼84% of the overall features, as shown in Table 4. In this section, Error rate (%) 10-fold CV Figure 8. Comparison of error rates for detecting AD patients (Logit: previous section, RF: performance using Random Forest [11]) the features are selected and their performance is examined. Discrimination is made using a logistic regression function with the features of PLRs, as was done in the previous section. Since a logistic function is introduced as a link function, the probability of AD patients can be calculated. Here, performance is evaluated as a binary classification in order to optimise feature selection from all of the features as was done in the previous section. The LASSO technique calculates that some of the coefficients are driven to zero, thus making feature selection easy using the algorithm. The optimisation of selected variables and their weights for the prediction function 8 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6 using logistic regression is calculated using the LASSO procedure [14,15] with a package 'glmnet' [30]. An example of optimisation progress for a condition with a stimulus 10 cd/m 2 blue light (b10) is illustrated in Figure 6. The horizontal axis indicates parameter λ, and the vertical axis indicates weights of variables (β). As the figure shows, weights for all variables converge to zero along with the parameter λ. This indicates a variable selection procedure. When a certain value of λ is optimised, the appropriate weights for the variables are extracted from the features of PLRs.
Results of the optimisation with LASSO. Optimisation was conducted using AIC values which can be calculated using an index "sAIC" [31]. The optimised functions for each chromatic stimulus and for combinations of conditions, such as a blue light condition consisting two light intensities of 10cd/m 2 and 100cd/m 2 , were calculated. The optimised weights (βs) at the minimised AIC are summarised in Table 6. Since the minimum parameter λ is produced by calculating the 10-fold cross validation (CV) value, the optimised weights are influenced by variances during the procedure. Another set of AICs and CV values at the minimum λ is attached to the table.
The selected variables and their weights are summarised in Figure 7. The bars indicate the weights of the 9 conditions, while the variables which have have been omitted have no weights. In comparing the number of variables, more variables were selected for conditions of 100cd/m 2 than for condition of 10cd/m 2 . Some of the selected variables concern factors for time and speed of pupillary change. Also, the ages of participants for blue and red light conditions were selected when the variable "age" is added to the selectable variables.
In evaluating AIC as a model of fitness for single light conditions (b10, b100, r10, r100, w10, and w100), the minimum AIC occurs under conditions "w100" and "b10", as a result of optimisation using AIC. The area under the curve (AUC) is 1.00 for both conditions when a different method of classification performance is used. The CV rates are zero for all conditions using this procedure. When the optimum solution is based on the minimum of λ, the condition for "b10" is the lowest of all AICs for all of the single light conditions. The CV rates for "w10" and "w100" are the lowest for any of the conditions.
The evaluation was extended to converged conditions that converged (b10&b100, r10&r100, and w10&w100). The CV rates mostly decrease though the AICs increase with the number of variables, however. The 10-fold CV rates are summarised in Figure 8. The error rates using logistic function with factor scores from the previous section and using the random forest technique [11] are also displayed at the bottom of the bar graphs. As the figure shows, the error rate for the converged condition "b10&b100" is the lowest, and this rate is significantly lower than the prediction performance using logistic regression with factor scores and Random Forest technique.

Discussion
When logistic regression analysis was introduced, the probability of classifying AD patients was calculated using features extracted from PLRs. Prediction accuracy improved to nearly perfect using generated factor scores or some specific features of PLRs. The results confirm that some variables concerning factors for time and speed of pupillary changes contribute to the prediction of the probability of diagnosing patients with AD. In comparing the prediction performance with the performance in previous studies, the accuracy reaches to the levels of comparable without considering the sample size. However, as performance is based on the surveyed data set in this paper, a more generalised procedure may be needed using larger clinical samples. Before gathering increased numbers of samples, experimental conditions for stimulus using different colours of light and at varying intensities, and the validity of various features should be examined. These points are discussed as follows.
As shown at the end of the previous section, the prediction performance of the proposed procedures in our study should be compared and discussed. In some cases, performance using the LASSO technique is higher than for other techniques. When LASSO is applied to the converged condition "b10&b100", the error rate for 10-fold CV values can be reduced to 0.01. Also, some of the conditions can produce error rates lower than 0.10, as shown in Table 6. The results suggest that an optimisation technique can reduce the error rate of prediction performance for AD patients using PLR features. Also, optimisation 9 EAI Endorsed Transactions on Pervasive Health and Technology 09 2020 -12 2020 | Volume 6 | Issue 24 | e6 of prediction using the LASSO technique differs from other techniques. Conventionally, the contribution of the colour of the light, in particular blue light, is key to accurate prediction. Performance using PLR features of blue lights is the highest, while performance using features of white lights with conventional techniques is not high at all. As white light covers a wide range of light wavelengths, which stimulate ganglion cells on the retina, it is thought that the specific responses may not be detected. When the LASSO technique is introduced, error rates can also be suppressed, as is shown in Figure  8. Therefore, the technique can select the significant variables which are capable of predicting AD patients.
In order to confirm the weight selectivity which was hypothesised, weight patterns of three procedures for variables with the "b100" condition are compared in Table 7. Here, weight patterns using LASSO, importance values using Random Forest, and factor loadings using Bayes technique with factor scores are compared. The last factor loadings depend on whole responses, so they are not specific values with the "b100" condition.
As Table 7 shows, selected variables using LASSO have certain weights using the two other techniques, but some of the non-selected variables also have large weights. Regarding prediction performance using error rates, the LASSO technique may be able to select significant variables, as its performance is the highest. It may be a benefit of the technique, as the effectiveness of the LASSO technique has been confirmed. However, more detailed analysis of means to improve performance is needed in order to establish a more robust procedure for predicting AD patients will be a subject of our further study.

Conclusion
This paper presents a procedure for predicting the probability of diagnosing AD patients using features of PLR, which respond to the activities of ipRGCs.
Three colour lights at two levels of brightness were illuminated for 1 second, and pupil light responses were observed. 15 features were extracted from each PLR, and two factor scores were calculated using a factor loading matrix. In order to predict AD patients using these feature variables, two techniques were introduced: logistic regression and LASSO technique for optimisation. The following results were obtained.
1. There are significant differences in some features between AD patients and control group subjects, in particular for the b100 condition. Also, for a few features for white light there are significant differences between age levels.
2. Logistic regression analysis was introduced to discriminate AD patients from the control group using two factor scores in response to chromatic stimuli. Performance was evaluated using the indices of the fitness of equations. In the results, the performance for b100 was the highest.
3. The MCMC technique was introduced to estimate the parameters of the regression functions. The model provides a distribution of probability for AD patients and the control group. 4. In order to predict AD patients using the extracted features of PLRs, LASSO technique was introduced to logistic regression analysis. As a result, the error rate for prediction can be reduced to 0.01 using a converged condition of "b10&b100" where significant variables of features are selected.
5. In comparing weight patterns for feature variables of PLRs, the features which contributed to prediction were discussed.
The validity of the probability estimations should be confirmed using the PLR data of existing and additional new patients. This will be a subject of our further study.