Remodeling the light-adapted electroretinogram using a bayesian statistical approach

Abstract

Objective

To present a remodeling of the electroretinogram waveform using a covariance matrix to identify regions of interest and distinction between a control and attention deficit/hyperactivity disorder (ADHD) group. Electroretinograms were recorded in n = 25 ADHD (16 male; age 11.9 ± 2.7 years) and n = 38 (8 male; age 10.4 ± 2.8 years neurotypical control participants as part of a broad study into the determining if the electroretinogram could be a biomarker for ADHD. Flash strengths of 0.6 and 1.2 log cd.s.m^− 2 on a white 40 cd.m^− 2 background were used. Averaged waveforms from each eye and flash strength were analyzed with Bayesian regularization of the covariance matrices using 100 equal length time intervals. The eigenvalues of the covariance matrices were ranked for each group to indicate the degree of complexity within the regularized waveforms.

Results

The correlation matrices indicated less correlation within the waveforms for the ADHD group in time intervals beyond 70 msec. The eigenvalue plots suggest more complexity within the ADHD group compared to the control group. Consideration of the correlation structure between ERG waveforms from different populations may reveal additional features for identifying group differences in clinical populations.

Introduction

The electroretinogram (ERG) assesses retinal function and consists of an initial negative a-wave that is dominated by hyperpolarization of the cone photoreceptors with the subsequent positive b-wave formed by contributions of inner retinal neurons [1,2,3]. The clinical interpretation of the ERG waveform typically involves measures of the amplitudes and time to peaks of the a- and b-waves as well as the higher frequency Oscillatory Potentials and the Photopic Negative Response (PhNR) [4, 5]. Group differences are usually compared using the 95% confidence intervals between groups to compare the differences between the main peaks and troughs of the ERG [5, 6]. This re-analysis of group clinical data uses a Bayesian remodeling of the ERG time-domain trajectory to offer a new perspective on the analysis of the ERG waveform by evaluating correlations between regions of the continuous time series. The ‘trajectory’ is used to denote the amplitude change with time as a vector with magnitude and direction that varies with the shape of the ERG waveform over time.

Analyzing the ERG using signal analysis to decompose the waveform with discrete [7,8,9], continuous wavelet transforms [10, 11], Fourier analysis [12] and variable frequency complex demodulation [13, 14] have been employed to reveal ERG characteristics to inform disease classification [15, 16]. In addition to these methods, that decompose the raw ERG signal into different frequency or scale contributions, other mathematical approaches based on statistical modeling aim to provide insights into the shape of the ERG using a registered time series to analyze the trajectory of the ERG [17]. This later approach drew on aspects of Functional Data Analysis (FDA) where the time series data was treated as an element in the space of possible functions [18]. Broadly speaking, this time-domain analytical approach, which is commonly adopted in fields such as economics, biological signals and meteorology [19] has been neglected in ERG studies. The use of traditional stationary points such as the peak amplitude at a unique (x, y) co-ordinate negates the functional non-stationary y = f(x) properties of the waveform’s trajectory.

Building on the principles of FDA [18], this remodeling of the ERG investigated the correlations within the regularized ERG waveform using detailed Bayesian modeling in the time domain. The ERG waveform was treated as the dependent variable with the ADHD/Control group indicators as explanatory variables to evaluate regions of interest in the ERG waveform shape that differed between the groups. This approach allowed for a broader analysis of the waveform characteristics without relying on a pointwise comparison of the main a- and b-wave peaks. This offers the possibility of identifying more subtle differences between clinical groups or subtypes by examining the correlations of the functions that define the regularized space in the autocorrelation function. In this report the main aim was to demonstrate the Bayesian remodelling approach using recorded ERGs from a control and ADHD group as a exemplar. Classification between the groups is beyond the scope of this preliminary methodological report and will be the focus of future work encompassing other neurodevelopmental disorders.

The remodeling of the ERG waveforms involved the construction of correlation and covariance matrices based on uniform divisions of a 100 × 100 matrix through the waveform. From the rows of the matrix, the interrelationships between data points could be more readily visualized to determine if there were aspects of the waveform that varied between the groups. The approach presented was based on Bayesian analysis using regularization of the (high dimensional) covariance matrices [20], with posterior estimates of both mean ERG time-trajectories and selected features of the covariance or correlation matrices. The posterior estimates were obtained via Hamiltonian Monte Carlo simulations performed in ‘Stan’ [21]. To elucidate uncertainty, the 95% credible intervals for the time-trajectories of ADHD/Control groups were used to identify the time intervals of interest where the credible intervals for the groups did not overlap. The idea of the approach is to provide fine insight into the correlation structure of ERG signal values at different time points. This cannot be done directly (as the sample size is not large enough to estimate full correlation matrix unambiguously). Instead, we employ a fully formalized Bayesian regularization approach to estimate the covariance structure. This analysis provides an example of how a time-domain statistical analysis could complement other existing methods which may be important, especially when acknowledging prominent non-stationarity ERG trajectories and their complicated functional shapes.

Methods

Participants

Participants were recruited through social media and local community psychiatric services for children and adolescents. The sample included n = 25 ADHD (16 male) and n = 38 (8 male) typically developing controls with age (mean ± SD; range) ADHD (11.9 ± 2.7, 6.2–16.2) and control (10.4 ± 2.8, 5.0–15.0) years. Individuals with a diagnosis of ADHD hyperactive/impulsive (n = 23) or inattentive subtype (n = 2) that met DSM-5 [22] diagnostic criteria as confirmed by medical report were included in the sample. Exclusion criteria for both groups included amblyopia, strabismus, retinal dystrophy, epilepsy, refractive error > ± 6.00 equivalent spherical refraction, or systemic disease that may affect retinal function such as diabetes. For the ADHD group participants were also excluded if they had a co-diagnosis of autism spectrum disorder or obsessive-compulsive disorder and for the control group if there was a family history of a neurological disorder such as autism spectrum disorder, schizophrenia or bipolar disorder. Participants (n = 10) abstained from any medication use 24 h prior to testing with methylphenidate that increases dopamine shown to increase the ERG amplitude in some individuals [23]. Medications (with some in combination- see Supplementary Data sheet for details) included methylphenidate (Ritalin or Concerta) (n = 8), melatonin (n = 3), dexamphetamine (n = 2), lisdexamfetamine (Vyvanse) ( (n = 3), alpha-agonists for hypertension (Clonidine, guanfacine or Catapres) (n = 6), D2 antagonists (risperidone) (n = 1), and a selective serotonin reuptake inhibitor (sertraline) (n = 1). The participants dosing regimen varied with symptom severity, but all had been using medications for a minimum of six months prior to testing. Informed written consent to participate was obtained from the parent/guardian prior to testing for individuals aged under 16 years of age or from the participant if aged greater than 16 years of age.

Electrophysiology

ERG recordings were obtained from flash strengths of 0.6 and 1.2 log phot cd.s.m^− 2 presented to the right then left eye of the participant with two replications for each eye/strength combination. One representative waveform was selected from each eye/strength combination based on the quality of the recording (minimal baseline drift or noise) for analysis. The white flash was presented at 2 Hz on a 40 cd.m^− 2 white background with 30 averages obtained to generate the reported averaged ERG waveform with the RETeval (LKC Gaithersburg, MD, USA) device. Skin electrodes were used and placed 2–3 mm below the lower lid. Pupils were undilated with the RETeval automatically adjusting the flash strength to vary with pupil diameter in real time. Waveforms were filtered 0.1–300 Hz and artefacts were automatically rejected if they fell outside the upper or lower quartile of the average. For further information on the recording methods see Constable et al., (2022) [9].

Dataset is available at the public repository: https://figshare.com/s/6e3b93a0866234c6fff4.

Statistical analysis

The ERG signal trajectory was modelled in the time-domain with one functional measurement (ERG profile) per combination of participant, eye (Left or Right) and for each flash strength (0.6 or 1.2 log cd.s.m^− 2). A single measurement profile was the set of time (msec) and the ERG signal amplitude (µV) values. The reported raw times were sampled at 1.95 kHz but had different absolute placements for different measurements from each recording which meant that the data from different measurements were not directly comparable for direct analysis and called for a functional approach (typically, there was a difference in the reported time between recordings of ± 0.01 msec and so there was not an exact match between traces in the time domain). The raw average waveform data was preprocessed by interpolation to the same finely spaced timings for each waveform. This was done by fitting a Generalized Additive Model with time as the explanatory variable using a thin-plate smoother [24] with large k (reduced rank version of the thin plate regression spline) in the mgcv package [25] which were then evaluated at the same dense mesh for each measurement (100 equal lengths placed time values in the interval [-50.24 to 169.34] msec).

The preprocessed functional data were then modelled in the Bayesian way with a multivariate t distribution as the observation model. The multivariate t was fitted separately for the four combinations of eye (Left and Right) and flash strength (0.6 or 1.2 log cd.s.m^− 2). The multivariate t (instead of multivariate normal) was used as a precaution against occasional outlying values. To achieve robustness against the outliers as a part of the observational model, a very low degree of freedom (3) was chosen that still allowed for the existence of second moments [26].

To model the profile flexibly, a flexible covariance matrix was required whose structure was not available a priori. A completely free (general positive definitive) covariance matrix estimation was not feasible here (with 4950 entries, it would have required an impractically large sample size). Therefore, the ERG profile was regularized [27] using the Lewandowski-Kurowicka-Joe (LKJ) distribution [20] prior with the eta parameter equal to 1 (corresponding to uniform selection from the class of legal covariance matrices) for the correlation matrix. This choice was motivated by both flexibility and regularization with respect to the finite information available in the observed data. Other priors (for multivariate mean and for the scaled versions of standard deviations) were uninformative but proper (normal, respectively half normal with large standard deviations).

The Bayesian model was fitted in the probabilistic modelling environment ‘Stan’ [21] using the computationally highly effective Hamiltonian Monte Carlo approach. R was used for data manipulation and pre-processing with the rstan (https://mc-stan.org/rstan/) interface of ‘Stan’. The high-dimensional LKJ prior could not be performed by a default approach in ‘Stan’, so an onion method was used for the implementation (https://discourse.mc-stan.org/t/using-the-onion-method-in-occupancy-hierarchical-model/24901) published on ‘Stan Discourse’.

Results

Figure 1 shows the posterior mean and 95% credible intervals of the ERGs for the ADHD and typically developing control groups for each eye at 0.6 and 1.2 log phot cd.s.m^− 2. The posterior mean is a standard summary of posterior distribution obtained as the centerpiece of Bayesian inference and not the expected value or arithmetic mean of an observable quantity. The main differences based on non-overlapping regions are at time points beyond 70 msec with the control group having a more negative waveform profile. The main group differences appear in the region beyond 70 msec with separation of the credible confidence intervals between the ADHD and Control groups. See the Additional file 1 for the posterior mean and confidence intervals for the amplitudes of the a- and b-waves and additional plots of the correlation and covariance matrices for each group and for each eye and flash strength.

Fig. 1

Posterior mean ERG trajectory for the Control and ADHD groups together with 95% credible intervals obtained from the Bayesian model. For time points beyond 70 msec the Control Groups Waveform is lower than that for the ADHD group with this difference being more pronounced at the 1.2 log cd.s.m^− 2 strength indicated by the boxed area where there is no overlap of the 95% credible confidence intervals in the right and left eyes

Correlation matrix plots for the 100 equal time intervals (See Supplementary Data Sect. 7.1) that spanned the measured ERG signal from − 50.24 to 169.34 msec. Each time interval corresponds to a row number with for example row 37 incorporating data within the time interval (29.60 to 31.82 msec) and row 52 incorporating data within the time interval (62.88 to 65.09 msec). Selected time intervals that encompassed the region including the a-wave (row 29, t = 11.86 to 14.08 msec), b-wave (row 37, 29.60 to 31.82 msec), PhNR (row 56, 71.75 to 73.97 msec) and at row 80, (124.98 to 127.20 msec) where there were noticeable group differences based on the 95% credible confidence intervals indicated by the red boxed areas in Fig. 1. The correlations for these selected rows are shown in Fig. 2 where the correlations are greatest closest to neighboring time points at adjacent rows. As distance from the row increased the correlations between data declined. Here it is important to observe the differences in how these correlations differ between the groups in each row with in general less correlation between data points within the ADHD group compared to the Control group at time points distant to b-wave in rows 51 and 80. Figure 2C and G show correlations at time ~ 72 msec, the waveform data of the ADHD group preceding this time point shows less correlation compared to the Control group. Similarly for Fig. 2D at t = 124.98 msec there is a large difference in the way the data correlates at this point with the neighboring data with the ADHD group showing weaker correlation within the time series. The main discrepancies in the correlations of the time series trajectories are highlighted by the boxed regions shown in Fig. 2.

Fig. 2

Mean posterior correlation (y-axis) against the 100 equal length time divisions (x-axis) of the right eye at the two flash strengths of 0.6 (2 A-D) and 1.2 log phot cd.s.m^− 2 (2E-H) for the Control and ADHD groups. Row 29 (A, E) is at the beginning of the time division beginning at (t = 11.9 msec), for the a-wave with a similar pattern of correlation across the 100-time divisions which was also observed for rows 37 (B, F) at the time interval beginning at t = 29.6 msec for the b-wave. In contrast for row 56 at time division beginning at 71.7 msec and corresponding to the Photopic Negative Response (PhNR), there was less correlation with the preceding time divisions in the ADHD group as indicated by the red boxed areas in Fig. 2C and G. This pattern of a reduced time division correlations in ADHD was more prominent at row 80 at time division beginning at 124.9 msec and is indicated by the purple box in 2D for the 0.6 log phot cd.s.m^− 2 strength. See Additional file 1 for left eye exemplars

Figure 3 illustrates the eigenvalues of the mean posterior correlation matrix of the ERG waveforms data for the right and left eyes comparing the ADHD and Control groups at each flash strength. The decline of the eigenvalues is very fast showing a relatively simple correlation structure of the (prior-regularized) correlation estimate (with only about 30 eigenvalues dominant). To note that whilst the overall shape of the eigenvalue plots were similar for both groups, the Control group plot tends to be larger than the ADHD group for larger eigenvalues that correspond to simpler, more regular behavior. In contrast for the smaller eigenvalues that correspond to more complex and less regular behavior the group plots demonstrating that the ADHD group tends to feature more complex or irregular behavior.

Fig. 3

Ranked eigenvalues of the mean posterior correlation matrices for the right and left eyes of the two groups for each flash strength. Larger eigenvalues are present in the control group in the first few ranks with more frequent higher eigenvalues present in the ADHD group in ranks greater than 50 indicating a more irregular structure in the correlations for the ADHD group and more irregular behavior of the time series

Discussion

This time-domain approach to ERG signal analysis provides a different perspective to the more traditional point-wise comparison of the mean values of amplitude and time to peak of the main features of the ERG waveform [5, 6]. In this approach, a comparison of the correlations of the mean trajectories in the time-domain elaborates where the two groups differ most in terms of mean magnitude of ERG signal, how smooth the behavior of the trajectories is (how much correlated or redundant in time sense they are) and whether the (inherently nonstationary) correlation properties of the ERG signal differ more in some time regions. Both the posterior mean and correlation/covariance behavior differences might be of interest in future constructions of e.g. variants of quadratic discriminant analysis weighting different time regions differently to achieve better classification performance [28,29,30].

The first step in the analysis was to remove the influence of outliers so that the time-domain features of the ERG waveforms could be compared by using time as the explanatory variable. This enabled a comparison of the regions with overlap and non-overlap in the ERG that added additional insights to the group comparisons beyond the a- and b-wave maxima and minima. In this sample, the most distinguishing features between these groups was observed in the time points beyond 70 msec where retinal ganglion cells influence the shape of the ERG [31]. This pattern of differences between these groups may support previous work that found greater background retinal noise associated with the pattern ERG [32].

The main point of comparing the correlations between groups is that they depict where the non-stationary or variable nature of the ERG trajectories are evident which may be used as a potential feature with which to identify differences in the ERG’s trajectory between clinical populations. In this sample, the behavior of the correlations between data points varies for each group and for the ADHD group the correlation/inter-relationships between data points of the ERG are lower than in the Control group. These local variations may be used to infer a difference in the groups and may support classification models in larger samples- or help to identify subsets within a heterogeneous population. These observations serve to illustrate the potential of comparing the degree of correlation of time points with those close and far given that the shape of the waveform is a dynamic interaction arising from the underlying excitatory and inhibitory neural pathways of the retina [1,2,3, 9].

The eigenvalues of the mean posterior correlations explain the main overall multivariate variability which declined rapidly in the groups. This means that less of the variability was explained by the higher order eigenvalues. There is a crossing (at approximately rank 30) of the eigenvalues with the Control group dominating the ADHD eigenvalues in the lower order (first ranks) but the ADHD group dominating the Control eigenvalues for the higher ranks. This hints that the Control group ERG has a more regular, less complex nature whilst the ADHD group ERG is more complex and more irregular.

Limitations

This report is limited to two flash strengths and two groups and a small sample of participants. More detailed analyses to ascertain the sensitivity and specificity of the group difference based on the differences in the correlation matrices and the eigenvalues compared to standard stationary markers would support the use of this approach. Interpretations of the underlying physiological processes involved in the variations between the groups with respect their correlations would also be speculative.

Conclusion

This novel remodeling of the time trajectories of the ERG provides a new perspective on how the correlation structure of the time series behaves. The pre-processing using a Bayesian approach, with a multivariate t distribution as the observation model, enabled explicit regions of interest where the groups waveforms showed time intervals of non-overlap between the credible confidence intervals. These were distant to the main a- and b- wave peaks, that are generally used to compare group differences in neurodevelopmental disorders [33,34,35]. This analysis may be useful when interpreting the rise and fall of the b-wave or slope and turning point of the a-wave to gain insights into the structure of the waveform before and after the main peaks that are shaped by underlying neural pathways [1,2,3, 31] and provide information about the location of significant fluctuations in the shape of the ERG waveform to inform interpretations of the contributions to the waveform in retinal and neurological disorders [36, 37]. This approach could also benefit analysis of waveforms whose nature may be variable owing to cortical topography such as visual evoked potentials to identify regions of interest in clinical populations [38]. Similar analyses may also benefit changes in the pattern ERG which shows early changes in glaucoma that may be analyzed using a Bayesian remodeling approach outlined in this report [31, 39].