Pupillary Complexity for the Screening of Glaucoma

Glaucoma is a leading cause of irreversible visual field loss. The early detection and diagnosis of the disease are therefore necessary to prevent blindness. Pupillary light responses are an interesting new technique for the detection of glaucoma. However, the analysis of pupillary signals has been associated with manual supervision or involved high computational costs. The present paper is to propose an analysis framework to automatically investigate changes in the complexity of pupillary signals under ambient light conditions for the screening of glaucoma. In this work, pupillary data of 13 glaucoma patients, 13 age-matched controls, and 11 young controls were recorded at the light intensity of 100 cd/m2 using a commercial eye tracker. The pupillary complexity of the participants was analysed using Higuchi’s fractal dimension, permutation entropy, and conditional entropy. We found that there was a statistically significant difference in the pupillary complexity between glaucoma patients and control groups (P < 0.0001). Specifically, the difference was more pronounced when using the fractal dimension measure. These results confirm the potential of using pupillary complexity for the screening of glaucoma using commercial devices.


I. INTRODUCTION
Glaucoma is characterised by an optic neuropathy associated with loss of the visual field. The disease is the second largest leading cause of irreversible vision loss worldwide. Glaucoma affected over 60 million people in 2010 [1]. This number will increase to more than 111 million in 2040 [2]. Early diagnosis and treatment of glaucoma play a vital role in the prevention of vision loss for millions of people.
The common glaucoma tests are tonometry, gonioscopy, ophthalmoscopy, perimetry, optical coherence tomography, and pachymetry. Among these, standard automated perimetry is the most popular method for detecting glaucoma [3]. Optical coherence tomography can identify the disease accurately [4]. These tests are generally available only at ophthalmology clinics.
The discovery of melanopsin expressing intrinsically photosensitive retinal ganglion cells has helped to explain the regulation of pupil size [5]. This helps to widen the The associate editor coordinating the review of this manuscript and approving it for publication was Yongming Li . understanding of the retinal and optic nerve disorders, including glaucoma, through pupillometry recordings [6], [7].
Abnormalities in pupillary light reflex (PLR) have been reported in patients with glaucoma. The melanopsinmediated post-illumination pupil response (PIPR) was used to detect melanopsin dysfunction in patients with early glaucoma [8]. The PLR was measured using 1-second pulses of blue lights (λ max = 464nm) and red lights (λ max = 658nm). The post-illumination pupil response to light stimuli has been found to alter in moderate and advanced stages of glaucoma using 10-second light stimuli of 488 nm (blue) and 610 nm (red) [9]. There was a significant decrease in PIPR of glaucomatous patients compared to age-matched controls using two wavelengths (470 nm and 623 nm) [10].
Open-angle glaucoma is commonly bilateral and asymmetric. The difference in the responses between the two eyes could be used as a biomarker for glaucoma [11], [12]. Relative afferent pupillary defect (RAPD) was detected in glaucomatous patients with 66.7% sensitivity and 82.9% specificity using the swinging flashlight test [13]. In [3], 0.2 seconds of coloured light stimuli of different wavelengths (red, green, VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ yellow, and blue) were used to study the difference in RAPD between glaucoma patients and controls. The above studies show the potential use of the pupillary response for detecting glaucoma. However, these studies are highly dependent on protocols exploited for chromatic pupillography. The studies were involved with the manual selection of dark adaption time and pulse-width, and different wavelengths of light stimuli [14], [15]. In addition, the analysis needed manual supervision. Hence, there is a need for automatic recording and analysis of pupillary data for assisting in the diagnosis of glaucoma that can work under ambient light conditions. Non-luminance-mediated changes in pupil size have been used extensively as markers of arousal and cognitive effort [16]. The oscillations in the pupil diameters, referred to as hippus, were found to originate from the parasympathetic nervous system [17]. Therefore, changes in these oscillations may reflect the dysfunction of neural pathways in the visual system. The use of power spectral density and other spectral features were found to be unsuitable for the analysis of such recordings [18], [19]. It has been reported that sample entropy to estimate the complexity of pupillary data can differentiate between glaucoma and healthy eyes [20], [21]. It was observed that the pupillary complexity was higher in glaucoma eyes than in healthy eyes.
Signal complexity, a concept related to the irregularity of a signal, can be investigated using a number of non-linear algorithms. Complexity measures are helpful in detecting changes in the dynamics of biological systems that can be associated with physiological events [22]. For analysing biomedical signals, entropy-based algorithms have been used extensively. There are various types of entropy-based complexity parameters for real data: approximate entropy, sample entropy, and permutation entropy. Although approximate entropy and sample entropy measure irregularity, they give high values to random data and do not imply an intrinsic physiological complexity [23], [24]. In addition, approximate entropy and sample entropy are very sensitive to its input parameters: similarity criterion (r) and data length (N ) [25].
We propose the analysis of the complexity of pupillary signals in patients with glaucoma using Higuchi's fractal dimension (FD), permutation entropy, and conditional entropy in the present study. These complexity indices have a lower computational cost than sample entropy (∼ O(N ) vs.
∼ O(N 2 )) [26], [27]. Permutation and conditional entropy parameters are robust to dynamical and observational noise and less dependent on the data length [28], [29]. These entropy measures have been successfully used to analyse electroencephalogram (EEG) signals and electrocardiogram signals [26], [30]. In an EEG-based hypoglycemia detection study, Higuchi's FD was found to provide better classification results of hypoglycemic episodes compared to entropy-based algorithms [26]. We aim to develop a framework for the endto-end investigation of pupillary recordings under ambient light conditions, including pre-processing, signal conditioning, unsupervised segmentation, and feature extraction of the recordings. In this pioneering work, the FD, permutation entropy, and conditional entropy in pupillary recordings of glaucoma patients, age-matched healthy controls, and young controls were computed and compared to identify the differences between glaucoma and healthy eyes. The computational cost of these was recorded for determining the computerised and real-time implementation.
The rest of this paper is structured as follows. Section II demonstrates the study protocol, pupillary signal preprocessing, complexity analysis, and statistical analysis. We report the results of this study in Section III. This is followed by Section IV, which provides a discussion on the obtained results. Finally, we conclude the entire work in Section V.

A. STUDY PROTOCOL
This study was conducted in accordance with the ethics approval from the RMIT Human Research Ethics Committee. Informed consent was obtained from all participants prior to data recording. The protocol and data of the study have been reported earlier in [20] and [21].
Thirty-seven participants volunteered for this study: thirteen glaucoma patients, thirteen age-matched controls (aged 66.231±8.097 years) and eleven young controls (aged 30.455±2.903 years). Glaucoma patients had early-stage, bilateral, open-angle glaucoma who were recruited from Essendon Eye Clinic & Laser Centre (average age = 70.154±8.160). All patients were being treated for intraocular pressure (IOP), and their IOP was in the range of 17 -18 mm Hg. All participants were asked not to consume caffeine for at least 1 hour before data recording. The demographic details of the participants are shown in Table 1. Participants' pupil diameters were measured using a commercial eye tracker, GP3 (Gazepoint, Canada), with a sampling rate of 60 Hz. The device was attached to a computer screen and placed approximately 65 cm from the participant's eyes. Before each experiment, the GP3 device was calibrated using the provided company's software and protocol. During the experiment, scale factors and pupil diameters were automatically recorded by the eye tracker software. The scale factor reflected changes in the distance between the eyes and the device. A scale value of 1 represents the distance at the time of calibration. When the participant is closer to the eye tracker, the scale value is less than 1, and when the participant is further away, the scale value is greater than 1.
The experiments for pupillary data collection were conducted during regular office hours (9:00 AM to 5:00 PM), in a standard clinic room and with ambient light intensity near the eyes of participants being 100 cd/m 2 . Participants were requested to focus on a blue point target at the centre of the screen for 120 seconds. During the experiment, they were sitting comfortably on a fixed chair. The experimental setting is shown in Fig. 1.

B. PRE-PROCESSING
The pupillary signals were multiplied by scale factors to correct movement artefacts. The next step of pupillary data analysis was to remove outliers using Hampel filters (sliding window size = 120 and number of standard deviation = 3). After that, a 0.02 to 4 Hz Butterworth band-pass filter was used to extract pupillary information since the typical frequency range of pupillary activity is less than 4 Hz [31], [32]. The selection of the low-cut frequency of 0.02 Hz for the filter was to remove the baseline wander and DC in the recorded signals [17]. Fig. 2 shows raw and filtered pupillary signals from the left eye in the first 60 seconds of a representative participant.
In this study, the recording during the transition period was removed automatically. This period was determined using the time stamp from the eye tracker. We only analysed signals when the eyes were in the steady-state, in which the eyes adapted to the light intensity when looking at the screen. Thus, the initial 10 seconds of pupillary recordings were removed. The amount of 10 seconds is more than the time to get the eye to reach its steady-state [33]. We also removed the last 10 seconds of the recordings to reduce any interference at the end of the experiment.
The resulting signals were divided into non-overlapped 5-second segments. This helped to increase the number of samples for the analysis. The complexity parameters were then computed for each of these 5-second segments using our proprietary software written in MATLAB language (MATLAB 2017b, The MathWorks, Natick, MA, USA).

C. HIGUCHI'S FRACTAL DIMENSION MEASURES
Higuchi's FD is a measure of signal complexity [34]. It is based on the relationship between the length of a curve L(k) and the scale k. Given a segment x of pupillary signals, the difference between samples at distance k can be computed as follows: L(k) is then computed as: If L(k) is proportional to k −D for k = 1,. . . , k max , the curve is fractal with dimension D. As a result, D is estimated as the linear coefficient of the regression line between log(L(k)) and log(k).

D. ENTROPY MEASURES
The permutation entropy (PE) of order d ∈ N and delay τ ∈ N of a segment x t is calculated as: where t = 0 . . . N − 1, N ∈ N and is the relative frequency of ordinal patterns π in the time series and 0 ln 0 is defined by 0. The vectors v(x t , x t−τ , . . . , x t−dτ ) are delay vectors.
An alternate is conditional entropy (CE) which is akin to permutation entropy and has been shown to provide better performance in similar applications [35]. Thus, for this study, conditional entropy was investigated.
The conditional entropy of a segment x t is defined as: in which p τ π is defined in (4) and is the relative frequency of ordinal patterns π 1 , π 2 in the time series.

E. STATISTICAL ANALYSIS
In this study, we used t-tests to identify the significant differences in pupillary complexity between the groups with the hypothesis of normal distribution of samples. T-tests are suitable for two-group comparison and were conducted to investigate the difference between i) glaucoma patients (G) and age-matched controls (H), ii) glaucoma patients and young controls (Y), and iii) age-matched controls and young controls. Features having P-values < 0.05 are considered statistically significant. In addition, Cohen's d effect size was computed to discuss the differences between the groups. Logistic regression is a binary classifier and is routinely used in bio-signal analysis and, in this study, was used to classify glaucoma and healthy eyes. The extracted complexity parameters were the input for the regression model. Receiver operating characteristic (ROC) curves were plotted to describe the classification using the output of the logistic regression model with every possible cut-off value. The area under the curve (AUC) was computed.

III. RESULTS
The pupillary dataset in this study was composed of 260 samples of glaucoma eyes, 260 samples of the age-matched healthy eyes, and 220 samples of healthy eyes of young people. Each sample includes 5-second pre-processed segments of pupillary recordings from the left and right eyes. For the age-matched control group, 7 samples were removed from the dataset because the files were corrupted.
From each sample, three complexity features were computed: Higuchi's fractal dimension, permutation entropy, and conditional entropy. A range of values of the parameters for each feature was tried, and the ones that gave the best results were selected. For Higuchi's fractal dimension, k max = 12, for the entropy measures, the order d = 6, and the delay τ = 1. Finally, the dimensions of glaucoma eye, age-matched healthy eye, and young healthy eye subsets are 260 × 6, 253 × 6, and 220 × 6, respectively. Higuchi's fractal dimension of pupillary signals for glaucoma (G), age-matched (H), and young (Y) participants. The level of significance is represented by asterisks ( * * * * means P < 0.001 and * * * * * means P < 0.0001). Error bar: 95% CI. Fig. 3 shows that the fractal dimension of pupillary signals in glaucoma eyes is significantly higher than age-matched healthy eyes (right eyes, P < 0.0001). In contrast to the right eyes, the level of significant difference found from the left eyes was 0.05. Notably, the fractal dimension in glaucoma patients was significantly higher than in young people. There were also significant differences in the fractal dimension between age-matched controls and young participants (left eyes: P < 0.001 and right eyes: P < 0.0001).
Permutation entropy indices of glaucoma eyes and healthy eyes are shown in Fig. 4. There was a significant difference in the right eye complexity between glaucoma patients and agematched controls (P < 0.0001). However, it can be seen that there was no significant difference among the control groups: the young and the age-matched. Fig. 5 shows the results of conditional entropy of pupillary recordings. The conditional entropy indices calculated from the right eyes of the glaucoma patients were significantly higher than the control groups (P < 0.0001). Conditional entropy of the right eyes also shows that there is a slight difference in pupillary complexity between older controls and young controls (P = 0.06).
The relationship between left eyes and right eyes in the three groups is shown in Table 2. It shows that there are no statistically significant differences between the left eyes and the right eyes in older people (both glaucoma patients and agematched controls). Nevertheless, for young controls, the table reveals the complexity of pupillary signals calculated from the right eyes is significantly different from the left eyes. Among the three complexity measures, permutation entropy provided the highest significant level (P = 0.002).
The plots of the ROC curves obtained using the logistic classifiers for i) discriminating glaucoma eyes and age-matched healthy eyes (G-H) and ii) discriminating    and age-matched healthy eyes with an AUC of 0.65. When comparing glaucoma eyes and young eyes, the AUC increases to 0.75. From the G-H curve, the classification results of glaucoma eyes were 63% sensitivity and 60% specificity. From the G-Y curve, we achieved a sensitivity of 72% and a specificity of 63%.

IV. DISCUSSION
This study has investigated the pupillary complexity in three different groups: glaucoma, age-matched controls, and young controls with the aim to develop the outline for automatic analysis of the detection of glaucoma. The proposed framework for automatic feature extraction of pupillary signals was composed of signal conditioning (i.e. Hampel filters and bandpass filters) and complexity measurements. In this work, we used Higuchi's fractal dimension, permutation entropy, and conditional entropy to measure the complexity. We found that glaucoma patients had a higher complexity than other groups, irrespective of age. These findings are consistent with the results from sample entropy-based algorithms in previous studies [20], which, however, was not suitable for computerised analysis.
The computational time of complexity for a sample of pupillary signals using different algorithms is shown in Table 3. We used MATLAB's built-in time functions on the hardware with Intel Core i5-7200U CPU 2.71 GHz processor and 8 GB RAM for the estimation of time consumption. For comparison, we also calculated the computational time of sample entropy using the parameter settings described in [21]. While the table shows that Higuchi's FD required the lowest computational time (4.44±0.68 milliseconds), however, the difference may not be relevant for practical applications. What is important is that while earlier studies required manual segmentation, the proposed approach was performed unsupervised and hence suitable for computerised analysis. The selection of k max in Higuchi's fractal dimension depends on signal characteristics and applications [36]. Generally, it was based on the value that provided the best results for the discrimination between two or more features. For electroencephalogram-based applications, k max was set equal to 6 [37]- [39]. In the present study on pupillary signals, a range of k max was investigated to figure out the value that provided the best statistical results for the analysis. We found that k max = 12 was a reasonable value for investigating the complexity of pupillary signals using the fractal analysis. For the calculation of entropy in this study, the same judgement was used in the selection of parameters.
As pupillary signals are non-stationary [19], the selection of segment lengths may affect the results. In this study, 5-seconds segments of pupillary signals were used for the analysis. By removing the first 10-second and the last 10-second of each recording, we avoided any transitions in the signals. The duration of 10-second is enough for the eye to reach its steady-state [33]. The non-overlapped segments were extracted in this period of the steady-state response of the pupils.
The results of Cohen's d effect size calculation are shown in Fig. 7-9. It is worth noticing that all the effect size values for the difference between glaucoma eyes and age-matched healthy eyes were greater than 0.35. Especially, the effect size for the fractal dimension when comparing glaucoma eyes and healthy eyes of young participants was larger than 0.8. A benchmark of the effect size has three thresholds: small (d = 0.2), medium (d = 0.5), and large (d = 0.8) [40]. These give evidence that pupillary complexity could be used as a biomarker for glaucoma.   The correlation between Higuchi's fractal indices with the other entropy-based indices was highest in glaucoma patients. Among the complexity measures, the distribution of Higuchi's fractal indices was more symmetrical compared to those of PE and CE, which were left-skewed. The results from the effect size and correlation analysis show that, besides entropy-based algorithms, Higuchi's fractal dimension can be a good candidate for the analysis of pupillary complexity In this study, statistically significant levels of differences were found i) between the left eyes and the right eyes and ii) between older controls and young controls. Specifically, the level of significant difference when comparing glaucoma eyes and age-matched eyes was much higher with features  extracted from the right eyes. The results from Higuchi's fractal dimension showed significant differences in healthy eyes between older people and young people. In order to achieve a robust conclusion about these observations, a large number of participants is required. This is a pioneering study in which the data were collected using a commercial and portable eye-tracking device. The procedure for data collection does not require highly skilled personnel or even a special clinical setting. The proposed automatic analysis framework could be used in primary healthcare facilities for the screening of glaucoma.
There are some major limitations in the current study. Firstly, this is a pilot study in which the number of participants is not sufficient to investigate differences due to factors such as gender and ethnicity, and due to other eye diseases or other clinical factors. We performed binary comparisons: glaucoma subjects vs. age-matched controls, glaucoma  subjects vs. young controls, and age-matched controls vs. young controls. While the number of participants is low, the analysis indicates that there were statistically significant differences in pupillary complexity between glaucoma and controls. By extracting 20 samples from each eye of each participant, we achieved 260 samples of glaucoma eyes, 260 samples of age-matched healthy eyes, and 220 samples of young healthy eyes. These numbers satisfied the statistical power of 80%. However, we are aware that these samples are not independent, and this may affect the results. Therefore, a larger number of participants is required to validate the findings. Secondly, this work has not investigated the correlation between pupillary complexity and the severity of glaucoma. We only studied OAG patients with a treated IOP in the range of 17-18 mm Hg. Finally, a neurological interpretation of the results has not been discussed in this study. In real-life situations, changes in pupillary complexity may reflect other ophthalmological diseases. Further work will be required to investigate different eye diseases and also other confounding factors. VOLUME 9, 2021

V. CONCLUSION
Our pioneering study has investigated the difference in complexity of hippus between glaucoma and healthy eyes and suggested the framework for detecting glaucoma, which has the potential for automatic and unsupervised data recording and analysis. This study on a total of 37 participants shows the potential of pupillary complexity as a biomarker for glaucoma, which can be recorded using an inexpensive eyetracker, and the pupillary recordings can be automatically analysed. The fractal dimension of pupillary signals of glaucoma patients was significantly higher than the age-matched and young control groups. The same results were found using permutation entropy and conditional entropy. This shows that pupillary data has the potential for detecting glaucoma.
Future development of the current work will be related to investigating the effect of progression of the disease and the visual acuity, differences between gender and ethnic groups and identifying the effect due to other eye diseases and confounding factors. There is also the need for optimisation of the proposed framework on a broader pool of participants.

ACKNOWLEDGMENT
Our studies in this article were conducted with the assistance of the Essendon Eye Clinic & Laser Centre, Melbourne, Australia.