Visualizing the dynamic change of Ocular Response Analyzer waveform using Variational Autoencoder in association with the peripapillary retinal arteries angle

The aim of the current study is to identify possible new Ocular Response Analyzer (ORA) waveform parameters related to changes of retinal structure/deformation, as measured by the peripapillary retinal arteries angle (PRAA), using a generative deep learning method of variational autoencoder (VAE). Fifty-four eyes of 52 subjects were enrolled. The PRAA was calculated from fundus photographs and was used to train a VAE model. By analyzing the ORA waveform reconstructed (noise filtered) using VAE, a novel ORA waveform parameter (Monot1-2), was introduced, representing the change in monotonicity between the first and second applanation peak of the waveform. The variables mostly related to the PRAA were identified from a set of 41 variables including age, axial length (AL), keratometry, ORA corneal hysteresis, ORA corneal resistant factor, 35 well established ORA waveform parameters, and Monot1-2, using a model selection method based on the second-order bias-corrected Akaike information criterion. The optimal model for PRAA was the AL and six ORA waveform parameters, including Monot1-2. This optimal model was significantly better than the model without Monot1-2 (p = 0.0031, ANOVA). The current study suggested the value of a generative deep learning approach in discovering new useful parameters that may have clinical relevance.

different between the two eyes 7 . Furthermore, we previously reported that AL may increase in adulthood 11 , in line with a previous paper 12 . In line with this, we previously reported that the correlation between AL and the cpRNFL peak angle (r = −0.49) or the PRAA (r = −0.38) was moderate 7 .
Ocular biomechanical properties can be measured by using the Ocular Response Analyzer (ORA; Reichert Inc., Depew, NY, USA) and Corvis Scheimpflug Technology (CST; Oculus, Wetzlar, Germany). By analyzing CST biomechanical parameters, we previously reported that the ability to absorb the applied external energy (hysteresis) was significantly associated with myopic retinal stretch as estimated by the cpRNFL peak angle 13 . Another study reported that the maximum deformation amplitude as measured by using CST was associated with the size of β-zone parapapillary atrophy 14 . With the ORA, it is possible to evaluate corneal biomechanical properties in detail by analyzing the recorded waveform rather than the single parameter of Corneal hysteresis (CH). For instance, we have reported that parameters extracted from the ORA waveform were more strongly correlated with glaucomatous visual field progression compared to CH 15 . By using this approach, we recently reported that the corneal biomechanical properties described by parameters extracted from the ORA waveform were significantly related to myopic retinal deformation 16 . Currently, in the ORA 37 parameters are shown, which are derived from the two ORA response wave peaks, as implemented by the manufacturer. However, in addition to these parameters, it is likely that analysis of the ORA waveform by means of advanced algorithms (e.g., machine learning) may be useful to capture further aspects of ocular biomechanics.
In machine learning, including deep learning, discriminative and generative models are used 17 . Several studies suggested the usefulness of a discriminative deep learning approach in Ophthalmology, for example for diagnosing glaucoma from a fundus photography [18][19][20][21][22] and from optical coherence tomography (OCT) 23 . Variational Autoencoders (VAEs) are a type of deep learning approach that allows powerful generative models of data 24,25 . However, so far the usefulness of generative deep learning in Ophthalmology has not been assessed. VAEs consist of an encoder, a decoder, and a loss function. The input data is first processed using the encoder (a neural network), represented as a multidimensional probability density in a latent space, and then reconstructed by the decoder (a neural network). VAEs have demonstrated remarkable generative capacity and modeling flexibility, especially with imaging data 24 . Indeed VAEs have been used for various purposes, such as anomaly detection (for example, in Electrocardiograms 26 ), clustering, and in particular, noise filtering 27 . Another feature of the VAEs approach is that it enables visualization of the dynamic change of input data by gradually shifting the latent variables, which may be helpful to understand its characteristics 24 .
The aim of this study was to assess whether possible changes in ORA waveform associated with changes in PRAA could be detected by using new ORA waveform parameters extracted from the ORA waveform by means of a VAE approach.

Methods
Study population. The study protocol was approved by the institutional review boards of the University of Tokyo Hospital, the University of Hiroshima Hospital, and the Tsukazaki Hospital and adhered to the tenets of the Declaration of Helsinki. Informed consent was obtained from each subject.
The sample of eyes included in the current study was the same as in our previous study 16 . ORA measurements were conducted in forty-nine normal eyes from 47 subjects. Inclusion criteria were: 1) no pathological findings by slit-lamp microscopy, ophthalmoscopy, and/or OCT; (2) best-corrected visual acuity ≤0.1 LogMAR (logarithm of the minimal angle of resolution); and (3) intraocular pressure ≤21 mmHg, as measured by using Goldmann applanation tonometry. Exclusion criteria were: (1) known ocular diseases such as glaucoma, staphyloma, and optic disc anomalies; (2) systemic diseases such as hypertension and diabetes; (3) the presence of visual field defects; and/or (4) a history of refractive or intraocular surgery, including cataract surgery.
Measurements of AL and SERE. The AL was measured by using an optical biometer (OA-2000; Tomey, Nagoya, Japan). Three subsequent measurements were taken and the average value was used as AL measure. The SERE was measured by using the Topcon KR8800 autorefractometer/keratometer (Topcon, Tokyo, Japan). peripapillary retinal arteries angle (pRAA). The methodology to measure the PRAA has been reported in detail elsewhere 16 . Briefly, optic disc color fundus photographs were obtained by using either an OCT (OCT-2000 ® ; Topcon, Tokyo, Japan) or a fundus camera (TRC-50DX ® ; Topcon). ImageJ software (https://imagej.nih. gov/ij/; provided by the National Institutes of Health, NIH, Bethesda, MD) was used to draw a 3.4-mm-diameter peripapillary scan circle on the obtained fundus photographs, and the PRAA was calculated as the angle between the radia crossing the points at the intersection between the 3.4 mm-diameter peripapillary scan circle and the supratemporal/infratemporal major retinal arteries (as shown in Fig. 1). Magnification effects of the camera were corrected by using the Littmann's formula 28 . Analysis of the oRA waveform. The ORA waveform is recorded by using an electro optical system consisting of a collimated beam of infrared light that is reflected by the surface of the cornea onto a photodetector and represents the deformation of cornea (inward and outward movements) in response to a rapid air jet. The ORA waveform is composed of two peaks that correspond to the corneal inward and outward applanation events. The CH is defined as the difference in the two applanation pressure values and it represents the viscous damping of the corneal tissue 29 . The corneal resistant factor (CRF) is also calculated from the two applanation pressure values, but places greater emphasis on the first applanation pressure in order to give more information about the elastic properties of the cornea 29 . By using the ORA software (version 3.01), 37 waveform parameters can be obtained from the ORA waveform. The ORA waveform was recorded by using an ORA G3 model with the related PC software for waveform analysis. The ORA measurements were carried out three times with at least a 5-minute www.nature.com/scientificreports www.nature.com/scientificreports/ interval between each measurement, and the average of each obtained values were used in the analysis. All the measurement had a quality index of higher than 7.5 as recommended by the manufacturer.
Variational autoencoder. The structure of the VAE model used in this study is shown in Fig. 2. The encoder is a 1-layer neural network consisting of 400 units (one for each of the 400 ORA waveform observation points). This encoder is connected to 2 hidden layers consisting of 40 and 20 units, and is then represented by the mean and logarithmic variance-covariance matrix of a two-dimensional Gaussian probability density in the latent space. The decoder reconstructs the 400 units through a further 2 hidden layers and 1 output layer, which represents the reconstructed ORA waveform. This VAE model was optimized by maximizing the sum of the negative reconstruction loss obtained by using the current study's dataset, which is defined as the Kullback-Leibler divergence between the distributions of the differences between the input ORA waveform and reconstructed ORA waveform. Then, the reconstructed ORA waveform was analyzed in conjunction with the PRAA.
identification of new oRA waveform parameter(s). By visual analysis of the changes in VAE-reconstructed ORA waveforms as a function of the changes in PRAA, we identified the change in monotonicity (the retrogressive movement) between the first and second applanation peaks (Monot1-2) as a possible new ORA waveform parameter potentially sensitive to changes in PRAA. More specifically, Monot1-2 was estimated as the total length of retrogressive movement between applanation peak 1 and applanation peak 2 (Fig. 3).
Statistical analysis. The relationships between the PRAA and the 41 variables of age, AL, SERE, CH, CRF, 35 out of 37 ORA waveform parameters, and the newly proposed parameter Monot1-2 were evaluated by using a  www.nature.com/scientificreports www.nature.com/scientificreports/ two-step feature selection approach in view of the high number of variables (41), following our previous study 16 . Two out of 37 ORA waveform parameters (i.e., h11 and h21) were not used as they are proportional to h1 and h2, respectively. First, 20 candidate variables were selected using the least absolute shrinkage and selection operator (Lasso) regression. Then, model selection was used to identify the optimal model for the PRAA by using the second-order bias-corrected Akaike information criterion (AICc) index from 2 20 different patterns generated by using the 20 candidate variables. The Akaike information criterion (AIC) is a statistical measurement used in model selection 30 and the AICc is a corrected version of the AIC, which is able to provide an accurate estimation even when the sample size is small 31,32 . A decrease in AICc indicates an improvement in the model 33 and suggests that the variables selected through the model are significant 34 . It is worth noting that multivariate modeling such as the one here used can be useful for detecting patterns and characterizing data because it provides more control over potential confounders compared to univariate analysis 35 . The log-likelihood values of paired models (e.g., with vs without Monot1-2, multivariate vs univariate) were compared by using the analysis of variance (ANOVA) test.

Results
The demographic characteristics of participants and of the ocular properties of the studied eyes are shown in Table 1.
The change of the reconstructed ORA waveform using VAE associated with the decrease of PRAA is shown in Supplementary Video 1. Table 2 shows the summary descriptive statistics of the ORA parameters. Table 3 shows the results of univariate and multivariate linear regression between PRAA and the values of age, AL, CH, CRF,

Discussion
In the current study, we assessed the dynamic change of the ORA waveform in relation to retinal deformation as estimated by PRAA in a sample of 49 eyes from 47 participants, using VAE approach, a deep learning generative model. This study highlighted that a novel parameter extracted from the ORA waveform (Monot1-2) may be used to generate multivariate models of the PRAA that are more accurate than models without Monot1-2.

Variables
Mean ± Standard Deviation (Range) www.nature.com/scientificreports www.nature.com/scientificreports/ The optimal model for the PRAA in the previous study was as follows: PRAA = 68.6 -3.0 × AL -3.1 × CRF + 9.5 × Aindex + 1.8 × w2 + 0.40 × slew1, when the new ORA waveform parameter of Monot1-2 was not included 16 . Such waveform parameters represent a quick response of the cornea to external forces and suggest a soft cornea 36,37 . In the current study, eyes with high myopia were added, and the optimal model for PRAA was PRAA = 238. In the current study, AL was negatively associated with PRAA, similarly as in our previous study 16 . Besides, the ORA waveform parameters selected in the current study indicated that decrease in the size of the area of applanation (peak1, aspect2, uslope11, and h1), increase in maximum single length of the outside line segments of peak1 (mslew1), and increased lability in restoration phase/peak2 (Bindex) were related with smaller PRAA 36 , which infer quick cornea response and soft cornea behavior in line with our previous study 16 . In addition to these parameters, the  www.nature.com/scientificreports www.nature.com/scientificreports/ current study showed that larger Monot1-2 was significantly related to smaller PRAA. This might be because eyes with smaller energy dissipation might have a greater amount of stored elastic energy during applanation 1 and 2, causing an increase in Monot1-2. Thus, results of the current study also seem to suggest eyes with greater myopic retinal deformation may demonstrate decreased energy dissipation.
The previous reports show that small energy dissipation of an eye represented by CH is related to rapid progression of glaucoma [38][39][40] . This may be because of similar biomechanical properties of the cornea and sclera 41 , as they are made up of the same types of collagen 42 . In addition, eyes experience great changes in the intraocular pressure even in daily life events, such as postural change 43 , eye lid blinking 44 , ocular pulsatility due to ocular hemodynamics 45 , Valsalva maneuver 46 , and eye movements 47 . Reduced energy dissipation in the cornea may cause higher vulnerability to these stresses 41 . Furthermore, we recently reported that angioid streaks were significantly associated with the corneal biomechanical properties as measured by using CST 48 . Considering these facts, the newly proposed Monot1-2 may also be useful in these diseases, and this should be further investigated in a future study.
A limitation of this study is the cross-sectional design. Further validation would be needed by using a longitudinal research approach, in particular in young population. We investigated the relationship between ORA waveform parameters and PRAA; however, subanalyses in association with the severity of myopia were not carried out due to the sample size in the current study. Further investigation is preferable with the larger population.
In conclusion, in this study a new ORA waveform parameter was proposed by analyzing, by using VAE, the dynamic change of the ORA waveform in relation to retinal deformation as estimated by PRAA. This was the first study to demonstrate the value of generative deep learning models such as the one generated by VAE, in discovering new useful parameters that may be helpful in the clinical setting.