Functional decomposition of the human ERG based on the discrete wavelet transform

Abstract

The morphology of the electroretinogram (ERG) can be altered as a result of normal and pathological processes of the retina. However, given that the ERG is almost solely assessed in terms of its amplitude and timing, defining the shape of the ERG waveform so that subtle, physiologically driven, morphological changes can be systematically and reproducibly detected remains a challenging problem. We examined if the discrete wavelet transform (DWT) could meet this challenge. Normal human photopic ERGs evoked to a broad range of luminance intensities (to yield waveforms of various shapes, amplitudes, and timings) were analyzed using DWT descriptors of the ERG. Luminance-response curves that were generated using the various DWT descriptors revealed distinct (p < 0.05) luminance-dependence patterns, indicating that the stimulus luminance differently modulates the various time-frequency components of the ERG and thus its morphology. The latter represents the first attempt to study the luminance-dependence of ERG descriptors obtained with the DWT. Analyses of ERGs obtained from patients affected with ON or OFF retinal pathway anomalies were also presented. We show here for the first time that distinct time-frequency descriptors can be specifically associated to the function of the ON and OFF cone pathway. Therefore, in this study, the DWT revealed reproducible, physiologically meaningful and diagnostically relevant descriptors of the ERG over a wide range of signal amplitudes and morphologies. The DWT analysis thus represents a valuable addition to the electrophysiologist's armamentarium that will improve the quantification and interpretation of normal and pathological ERG responses.

Introduction

The electroretinogram (ERG) represents the biopotential generated by the retina in response to a light stimulus. To date, the ERG remains the only mean to objectively and noninvasively assess the functional integrity of the retina. Although interpretation of the ERG relies mostly on peak-time and amplitude measurements (i.e., usual time domain [TD] measurements) of the a- and b-waves (McCulloch et al., 2015), it is well known that retinal disorders can also alter the morphology of the signal, a feature of the ERG that is not objectively measured. This is best exemplified with the ERG of patients affected with congenital stationary night blindness (CSNB) which is often qualitatively reported as having a square-wave–like a-wave and truncated b-wave morphology (Heckenlively, Martin, & Rosenbaum, 1983; Lachapelle, Little, & Polomeno, 1983; Miyake et al., 1997). This morphological feature is pathognomonic for CSNB, but less obvious changes of the ERG morphology are likely to be unnoticed and/or ignored due to the nonexistence of more objective ERG shape–dependent measurements. Consequently, the present work addresses the following question: Is it possible to objectively quantify ERG signals of various shapes so that, physiologically driven, morphological changes can be systematically and reproducibly detected and used as diagnostic clues? 

Intuitively, the morphology of a given ERG waveform is most likely determined by how its different frequency components combine, in energy and latency, to generate the waveform. While different frequency components of the ERG signal can be identified using Fourier analysis (such as the fast Fourier transform [FFT]), they cannot be localized in time because the FFT does not include a temporal resolution (Breslin & Parker, 1973; Gur & Zeevi, 1980; Poppele & Maffei, 1967). Both spectral and temporal resolutions are needed to accurately define the properties (e.g., morphology) of transient time-varying signals (such as the ERG) and, fortunately, time and frequency resolutions can be optimally obtained with the use of continuous (CWT) or discrete (DWT) wavelet transforms (Mallat, 2009). The growing use of these modern approaches in the fields of auditory, somatosensory, and event-related evoked potentials has led to significant diagnostic improvements. For example, early work by Thakor, Guo, Sun, and Hanley (1993) revealed that the wavelet decomposition was well suited to characterize and rapidly detect complex changes in the shape of somatosensory evoked-potentials that resulted from cerebral hypoxic injury. Similarly, Demiralp and Ademoglu (2001) also went beyond traditional TD amplitude and latency measurements and beyond Fourier analysis by processing an extensive characterization of event-related potential (ERP) signal morphologies through the use of the wavelet transform. They also demonstrated that the waveform morphology of the ERP signal was carrying important information that could be segmented in multiple functional components by the wavelet decomposition (Demiralp & Ademoglu, 2001). In another study, Quiroga and Garcia (2003) improved the detection of single-trial ERPs using a wavelet-denoising approach that selectively reconstructs the most meaningful DWT coefficients within precise frequency bands and time windows (Quiroga & Garcia, 2003). Given the better performances of their algorithm, compared to more common filtering approaches, their technique subsequently inspired several investigators. 

In contrast, fewer attempts have been made to evaluate the usefulness of wavelet techniques to the ERG. So far, wavelets were sporadically applied to the analysis of pattern ERG (PERG), multifocal ERG (mfERG) and full-field scotopic ERGs of human and rats. In rats, CWT analysis revealed that two frequency components (i.e., 70–80 Hz and 120–130 Hz) contributed to the genesis of the oscillatory potentials (OPs; Forte, Bui, & Vingrys, 2008). This study also showed that the energy level, the frequency and the latency of the two OP frequency bands had distinct luminance-response (LR) functions, suggesting that they might be evoked by different retinal elements/mechanisms. The CWT also revealed that the scotopic a-wave of normal human subjects was composed of three frequency components (20, 140, and 180 Hz) and that the higher frequency component was absent in achromates, suggesting that the CWT could also help in the diagnosis of photoreceptoral diseases (Barraco, Persano Adorno, & Brai, 2011). Moreover, DWT coefficients were shown to be superior to traditional TD measures in segregating normal and pathological PERG waveforms using principal components analysis (Rogala & Brykalski, 2005). It was also shown that it is possible to approximate the TD a-wave, b-wave, and at least one OP using the inverse DWT of the scotopic ERG waveform (Varadharajan, Fitzgerald, & Lakshminarayanan, 2007). Similarly, Miguel-Jiménez, Boquete, Ortega, Rodriguez-Ascariz, and Blanco (2010) used the DWT decomposition to reconstruct mfERG waveforms into different frequency bands and suggested that this approach could potentially be more sensitive to detect changes in glaucoma patients, compared to Humphrey visual field tests. As indicated above, wavelet analysis can reveal subtle (and possibly diagnostic) ERG changes that are almost impossible to appreciate with TD measures. 

Of interest, we previously compared the FFT, CWT, and DWT with the more traditional TD measures (amplitude and peak time of the a- and b-waves) in their abilities to analyze the ERG signal and concluded that the DWT offered significant advantages (Gauvin, Lina, & Lachapelle, 2014). With its simplified scalograms and its predetermined time-frequency borders assigned to each ERG component, the DWT significantly eases the identification of relevant descriptors and highlights additional frequency components. However, in the latter study of Gauvin et al. (2014), analysis was limited to the suprathreshold photopic ERG response, which was previously shown to be the most complete (optimal) cone-mediated ERG response (Lachapelle, Rufiange, & Dembinska, 2001). Of course, it is also well known that flash intensity and/or pathology may reduce or remove some of the components, which when added together, make this optimal photopic ERG signal, an aspect of photopic ERG genesis that was not explored in our first paper, thus preventing us from evaluating if the use of DWT descriptors can be generalized to detect changes in ERGs, irrespective of their morphology, amplitude, and timing. 

As formerly demonstrated, the LR function of the human photopic ERG b-wave is rather unique in that, with a progressive increase in stimulus luminance, the amplitude of the b-wave first increases, reaches a maximal value and then decreases before reaching a plateau with the brightest intensities. This change in amplitude is always accompanied by variations in the overall ERG waveform morphology. This phenomenon, first described by Wali and Leguire (1992) as the photopic hill (PH), was previously explored by the current authors (Garon et al., 2010; Rufiange et al., 2003; Rufiange, Dumont, & Lachapelle, 2005; Rufiange, Rousseau, Dembinska, & Lachapelle, 2002) and others (Hamilton, Bees, Chaplin, & McCulloch, 2007; Kondo et al., 2000; Ueno, Kondo, Niwa, Terasaki, & Miyake, 2004). Normal ERG responses evoked to 21 stimulus intensities that covered a range of 5 log-units of intensity (to generate the PH) were therefore used herein as a method to modulate the morphology, amplitude, and timing of the ERG response. 

The current study aims to characterize and compare the luminance-dependence of the various DWT descriptors that were used to analyze the morphologically different ERGs known to compose the normal PH. We hypothesize that the ERG morphology will be captured by the DWT descriptors and that PH-like LR functions will also characterize LR functions obtained using DWT descriptors of the ERG that are concomitant with the b-wave. Moreover, we investigated the complexity of the ERG waveform using new DWT descriptors that we derived from wavelet variance analyses (WVAs). Finally, to assess if DWT descriptors appraise meaningful physiological attributes of the cone ERG, analyses of abnormal responses obtained from patients affected with known retinal pathway anomalies are also presented. 

Methods

Analyses were conducted on a total of 25 (16 women and nine men) normal subjects (30.5 ± 8.1 years old) who had signed an informed consent form previously approved by the Institutional Review Board of the Montreal Children's Hospital. Experiments were conducted in accordance with the Declaration of Helsinki. This study was limited to the photopic ERG for the uniqueness of its LR function (the so-called PH), with its four characteristic phases (i.e., ascent, maximal value, descent, and final plateau phases) and ensuing different cone ERG morphologies (Garon et al., 2010; Hamilton et al., 2007; Kondo et al., 2000; Rufiange et al., 2003; Rufiange et al., 2005; Rufiange et al., 2002; Ueno et al., 2004; Wali & Leguire, 1992). 

Preparation of subjects and ERG recordings

According to a previously used protocol of ours (described in Garon et al., 2010; Rufiange et al., 2003; Rufiange et al., 2005; Rufiange et al., 2002), ERG signals were recorded (LKC UTAS-E-3000 system; LKC Technologies Inc., Gaithersburg, MD) with both eyes dilated (Tropicamide 1%) using a Dawson, Trick, and Litzkow (DTL) fiber electrode (27/7 X-Static silver-coated conductive nylon yarn, Sauquoit Industries, Scranton, PA) placed in the inferior conjunctival bag, with reference and ground electrodes (Grass gold-plated cup electrodes filled with Grass EC2 electrode cream; Grass Technologies, Warwick, RI) pasted at the external canthi and forehead, respectively. Photopic ERGs (bandwidth: 1–500 Hz; amplification: 20,000×; attenuation: 6 dB; sampling frequency: 3413.33 Hz) were recorded against a broadband white (color temperature: 6500 K), rod-desensitizing, background light of 30 cd.m⁻². Photopic hills were obtained in response to graded intensities of stimulation (ranging between −0.8 and 2.64 log cd.s.m⁻² in 14 steps of ∼0.26 log-units; average of 10 flashes per intensity; flash duration: 20 μs; white light; interstimulus interval: 1.5 s; prestimulus baseline: 20 ms) in 15 of the 25 subjects. ERGs evoked to seven dimmer flash intensities (ranging between −2.23 and −1.00 log cd.s.m⁻² in 0.2 log-unit steps; average of 50 to 300 flashes; flash duration: 20 μs; white light; interstimulus interval: 0.3 s; prestimulus baseline: 20 ms) were obtained from the other 10 subjects. Background light and integrated flash luminances were measured with a research radiometer (IL1700; International Light, Newburyport, MA). ERGs from both eyes were averaged to yield a single waveform and imported in MATLAB R2014a software (Mathworks, Natick, MA) for further analyses. 

Selection of pathological ERGs

Depending on the nature of the disease process, retinal disorders can alter the amplitude and/or the peak time of the photopic and/or scotopic ERGs. Moreover, retinopathies can also remarkably affect the overall morphology of the ERG waveform. For the sake of the current work, we chose to limit our analysis of pathological signals to photopic (flash intensity: 0.64 log cd.s.m⁻²; rod-desensitizing background light: 30 cd.m⁻²) ERGs obtained from patients (n = 20) affected with Type-1 CSNB (n = 10) or congenital postreceptoral cone pathway anomaly (CPCPA; n = 10). These patients were selected for two main reasons. Firstly, the selected ERG waveforms show, upon visual inspection, strikingly different morphological features that, we claim, will be captured by the DWT descriptors. Secondly, the functional anomaly of patients affected with CSNB is known to specifically reside with the ON-pathway (based on electroretinographic (Langrova et al., 2002; Miyake, Yagasaki, Horiguchi, & Kawase, 1987; Quigley et al., 1996) and molecular (Bech-Hansen et al., 2000; Dryja et al., 2005; Gregg et al., 2007; Pusch et al., 2000) findings and complementary, the functional anomaly of patients affected with CPCPA is believed to lie on the OFF-pathway (Garon et al., 2014; Lachapelle et al., 1998). We hypothesize that retinal conditions affecting the ON (CSNB) and OFF (CPCPA) pathways will affect different DWT descriptors. 

TD analysis of the ERG

The amplitude of the a-wave was measured from the baseline to the trough of the a-wave and that of the b-wave from the trough of the a-wave to the peak of the b-wave. The amplitude of each OP was measured from trough to peak and summed to obtain the sum of OPs (SOPs). We also computed the signal-to-noise ratio (SNR) of selected ERGs as the ratio of the b-wave amplitude over the amplitude of the baseline noise (i.e., maximal peak-to-peak amplitude measured within the first 20 ms of recording prior to flash onset). 

Computation of the DWT

The DWTs of all ERG responses were obtained using the Matlab R2014a application (computation detailed in Appendix A). For each DWT, we displayed eight levels (i.e., eight frequency bands) of decomposition, each with a distinct central frequency (Level 1 = 1280 Hz; Level 2 = 640 Hz; Level 3 = 320 Hz; Level 4 = 160 Hz; Level 5 = 80 Hz; Level 6 = 40 Hz; Level 7 = 20 Hz; and Level 8 = 10 Hz; see Figure 1A), where each band quantified the contribution of a range of components oscillating around the central frequency (for example: 20 Hz = 20 ± 20/3 Hz, 40 Hz = 40 ± 40/3 Hz, and so on). 

Figure 1

 

Identification of selected DWT descriptors. (A) DWT scalogram of a normal photopic ERG (B) evoked to a 0.64 log cd.s.m⁻² stimulus. The most prominent oscillatory component of the ERG appears as the darkest red rectangle (highest energy; see the calibration color bar) in the region of the scalogram where it is located (here centered at 40 Hz [frequency] and between 17.5 and 36.25 ms [time]), while the absence of oscillating components, at any given location, appears as the darkest blue rectangle (no energy). Almost an infinite number of descriptors can be calculated from the DWT scalogram. For example, the six LWM descriptors included in this study are identified (20a, 40a, 20b, 40b, 80ops, and 160ops) and delimited by white borders at (A). As shown, some LWM descriptors include a single wavelet coefficient (i.e., rectangle) within their borders, such as 20a, 40a and 20b, while others comprise several coefficients, such as 40b (two coefficients), 80ops (five coefficients), or 160ops (10 coefficients). For quantification purposes, LWM descriptors that included more than one coefficient were defined as the maximal or averaged value of included coefficients (see Table 1 for details). (B) Normal photopic ERG evoked to a 0.64 log cd.s.m⁻² stimulus. “a” and “b” identify the a- and b-waves and “OPs” identify two oscillatory potentials (OP2 and OP3). (C) The original ERG waveform (blue curves) is accurately (Pearson's score = 0.9853) reconstructed (i.e., inverse DWT; black curves) by summating the four DWT levels that included our LWM descriptors (20, 40, 80, and 160 Hz). (D–E) Wavelet variance analyses also allows for the calculation of robust DWT descriptors. For example, the blue curve of (D) represents the variance of the wavelet coefficients (i.e., expressed as the standard deviation) measured at each of the eight DWT levels (identified as Levels 1 to 8 in the DWT scalogram of [A]) and the red line represents the linear fit (R² = 0.93) between Levels 2 and 5. The difference between the variance measured at Level 6 (50.13) and its prediction from the linear fitting (17.6) defines the Δ-variance descriptor (32.53 in this example). Furthermore, the natural logarithm of the wavelet variance data points shown in (E) was fitted by a linear model between Levels 2 and 6 (as represented by the red line; R² = 0.91). The Hölder exponent (1.19 in this example) is defined as the slope of this linear fitting.

Table 1

 

Definition of the local wavelet maxima (LWM) descriptors. Notes: Column 1: The name by which the LWM is referred to in the text. Column 2: Time window within which the descriptor is localized. Column 3: The frequency band of the LWM. Column 4: How the descriptor is computed. Note that for 160ops, the first maximal coefficient is computed by taking the maximal value of the first two coefficients, the second maximal value is computed by taking the maximal value of the third and fourth coefficients, and so on until the fifth maximal value, which is computed by taking the maximal value of the last two coefficients. Column 5: Extent of the translations of the electroretinogram (ERG) that is used to optimally localize each LWM descriptor. Negative and positive values indicate a translation of the ERG to the right and to the left, respectively. Translations are achieved with increments of 1 ms.

 

DWT quantification of the ERG waveform with local wavelet maxima descriptors

DWT scalograms were obtained by plotting the absolute value of the wavelet coefficients on a dyadic time-frequency grid. For this study, we have considered the following six LWM descriptors (identified in Figure 1A and defined in Table 1), namely 20a, 20b, 40a, 40b, 80ops, and 160ops (see Appendix A for details on how they were derived). As indicated in Table 1 (and illustrated in Figure 1A), each LWM descriptor quantified the energy of local oscillations of the ERG signal within a precise and predetermined time-frequency window. Due to their temporal positions (compare Figure 1A, B), we previously suggested (Gauvin et al., 2014) that the 20a, 40a, 20b, and 40b DWT descriptors were most probably associated with the a- and b-waves of the ERG signal, respectively. This claim will be further investigated herein by assessing if strong correlations exist between the LR functions of these DWT descriptors and that of the a- and b-wave amplitude, respectively. Similarly, due to their frequencies (80–160 Hz) and temporal positions, the 80ops and 160ops descriptors were associated to the OPs. In this study, we needed the OPs descriptors to quantify all OPs that were included in a given ERG waveform, irrespective of stimulus intensity or health status of the retina. Consequently, as shown in Table 1, we computed the 80ops and 160ops descriptors as the average of five coefficients in order to obtain a global measurement of the OPs. Note that the white borders delimiting the 160ops descriptor (Figure 1A) included 10 coefficients so, in order to have an unbiased average (i.e., to have the same number of coefficients to average for both the 80ops and 160ops), we computed the 160ops as the average of five coefficients (see details in Table 1 and its caption). As shown in Figure 1C, the summation of the four frequency bands (i.e., 20, 40, 80, and 160 Hz) that included the six LWM descriptors considered in the present study allows us to reconstruct a synthetic ERG waveform (i.e., inverse wavelet transform; last black tracing of Figure 1C) that explain 98.53% (as per Pearson coefficient) of the original ERG (blue tracing), indicating that the time-frequency components that we identified are the most (if not the only) important contributors to the genesis of the photopic ERG waveform. Finally, given that variation in the peak times of the different ERG components is expected and that the DWT is not a shift-invariant transform (Guo, 1995), we calculated each LWM as the maximum value obtained while shifting (i.e., translating) the complete ERG waveform to the left and right directions of the time axis. As reported in Table 1 (column 5), the magnitude of the translation was limited to a given range to the left and right (i.e., corresponding to positive and negative translation values in Table 1). The range of the translation was selected by trial and error in order to conservatively cover the expected variation of photopic ERG peak times (i.e., ∼ 5 ms/log-unit increment of the stimulus intensity, as estimated from Garon et al., 2010). The above-mentioned translations thus allowed an optimal (i.e., maximal) measurement of the LWM more independently of ERG peak times. Full translation details and demonstrations are reported in Appendix A. 

DWT quantification of the ERG waveform with WVA

WVA represents the variance of all wavelet coefficients obtained at each level (i.e., frequency band) of a DWT (Percival, 1995). As shown in Figure 1D, WVAs are reported as a plot correlating the variance values (expressed as the standard deviation in this study) with the corresponding DWT frequency level (Gallegati, 2008; Park & Willinger, 2000). As exemplified, the variance of a typical ERG evoked at 0.64 log cd.s.m⁻² starts at near-zero values (Levels 1 to 3), rapidly increases (steepest increment between Levels 5 and 6) to reach a maximum at Level 7, which is then followed by a final decrease in the value of the wavelet variance. This graphic representation allowed us to identify the delta-variance (Δ-variance) descriptor, the computation of which is presented in Appendix A and illustrated in Figure 1D. The Hölder exponent (also termed scaling exponent) represents another descriptor that can be computed using the WVA of DWT (Abry, Flandrin, Taqqu, & Veitch, 2002). This descriptor, which is amplitude-independent, characterizes the irregularity (or roughness) of the signal, a feature often associated to the complexity of a given waveform (Bishop, Yarham, Navapurkar, Menon, & Ercole, 2012; Sen, Litak, Kaminski, & Wendeker, 2008). The computation procedure of the Hölder exponent is detailed in Appendix A and illustrated in Figure 1E. 

Selection of mother wavelets

The selection of a wavelet is mainly related to the choice of an optimal balance between temporal accuracy and spectral resolution to represent the fluctuations in the signal. This trade-off is controlled by the number of vanishing moments of the wavelet. High number of vanishing moments implies a better spectral resolution, but lesser temporal accuracy and vice-versa (Strang & Nguyen, 1996). In this study, we defined the local wavelet maxima (LWM) descriptors of the ERG to analyze features of the signal that are well-localized in time (such as the a- and b-wave features). With zero vanishing moment, the Haar wavelet provides the best temporal accuracy and hence, the Haar wavelet coefficients have the ability to optimally identify local features of a signal (Daubechies, 1992). Furthermore, the Haar wavelet has the shortest wavelet filter length, making its coefficients less sensitive to cross-contamination from neighboring oscillations (such as that of the b-wave on the a-wave). Given the above properties, we opted for the Haar wavelet to analyze the LWM descriptors of the ERG. In contrast, the WVA descriptors analyze global scaling characteristics of the signal. The symmetric Daubechies wavelet (Lina & Mayrand, 1995) set with two vanishing moments (or SDW2) also has a compact support, but because of its shape and added vanishing moments (i.e., better spectral resolution), it has the property to specifically encompass the entire ERG waveform. This makes this wavelet blind to low-order polynomial trends that are nonspecific to the ERG waveform, thus leading to a sparse representation of the genuine fluctuations present in the signal. The latter allows a more robust estimation of wavelet variance at each level of the DWT (Abry et al., 2002). The SDW2 wavelet was thus selected to quantify the WVA descriptors. 

Statistical analysis

Statistical analysis of the LWM descriptors was performed using two-way within-subjects analysis of variance (ANOVA; between −0.8 to 2.64 log cd.s.m⁻²; n = 15 subjects) followed by post hoc Bonferroni-paired t tests for multiple hypotheses testing with repeated measures. Statistical analysis of the wavelet variance descriptors (Hölder and Δ-variance) was obtained using one-way (stimulus intensity) within-subjects ANOVAs (between −0.8 to 2.64 log cd.s.m⁻²; n = 15 subjects). Multiple comparisons (based on Bonferroni-paired t test for repeated measurements) were used to compare the values of selected pairs of means obtained under different stimulus luminances. The pathological ERG groups (CSNB and CPCPA) were compared to control and between them using unpaired two-sample t tests. The coefficient of variation (CV) of selected group data was computed as the standard deviation divided by the mean and multiplied by 100. Finally, Pearson correlation coefficients (i.e., r for rho) were used to compare the similarities between ERG waveforms or between LR functions. The significance level of each test was fixed at 0.05. 

Results

Figure 2A illustrates representative photopic ERG responses, evoked to progressively brighter stimuli, which were used, in the present study, in order to generate the typical PH (Figure 2B). The corresponding DWT scalograms are shown in Figure 2C. 

Figure 2

 

Representative ERG responses and associated DWT scalograms. (A) Composite ERG responses (i.e., obtained from arithmetic average of ERG responses obtained from 15 normal subjects) evoked to each of the 14 progressively brighter stimuli. “a” and “b” identify the a- and b-waves of the ERG and the vertical arrow indicates stimulus onset. (B) Luminance-response function of the b-wave amplitude obtained from the ERGs shown in (A). The black vertical line identifies the stimulus intensity (0.39 log cd.s.m⁻²), which yielded maximal b-wave amplitude. (C) DWT scalograms of the 14 ERGs shown in (A). The LWM descriptors considered in this study (20a, 40a, 20b, 40b, 80op, 160ops) are identified in the 0.39 scalogram where they are delimited with white borders. For each scalogram, the color-coding was normalized to the LWM of maximal energy (see color bar on the right-hand side of the bottom-right scalogram). The frequency bands are indicated on the right-hand side of the top-right scalogram. White arrows indicate delayed positions of the 40b descriptors compared to that of the previous stimulus intensities. (D) Luminance-response function of the total energy included within the white borders (i.e., ∑Energy: 20a + 40a + 20b + 40b + 80ops + 160ops). The black vertical line identifies the stimulus intensity (0.39 log cd.s.m⁻²), which yielded the DWT with maximal energy content.

The LWM descriptors considered in the present study (identified in the scalogram representing the ERG signal evoked to the 0.39 log cd.s.m⁻² stimulus in Figure 2C and delimited with white borders in the other scalograms) show that, irrespective of the intensity of the stimulus, the maximal energy of the signal is always contained within a time-frequency region delimited by these descriptors. Therefore, irrespective of the amplitude of the ERG signal, LWM associated with the a-wave (20a, 40a), b-wave (20b, 40b) and OPs (80ops, 160ops) remained quantifiable (as per color scale). Of interest, as shown in Figure 2D, when the total LWM energy (i.e., ΣEnergy = 20a + 40a + 20b + 40b + 80ops + 160ops) of each ERG waveform is plotted against the intensity of the stimulus, the shape of the resulting function is reminiscent of the PH obtained when only the amplitude of the b-wave is considered (as in Figure 2B)—both functions reaching their maximal values with the ERG evoked to the 0.39 log cd.s.m⁻² stimulation. This confirms that, as hypothesized, the PH-like shape (originally evidenced with TD measures of the b-wave; blue curve of Figure 2B) continues to remain a signature feature of the cone ERG LR function when the ERG response is quantified in the time-frequency domain. 

Luminance-dependence of the LWM descriptors

LWM descriptors were computed using the Haar wavelet. Figure 3A and B report the mean (±1 SD) photopic a- and b-wave LR functions obtained using the LWM descriptors weighing the a-wave (20a, 40a; Figure 3A) and b-wave (20b, 40b; Figure 3B) energy levels. As seen in Figure 3A, the 20a and 40a descriptors follow distinct LR patterns; the 20a descriptor first increases slowly (from −0.8 to 0.39 log cd.s.m⁻²) and then more abruptly (from 0.39 to 1.4 log cd.s.m⁻²) before reaching a plateau with the four brightest intensities. In contrast, the 40a descriptor appears to follow a logistic-like growth function. As shown in Figure 3D, the growth pattern obtained by summating the 20a and 40a DWT descriptors is almost identical to that obtained using TD measurements of the a-wave (Pearson's correlation coefficient = 0.9983). A similar decomposition of the b-wave (Figure 3B) reveals that the LR functions of the two frequency components (i.e., 20b and 40b) follow distinct PH-like shapes. First, the peak of the 20b is flat compared to the sharper peak of the 40b LR function, and second, while the maximal energy value for the 20b descriptor is reached with the 0.64 log cd.s.m⁻² stimulus, that of the 40b component is attained at 0.39 log cd.s.m⁻². Again, as shown in Figure 3E, the growth function obtained by summating the 20b and 40b DWT descriptors is nearly identical (Pearson's correlation coefficient = 0.9988) to that obtained using the TD measurements of the b-wave (the so-called PH shown in Figure 2B). Statistical analysis revealed that the interaction effects (frequency band × stimulus intensity) were significant for both the a-wave (df = 13; F = 37.16; p < 0.00001) and b-wave (df = 13; F = 27.28; p < 0.00001). Of interest, post hoc analysis revealed that, depending on the intensity of the stimulus used, the energy level concealed in the 40 Hz a-wave parameter (40a) was either significantly higher (p < 0.05) or equal (p > 0.05) to that concealed in the 20 Hz a-wave parameter (20a). Similarly, the 40-Hz b-wave parameter (40b) was either significantly higher (p < 0.05), equal (p > 0.05), or significantly lower (p < 0.05) than the 20-Hz b-wave parameter (20b). These significant intensity-dependent differences (indicated by the black asterisks in Figure 3A, B) in the energy level of the LWM descriptors indicate that the intensity of the stimulus significantly modulates the time-frequency composition of the ERG. As shown in Figure 3C, a similar PH pattern was also obtained with the higher frequency components of the ERG (80ops and 160ops), both attaining maximal values with the 0.39 log cd.s.m⁻² stimulus. Statistical analysis revealed that the interaction effect (frequency band × stimulus intensity) was significant (df = 13; F = 28.98; p < 0.00001). Post hoc tests indicated that irrespective of the intensity, the magnitude of 80 Hz OPs energy descriptor (80ops) was significantly higher (p < 0.05) than the 160 Hz (160ops) one. Finally, as shown in Figure 3F, the LR function obtained by summating the 80ops and 160ops descriptors covaried (Pearson's correlation coefficient = 0.9831) with that obtained using the TD measurement of the SOPs. 

Figure 3

 

Luminance-dependence of LWM descriptors. (A–B) Mean (±1 SD; only one direction shown for clarity) LR functions obtained using the LWM descriptors of the a-waves (i.e., 20a and 40a) and b-waves (i.e., 20b and 40b), respectively. Black asterisks located above and under the curves indicate when the 40 Hz descriptors are significantly (p < 0.05) higher or lower than those of the 20 Hz, respectively. (C) Mean (±1 SD) LR functions obtained using the LWM descriptors of the OPs (i.e., 80ops and 160ops). The black asterisk in the middle of the black line located under the curves indicates that irrespective of the intensity, the 80 Hz OPs energy (80ops) was significantly larger (p < 0.05) than the 160 Hz (160ops). (D–E) The normalized summation of the mean 20-Hz and mean 40-Hz LR functions of the a- and b-wave obtained using the LWM (black curves) are superimposed on the normalized mean a- and b-wave amplitude LR functions measured using the TD approach (green curves). (F) The normalized summation of the mean 80-Hz and mean 160-Hz LR functions (black curve) of the OPs is superimposed on the normalized LR function of the SOPs amplitude measured using the TD approach (green curve).

Luminance-dependence of the WVA descriptors

WVA descriptors were computed using the symmetric Daubechies wavelet. Represented in Figure 4A is the mean (±1 SD) LR function of the Hölder exponent. With progressively brighter stimuli, the value of the Hölder exponent first increases, reaches a maximal value (e.g., maximal complexity of the ERG waveform) at 0.64 log cd.s.m⁻², which is then followed by a gradual decrease with brighter stimuli. The main effect of the stimulus intensity was found to be significant (df = 13, F = 51.26; p < 0.00001). Of note, though the Hölder exponent is amplitude-independent, its LR function is reminiscent of the PH obtained with TD measure of b-wave amplitude, which is shown in Figure 4C. There are, however, important differences, such as (a) the peak value is reached at slightly brighter stimulus intensity (i.e., 0.64 compared to 0.39 log cd.s.m⁻²) and (b) the absence of a plateau effect with the brightest stimuli (compare Figure 4A, C). 

Figure 4

 

Luminance-dependence of WVA descriptors. (A) Mean (±1 SD; only one direction shown) LR functions of the Hölder exponent. (B) Mean (±1 SD) LR function of the Δ-variance descriptor. (C) Traditional LR function of the b-wave amplitude measured in the TD. Black vertical lines indicate the stimulus intensity (0.39 or 0.64 log cd.s.m⁻²), which yielded the maximal value.

Similarly, as shown in Figure 4B, the LR function of the Δ-variance parameter also adopts a PH-like pattern that reaches its maximal value with the 0.64 log cd.s.m⁻² flash and does not plateau with the brightest stimuli. This main effect of the stimulus intensity was also found to be significant (df = 13, F = 22.54; p < 0.00001). 

Of interest, although the LR functions of the two WVA descriptors (Hölder and Δ-variance) reached their peak value with the ERG evoked to the 0.64 log cd.s.m⁻² stimulus, the shapes of the resulting LR functions differed significantly; the peak of the Hölder exponent function being smoother compared to the sharper peak of the Δ-variance function. Supportive of the latter, post hoc analysis revealed that Hölder values immediately adjacent (i.e., evoked at 0.39 and 0.9 log cd.s.m⁻²) to the peak value (0.64 log cd.s.m⁻²) were not significantly different from this peak value (p > 0.05), while the peak value of the Δ-variance descriptor was significantly larger (p < 0.05) than that of the values from the two neighboring intensities. 

DWT analysis of ERGs of lower amplitudes

Figure 5A illustrates photopic ERGs evoked to flash intensities ranging between −2.23 and −1.0 log cd.s.m⁻², while the corresponding DWT scalograms are presented at the right-hand side of each waveform. As shown in Figure 5B, the amplitude of the a- and b-waves thus obtained grows progressively (a-wave: from 0.6 to 5.0 μV; b-wave: from 0.9 to 10.75 μV) with brighter stimuli to reach, in response to the −1 log cd.s.m⁻² stimulus, values that are slightly lower (a-wave: 5.0 μV; b-wave: 10.75 μV) than values obtained at −0.8 log cd.s.m⁻² (a-wave: 7.41; b-wave: 18.24 μV; e.g., Figure 2A, B). Similar to what was shown with the larger amplitude ERGs (see Figure 2C), DWT scalograms of the low-voltage ERGs also reveal that the maximal energy of the ERG signal remains localized in the time-frequency region that is delimited with the six LWM descriptors (white borders) defined above (i.e., 20a, 20b, 40a, 40b, 80ops, 160ops). Similarly, each DWT descriptor (Figure 5C through F) grows progressively with brighter stimuli to reach, in response to the −1 log cd.s.m⁻² stimulus a value slightly lower than the value obtained at −0.8 log cd.s.m⁻² (see Figures 3 and 4), confirming that DWT descriptors can be used to monitor the ERG over a wide range of amplitudes (i.e., in this study: from less than 1 μV to more than 100 μV). 

Figure 5

 

Analysis of low-amplitude ERGs. (A) Composite ERG responses (i.e., obtained from arithmetic average of ERGs obtained from 10 subjects) evoked to seven progressively dimmer intensities (in log cd.s.m⁻² given at the left of each tracings). The a- and b-waves are indicated as “a” and “b”; the vertical arrow (at top of first tracing) indicates stimulus onset. The horizontal and vertical calibration bars apply to each ERG, but for clarity, responses evoked at −1.81, −2.04 and −2.23 have been magnified 2, 4, and 8 times, respectively. Associated DWT scalograms in which the LWM descriptors have been surrounded by white borders are shown on the right-hand side of the waveforms. For each scalogram, the color-coding was normalized to the LWM of maximal energy. The frequency levels and the color bar are indicated on the top and bottom scalogram, respectively. (B) LR functions of the a- and b-waves' amplitude obtained from the ERG shown in (A). (C) Mean (±1 SD; only one direction shown for clarity) LR functions obtained using the DWT descriptors of the a-waves (i.e., 20a and 40a) and b-waves (i.e., 20b and 40b). (D) Mean (±1 SD; only one direction shown) LR functions obtained using the DWT descriptors of the OPs (i.e., 80ops and 160ops). (E–F) Mean (±1 SD; only one direction shown) LR functions obtained using the Hölder exponent and Δ-variance, respectively.

In addition, although intuitively one would expect more variability in the quantification of the different ERG parameters as the ERG response becomes noisier (i.e., lower SNRs; compare Figure 6A, B), it seems that this variability impacts more on TD parameters compared to DWT ones. This is best illustrated with the data shown in Figure 6C, comparing the CV of the a- and b-wave amplitudes measured using the TD approach (TD a-wave and TD b-wave) with their equivalent measurements obtained using the DWT (DWT 20a + 40a; DWT 20b + 40b) for ERGs evoked to −2.23 log cd.s.m⁻² (mean SNR in the TD: 1.62 ± 0.82) and 0.64 log cd.s.m⁻² (mean SNR in the TD: 19.65 ± 6.77). A 92% reduction in the SNR values, increased (i.e., seen from right to left in Figure 6C) the CV of the a- and b-wave measurements by 125% and 98%, respectively, using the TD approach, compared to 50% and 31% using the DWT approach. 

Figure 6

 

Analysis of ERGs with low and high SNR. (A–B) Photopic ERGs recorded from normal subjects and evoked to stimuli of 0.64 or −2.23 log cd.s.m⁻², respectively. The blue ERG (A1 and B1) shows the raw waveform and the superimposed red response (A1 and B1) represents the reconstructed ERG (i.e., obtained using the inverse DWT such as in Figure 1D). Finally, the black tracings shown below the ERG responses (A2 and B2) represent the baseline noise and residuals (i.e., obtained by subtracting the reconstructed ERG from the raw ERG). The noise amplitude (see vertical calibration bars) is about the same for both ERGs. In (A), the noise has a negligible effect (Raw ERG ≈ Reconstructed ERG) on the ERG measurements (since the SNR is high). However, as shown in (B), the noise has a significant corrupting impact on the accuracy of the measurements (since the SNR is low). Thus, at low SNR, the noise contributes to the variance of the measurements. (C) DWT (dashed lines: a-wave [20a + 40a] and b-wave [20b + 40b] energy) and TD (full lines; a- and b-wave amplitude) variability (expressed as the CV) of a-wave (blue curves) and b-wave (red curves) measured in ERG responses evoked at two different stimulus intensities (−2.23 and 0.64 log cd.s.m⁻²).

DWT analysis of ERGs from patients affected with ON or OFF pathway anomaly

Results presented above demonstrate that selected DWT descriptors follow significantly different LR functions when evoked to progressively brighter stimuli, such as the 40b energy (Figure 3B) being significantly higher than that of the 20b for the rising part of the PH (i.e., dimmest stimuli), and conversely, the 20b being significantly higher than the 40b for the descending part of the PH (brightest stimuli), thus suggesting that they might quantify distinct physiological processes (or pathways) of the retina. To further investigate this claim, ERGs of patients affected with known retinal pathway anomalies, specifically affecting the ON or OFF pathway, were analyzed using the DWT. 

This is best exemplified in Figure 7A, which shows representative ERG tracings obtained from patients affected with Type-1 CSNB (specific ON-pathway anomaly; Tracings 2, 3, and 4) or CPCPA (specific OFF-pathway anomaly; tracings 5, 6 and 7). As can be seen, the ERG waveform morphology of these patients is strikingly different from that of a normal control (Tracing 1) and between them (compare Tracings 2, 3, and 4 with Tracings 5, 6, and 7). As a result, several of their DWT descriptors differed significantly. This is best evidenced with the accompanying group data and statistics reported in Table 2. As shown, the Hölder exponent was the most affected DWT descriptor in both groups and was found to be as low as 6.6 and 8 SDs below the control value for CSNB and CPCPA, respectively. Of note, out of the 10 descriptors presented in Table 2, it is also the Hölder exponent that had the least variability in control subjects (CV of 4.14%), which could explain its superior sensitivity to pathological changes. Moreover, from a TD point of view, the amplitude of the a- and b-waves (TD a and TD b parameters in Table 2) of both patient groups was significantly reduced (apart from the normal a-wave amplitude found in CPCPA patients). However, these parameters were only the fifth and sixth most affected descriptors (for the TD b parameter of CPCPA and CSNB, respectively), suggesting that these traditional ERG parameters are less sensitive than DWT ones (maximum of 1.2 and 2.5 SDs below the mean for TD a and TD b, respectively, compared to a maximum of 8 SDs below the mean for the Hölder exponent). 

Figure 7

 

Analysis of ERGs from patients affected with ON- and OFF-cone pathway anomalies. (A) A representative normal ERG response (Tracing 1) evoked to a stimulus of 0.64 log cd.s.m⁻² is compared to representative pathological ERGs (also evoked at 0.64 log cd.s.m⁻²) obtained from patients affected with an ON-specific (CSNB: Tracings 2 [composite response obtained from the average of all patients], 3, and 4 [representative responses obtained in two different CSNB patients]), or an OFF-specific (CPCPA: Tracings 5 [composite response obtained from the average of all CPCPA patients], 6, and 7 [representative responses obtained in two different patients]) cone pathway anomaly. The a- and b-waves are indicated as “a” and “b.” The black arrow (at top of Tracing 1) indicates stimulus onset. The horizontal and vertical calibration bars apply to each ERG. Associated DWT scalograms in which the LWM descriptors have been surrounded by white borders are shown on the right-hand side of the waveforms. For each scalogram, the color-coding was normalized to the LWM of maximal energy. The frequency bands and the color bar are indicated on the top and bottom scalogram, respectively. White arrows indicate a preferred reduction of the 20b (scalograms of CSNB patients) or of the 40b (scalograms of CPCPA patients; note the delayed positon of the 40b in all cases) descriptors. (B–C) Group data showing the values (M and SD) of the 20b (B) and 40b (C) descriptors obtained in control subjects (black bars) and in CSNB (blue bars) and CPCPA (red bars) patients. (D) Values of the 40b-to-20b ratio obtained in control, CSNB, and CPCPA. Dashed line indicated a unitary ratio (i.e., 40b-to-20b = 1). Significant differences are indicated on each bar graph.

Table 2

 

Analysis of ERGs obtained from congenital stationary night blindness (CSNB), congenital postreceptoral cone pathway anomaly (CPCPA), and control groups. Notes: Column 1: Names of the descriptors and parameters. Column 2: Values (M ± SD) obtained in CPCPA group. Column 3: Values (M ± SD) obtained in CSNB group. Column 4: Percent difference (computed as [(maximal value – minimal value) / minimal value] × 100) between the two patient groups. Column 5: Values (M ± SD) obtained in control group. Bold: Value that is significantly different from control; Italic: Significant differences between CPCPA and CSNB; Bold + Italic: significant differences between CPCPA, CSNB, and Control. Values between square brackets (columns 2 and 3) indicate the number of standard deviations below or above the control mean (5th column) for that parameter, thus indicating the descriptions that were the most/least affected compared to control (highest value = most affected; lowest value = least affected). Values between parentheses (columns 5) indicate the coefficient of variation.

 

Significant differences were also found between the two patient groups, the most important ones being that of the 20b and 40b descriptors, which had a percent difference of 129.42% and 101.39% between the two groups, respectively. These DWT differences are better illustrated in the bar graphs presented in Figure 7B and C. As shown, the 20b (Figure 7B) and 40b (Figure 7C) descriptors were more specifically reduced (by 2.8 folds and 3.0 folds; p < 0.00001) in CSNB and CPCPA, respectively. As a result, while in control the 40b-to-20b ratio (i.e., 40b divided by 20b; shown in Figure 7D) was found to be almost unity (1.05 ± 0.06), that of CSNB patients was found to be of 2.01 ± 0.30 (i.e., 40b > 20b), and that of CPCPA of 0.43 ± 0.06 (i.e., 20b > 40b). Use of this ratio significantly reduced intersubject variability (CV of control is 5.71% for 40b-to-20b ratio compared to 18.18% and 14.92% for 20b and 40b, respectively) and significantly emphasized the effect size seen between the two patient groups (% difference = 367.44% for 40b-to-20b ratio compared to 129.42% and 101.39% for 20b and 40b, respectively). 

It is clear from the above that DWT descriptors can be specifically or differently affected by a given disease process and, consequently, that they could be used to highlight meaningful differences between various ERG responses. Clearly, the data shown in Figure 7 and Table 2 provide interesting demonstrations of how DWT analysis could complement the traditional analysis of the ERG by offering an alternative quantification, which could better reflect the (presumably) unlimited ways by which retinal function (as reflected with the ERG) can be modulated as a consequence of normal (Figures 3 and 4) and pathological processes (Figure 7). 

Discussion

Reproducible LWM, time-locked with the a-wave, b-wave, and OPs were identified in high- (>100 μV; Figures 3 and 4) and low- (<1 μV; Figure 5) voltage ERGs, a feature of diagnostic relevance especially if one wishes to use the ERG to monitor disease progression in severe degenerative retinopathies (such as Retinitis Pigmentosa), whose final outcome is often characterized by very low-amplitude ERGs, or even extinguished ERGs (Rispoli, Iannaccone, & Vingolo, 1994). Our results indicate that the LR functions that were generated using the DWT descriptors that quantified ERG components concomitant with the b-wave (i.e., 20b and 40b) presented with PH-like shapes that complemented the traditional TD measurement of the b-wave (TD b-wave curve in Figure 3E). A similar match was also found between the DWT and TD descriptors of the a-wave. The fact that we were able to mimic the characteristic PH, a signature trait of the cone ERG LR function (traditionally generated from TD measures), using DWT descriptors suggests that selected DWT descriptors appraised physiological attributes of the cone ERG that covaried (see Pearson coefficients in Figure 3) with those obtained with TD measures. 

As illustrated in Figure 8, there were also some important differences between the various DWT LR functions, such as the intensity of stimulation at which the maximal values were reached. DWT LR functions disclosed two maximal peaks, one at 0.39 log cd.s.m⁻² (maximal values for 40b, 80ops, 160ops, and ΣEnergy; Peak 1 in Figure 8) and another one at 0.64 log cd.s.m⁻² (maximal values for 20b, Hölder exponent and Δ-variance; Peak 2 in Figure 8). In a previous study, we also showed that the flash intensity required to reach the maximal b-wave amplitude and to optimally develop the photopic OP response (i.e., where OP2, OP3, and OP4 are fully developed) was (approximately) 0.3 and 0.6 log cd.s.m⁻², respectively (Lachapelle et al., 2001). The brighter intensity also generated the most complete response (in terms of number of detectable ERG components), a claim in accord with our finding that the Hölder exponent (descriptor of roughness or complexity of a waveform) also peaks at the brighter flash intensity (0.64 log cd.s.m⁻²). Moreover, as it can be seen in Figure 8, while some DWT descriptors (40b, 80ops, 160ops, and ΣEnergy) reached a plateau with the brightest stimuli (from 1.63 to 2.23 log cd.s.m⁻²), others either continued to decrease (Hölder and Δ-variance) or slightly increased (20b). 

Figure 8

 

Normalized (between 0% and 100%) LR functions derived from the DWT descriptors that presented with various hill-like shapes in the present study (see figure legend). The black curve (identified as TD b) shows the traditional LR function of the b-wave amplitude measured in the TD. The solid vertical lines (lines 1 and 2) indicate the 0.39 and 0.64 log cd.s.m⁻² stimulus intensities that evoked the maximal values (Peak 1 and Peak 2) of the descriptors. As seen in the insert, four DWT descriptors reached their maximal value at Peak 1 and three descriptors at Peak 2.

DWT analysis of the ERG also allowed us to further dissect the a- and b-waves into separate subcomponents, each oscillating at a specific frequency band centered at 20 and 40 Hz which, from the results shown in Figure 3A and B, appears to be distinct from each other. The latter indicates that the intensity of the stimulus significantly modulates the time-frequency composition of the ERG and consequently its morphology as evidenced with the different ERG waveforms shown in Figure 2A. Knowledge of the latter not only provides us with a more refined approach to ERG shape quantification, but also suggests that some retinopathies could, for example, preferentially affect one frequency band (or DWT descriptor) more than the others, thus permitting a more precise ERG-based segregation of different retinal disease processes. The latter was demonstrated in Figure 7 (and accompanying Table 2) where patients affected with ON- and OFF-pathway anomaly had a specific attenuation of the 20b and 40b descriptors, respectively. Analogous findings were also reported by Barraco et al. (2011), who used the CWT to compare the scotopic a-wave of normal and achromate subjects. They showed that the scotopic a-wave of normal subjects was composed of three frequency components (20, 140, and 180 Hz) and that the higher frequency component was severely reduced in achromates, further suggesting that some retinopathies can preferentially affect one frequency band more than the others (Barraco et al., 2011). 

The DWT also offered an alternative approach to assess the OPs without having to filter the broadband ERG signal. As shown with the scalograms of Figures 2B and 5A, OPs descriptors were computed as the average of five coefficients and, as a result, did include some low-energy (blue) coefficients, which contributed no more than those found outside the selected (white) boundaries. Of interest, despite this more inclusive approach, we nonetheless obtained very high correlations (r = 0.9831) between TD and DWT OPs measurements (see Figure 3F). Obviously, we might have obtained better fits by selecting which coefficients to include on a case-by-case basis (cherry-picking approach), a strategy that would have required some data interpretation from the user or more complex interventions from the algorithm, all of which could potentially lead to undesired artifacts, subjectivity and variability. Although our OPs descriptors (80ops and 160ops) do not, unlike band-pass filtering, permit the analysis of the photopic OPs individually, our results (see Figure 3C) did show that the two OP frequency bands differed in their respective energy level, suggesting that they might monitor different features of the high-frequency components of the ERG. Of interest, another wavelet study (Forte et al., 2008) previously revealed that rat OPs also comprised two frequency components (i.e., 70–80 Hz and 120–130 Hz). This study also showed that the two OP components had distinct LR functions, suggesting that they might be evoked by different retinal sources. Similarly, Zhou, Rangaswamy, Ktonas, and Frishman (2007) reported that the photopic OPs of primates also contained two distinct frequency bands, namely the slow OPs (presumably generated by the amacrine cells) oscillating at 80 Hz and the faster OPs (presumably generated by the ganglion cells) oscillating at 150 Hz (Zhou et al., 2007). Like our study, the slower OPs were of the highest energy level. In our study, the LR functions of the 80ops and 160ops DWT descriptors followed a PH-like function, suggesting an intimate tie linking the genesis of the b-wave with that of the OPs, a claim previously advanced using TD measures (Guite & Lachapelle, 1990; Lachapelle, 1987, 1990; Lachapelle & Benoit, 1994; Lachapelle & Molotchnikoff, 1986). 

Although both OPs descriptors (80ops and 160ops) were found to be significantly reduced (compared to control) in both CSNB and CPCPA patients (see Table 2), we did not find any significant difference between the 80ops or 160ops descriptors of the two patient groups. This contrasts with a study that was previously published by us in which we showed that these diseases were separable based on OPs measurements as patients affected with CSNB more specifically lose early OPs (OP2 and OP3), while patients affected with CPCPA preferentially lose the later OP (OP4; Lachapelle et al., 1998). However, in this previous study (Lachapelle et al., 1998) we looked at the amplitude of OPs individually, while herein, we measured the overall OPs energy (average of five coefficients; see Table 1). This is not to say that it is impossible to look at early or late OPs individually with the DWT. Indeed, there are almost an infinite number of descriptors that one can calculate from the DWT scalogram, and individual assessments of the early and late coefficients of the 80ops or 160ops descriptors (instead of averaging the five coefficients) could have been used to quantify the energy level of early and late OPs. Alternatively, the ratio between the 160ops and 80ops descriptors could have been used to detect differences between CSNB and CPCPA. Supportive of the latter, in CSNB and CPCPA patients, the value of the 160ops-to-80ops ratio (obtained by dividing 160ops by 80ops) is of 0.49 ± 0.09 and 0.72 ± 0.11, respectively—values that are significantly different (p < 0.00001) from each other. This result indicate that the low- (80 Hz) and high-frequency (160 Hz) OP bands can be differently affected by a given disease process and suggest a potential association of the low- (80 Hz; i.e., 80ops descriptor) and high-frequency (160 Hz; i.e., 160ops descriptor) OP descriptors to the OFF- and ON-pathway, respectively. In a recent study, Dimopoulos et al. (2014) reported that the light-adapted OPs were separated into two frequency bands and suggested the possibility of ON and OFF system representation. Data presented herein would support this claim. 

The Hölder exponent (a measure of waveform complexity) also varied significantly as a function of the intensity of the stimulus (Figure 4). Although signal complexity remains somewhat of an abstract concept, the Hölder exponent can be defined as an ERG descriptor that accounts for the number of distinct elementary components (or frequency bands) that compose the signal, as well as the local irregularities that these components add to the waveform when they combine together. The notion of complexity is thus related to signal morphology because signals that do not have the same components will necessarily be of distinct shapes. For instance, in the examples shown in Figure 1C, the Hölder exponent increases from 0.53 (20 Hz ERG) to a maximum of 1.2 (for the 20 + 40 + 80 + 160 Hz ERG). As can be seen in Figure 8, in the initial portion of the PH, the Hölder exponent (green curve) increases at a faster rate than the TD and DWT parameters, while all three parameters appear to be more similarly modulated in response to the brightest intensities. This was also seen with the lower amplitude ERGs (compare Figure 5D, E). However, since the Hölder exponent is amplitude-independent, we did not expect that it would follow the traditional PH shape. That is, of course, unless we postulate that the ERG waveform would become progressively more complex as its amplitude increases. Consequently, we believe that the faster increase of the Hölder exponent observed in the initial rising phase of the PH simply confirms previous reports that claimed that the photopic ERG signal was a composite potential where brighter stimuli added new components (i.e., thus increasing waveform complexity) to the threshold ERG waveform (Berson, Gouras, & Hoff, 1969; Lachapelle & Molotchnikoff, 1986). Of all the ERG descriptors (TD or DWT) considered in this study, the Hölder exponent was that which was shown to be the least variable (CV = 4.14% [Hölder] compared to 17.25% [b-wave amplitude]) for the ERG at 0.64 log cd.s.m⁻², suggesting that it could allow us not only to detect subtle pathological changes at disease onset, but also to facilitate the monitoring (and severity grading) of changing ERG morphologies as the disease progresses towards nearly extinguished responses. This claim is best supported with the gradual decline in the value of the Hölder exponent seen with progressively smaller ERGs in normal subjects (Figure 5E) as well as with the significantly attenuated values obtained from pathological responses (Figure 7; Table 2). 

Of interest, when analyzed with the DWT descriptors, both patient groups (CSNB and CPCPA) presented in Figure 7 and Table 2 had a unique DWT scalogram signature and were therefore statistically separable from each other using five of the eight DWT descriptors (values in italic in Table 2), but were not using TD descriptors. The latter indicates that DWT descriptors of the ERG can complement the diagnostic information obtained using TD descriptors. The DWT descriptors were also more affected (up to 8 SDs below the mean) than the traditional TD parameters (maximum of 2.5 SDs below the mean). The latter would suggest that DWT descriptors could also potentially detect subtle ERG anomalies even when the TD measures (amplitudes and timing) are still normal. Previous studies conducted on more specialized types of ERG signals further support the latter claim. For example, Miguel-Jimenez et al. (2010) used the DWT decomposition to reconstruct mfERG waveforms into various frequency bands and demonstrated the higher sensitivity of this approach to detect changes in glaucoma patients, compared to classical Humphrey visual field tests. Similarly, Rogola and Brykalski (2005) previously showed that DWT coefficients were superior to traditional TD measures in segregating normal from pathological PERG waveforms. Results presented in Figure 7 and Table 2 also indicate that certain disease can affect specific DWT descriptors. As shown, the photopic ERGs of CSNB patients suffered from a pronounced attenuation of the low-frequency b-wave component (20b) compared to the better preserved high-frequency b-wave component (40b), and conversely, CPCPA patients were characterized by opposite findings (40b specifically reduced compare to the better preserved 20b). Of interest, given that patients affected with CSNB and CPCPA have a specific anomaly of the ON and OFF retinal pathway (Bech-Hansen et al., 2000; Dryja et al., 2005; Garon et al., 2014; Gregg et al., 2007; Lachapelle et al., 1998; Langrova et al., 2002; Miyake et al., 1987; Pusch et al., 2000; Quigley et al., 1996) and that we found a specific attenuation of the 20b and 40b in CSNB and CPCPA, respectively, our results suggest that the 20b and 40b descriptors might therefore be more closely related to the ON (20b) and OFF (40b) cone pathway contribution to the ERG response, respectively. 

Previous studies of the PH phenomenon—using a variety of approaches such as early/late OPs analysis (Rufiange et al., 2002), long-flash ERG (Kondo et al., 2000), mathematic models (Garon et al., 2014; Hamilton et al., 2007), etc.—suggest that the initial rising phase of the PH is OFF-response–dominated while the descending phase/plateau is ON-response–dominated. The fact that we have found (see Figure 3) significantly higher values of the 40b (OFF) compared to 20b (ON) in the rising portion of the PH, and conversely, significantly higher values of the 20b (compared to 40b) in the descending portion and plateau of the PH, further supports that the 20b and 40b descriptors can be associated to the ON and OFF response, respectively. Furthermore, in another study conducted on monkeys, Ueno et al. (2004) used pharmacological blockages to isolate the ON and OFF components of the primate cone ERG and showed that the OFF-component was delayed with higher intensities. Of interest, we also found that the 40b component occurred at a delayed position with the six brightest intensities (see white arrows in Figure 2C). Based on the above, our findings indicated that the PH phenomenon of the human ERG effectively results from the combination of both ON- and OFF-pathway activity—the initial rising phase being characterized by a greater contribution of OFF-pathway (40b) energy and the descent phase and plateau being characterized by a greater contribution of ON-pathway (20b) energy. Based on the above, it is not surprising that the ERG response of maximal complexity (and of high amplitude; b-wave amplitude >100 μV) was reached at the intensity of 0.64 log cd.s.m⁻², where both the ON (20b) and OFF (40b) pathway are equally (and maximally) contributing to the response. 

Finally, although the DWT analysis of the photopic ERG did permit the detection of major differences between CPCPA and CSNB waveforms, it does not mean that the DWT is the only technique that is able to segregate the two signals. As aforementioned, OP measurements in the TD did allow us to distinguish the two retinal disorders (Lachapelle et al., 1998). Likewise, Fourier analysis (such as the FFT) could also have been used to distinguish the two abnormal ERG signals by determining the frequency of the fundamental component (see Supplementary Figure S1). However, one should keep in mind that the FFT suffers from significant drawbacks when compared to the DWT. For example, the absence of temporal resolution in the FFT prevents the a-wave and b-wave from being independently quantified (since the a- and b-waves oscillate at overlapping frequencies) and also precludes the detection of latency shifts of the ERG components (Gauvin et al., 2014). Of interest, in the current article, the DWT allowed us to detect significant differences between the a-waves of CSNB and CPCPA patients and detected delays in the 40b component of photopic ERGs evoked from CPCPA patients as well as in ERGs recorded from normal subjects in response to the brightest stimulus intensities. Indeed, both of these potential diagnostic features are inaccessible if frequency domain analysis of the ERG waveform is limited to the FFT. 

Conclusions

This study represents the first attempt to study the luminance-dependence of ERG components extracted with the DWT and the first application of wavelet analysis to study the ERG recorded from patients affected with retinal diseases impairing the retinal ON and OFF pathways. We characterized the LR functions of novel DWT descriptors of the human photopic ERG and demonstrated the usefulness of these descriptors in quantifying normal and pathological ERGs signals beyond what can be accomplished using the classical TD measures (e.g., a- and b-wave amplitude). Our results also suggest the possibility that a given retinopathy may more specifically impair one of the DWT descriptors, raising the possibility that alterations of one or more DWT descriptors (or combinations of) could be pathognomonic for a given disease process (as suggested from the data shown in Figure 7 and Table 2). The DWT thus offers an alternative approach to identify (early diagnosis), classify, and possibly stage the different pathophysiological processes that impair retinal function—information that will be highly relevant, especially when assessing diagnostically challenging cases where the retinal disorder has not yet impaired TD measures of the ERG. Of interest, while delays in the ERG response (i.e., a TD measure) remain a critical descriptor of pathological ERGs, as shown in Figure 7, significantly delayed ERG b-waves can also retard the temporal position of some b-wave descriptors (such as that of the 40b descriptor in the CPCPA ERGs shown in Figure 7A), demonstrating that latency shifts can also be detected (and thus measured) with the DWT scalogram. Finally, although our data strongly support the fact that ERG analysis is more complete in the time-frequency domain, this is not to say that analysis of the ERG in the TD should be abandoned. Rather we advocate the use of the DWT approach as a complement to TD measures of the ERG. 

Acknowledgments

This work was supported by grants-in-aid from the Canadian Institutes of Health Research (CIHR MOP-126082), the Fonds de recherche du Québec–Santé (FRQ-S) and its thematic research network (Vision Network) as well as by the CIHR ERA-132932 and FRQ-S JTC 2013, under the frame of E-Rare-2, the ERA-Net for Research on Rare Diseases. Thanks are due to Marie-Lou Garon, Marianne Rufiange, Julie Racine, and Allison L. Dorfman, who performed the ERG recordings. 

Commercial relationships: None. 

Corresponding author: Pierre Lachapelle. 

Email: [e-mail]. 

Address: Research Institute of the McGill University Health Centre, Montreal, Quebec, Canada. 

Appendix A: Methodological details

Calculation of the DWT

In order to localize the energy content of the ERG in both time and frequency we computed the DWT of each ERG as follows:  

 

DWT(j, k) represents the wavelet coefficients localized at discrete scales (indexed with j and corresponding frequencies in [(Fs / 2)2^−j, (Fs / 2)2^−j+1], where Fs is the sampling frequency) and discrete time k, x(t) designates the raw ERG time series, and ψ̄ denotes the complex conjugate of the mother wavelet. The DWTs were computed using the fast wavelet transform algorithm of Mallat (2009) implemented with MATLAB and Wavelab 850 routines (Buckheit, Shaobing, Donoho, Johnstone, & Scargle, 2005). Prior to DWT computation, each ERG was padded with 256 constant samples (by repeating the first and last value of the signal) on both sides of the response, in order to reduce edge effects (Torrence & Compo, 1998). With our parameters (i.e., 1,024 samples per padded signal and Fs of 3413.33 Hz), we were able to obtain 10 levels of decomposition per DWT. Only the first eight levels were displayed and the padding was discarded after computation (obtaining an eight-level time-frequency plan of 150 ms in length). Finally, several mother wavelets could have been used to extract the DWT descriptors described herein. However, one should be aware that use of different wavelets may affect the output of the DWT (such as the energy level), and hence, the same mother wavelet should be applied to normal subjects and patients for clinical decisions. In this study, we opted for the Haar wavelet to analyze the LWM descriptors of the ERG and for the symmetric Daubechies wavelet (set with two vanishing moments) to compute the WVA descriptors. Rationale for the use of these mother wavelets is indicated in the Methods section (see subsection “Selection of mother wavelets”).  

DWT Quantification of the ERG waveform with LWM descriptors

Six LWM (illustrated in Figure 1A and defined in Table 1) were considered in this study, namely 20a, 20b, 40a, 40b, 80ops, and 160ops. Furthermore, given that translation of the ERG response is expected (such as peak-time variation of ERGs) and that the DWT is not a shift-invariant transform (Guo, 1995), we computed, as previously suggested (Coifman & Donoho, 1995), the DWT of several translated versions (translation of the entire ERG waveform) of the same ERG signal in order to optimize the identification of the maximal LWM independently of the timing (peak time) of the a-wave, b-wave, and OPs. These translations allow the various components of the ERG to be optimally aligned under their associated wavelet boxes (such as that of the a-wave under the 20a and 40a boxes), thus compensating for the lack of shift-invariance of the DWT. Each LWM was thus calculated as the local maximum value obtained while shifting the ERG waveform to a given extent. The width of these translations is reported in Table 1 (column 5) for each LWM descriptor. Each translation was conducted with increments of 1 ms, and negative and positive translation values indicate displacement of the ERG signal to the right and left, respectively. As indicated in Table 1, the translation extent of the a-wave descriptors (20a and 40a) does not allow as much translation to the left (i.e., positive translation values) compared to that allowed for the b-wave descriptors (20b and 40b). In limiting leftward translations we prevent the alignment of the b-wave under wavelet boxes associated with the a-wave (20a and 40a) and therefore avoid the potential contamination of the b-wave on the 20a and 40a descriptors (see Supplementary Figure S2). One also notes (see Table 1) that the higher the frequency, the smaller is the magnitude of the translations. This is due to the fact that the higher the frequency is, the smaller the scale will be (i.e., width of the rectangular boxes in the scalogram), resulting in a higher temporal resolution of the DWT. Consequently, for smaller scales, minimal translations are needed to obtain optimal alignments. When the LWM descriptors are directly assessed without translation of the ERG responses, the resulting energy is considerably underestimated. This is better illustrated in Supplementary Figure S3, comparing the DWT LR function of the a- and b-waves obtained with (red curves) and without (blue curves) translations. Of note, although LR functions assessed without translations lead to an underestimation (i.e., inaccurate measure) of the energy, they are nonetheless reminiscent of those obtained with translations, in that they exhibit similar stimulus-dependent patterns. 

DWT quantification of the ERG waveform with WVA

Wavelet variances (expressed as standard deviation) were computed using the values of 256 wavelet coefficients for the first DWT level (i.e., 1280-Hz band) down to two coefficients on the last level (Level 8; 10-Hz band). Note that in the DWT scalogram, the wavelet coefficients are seen as the colored rectangles of various sizes (see Figure 1A; the 256 smallest rectangles are found on the first level and the two largest on the eighth level). The eight standard deviation values thus obtained were then plotted against their respective DWT level (see example in Figure 1D). This graphic representation allowed us to define the Δ-variance descriptor, which, as illustrated in Figure 1D, was defined as the difference between the variance measured at Level 6 and its prediction at the same level based on the linear fit between Levels 2 to 5. The linear fitting was achieved between Levels 2 and 5 as this portion of the curve permitted an optimal linear fit, while Level 6 is where the curve stopped behaving linearly (see major discontinuity between Levels 5 and 6 in Figure 1D). The Hölder exponent is an amplitude-independent descriptor used to assess the regularity (often associated to the complexity) of a signal. The DWT offers an optimal scheme to estimate this exponent. This is due to the fact that the variance of the wavelet coefficients d at a given level j varies as Var(d_j) ∼ 2^jH. Accordingly, the Hölder exponent H (also termed scaling exponent) can be estimated from the DWT variance-level plot, such as that illustrated in Figure 1D (Abry et al., 2002). Briefly, we process the natural logarithm of the variance-level curve (Figure 1D) to obtain the log-variance plot, such as the one shown in Figure 1E. This logarithm linearizes the curve (compare Figure 1D, E). The measurement of the Hölder exponent is then reduced to the calculation of the slope over the alignment region in the log-variance diagram. In this study, the linear fitting was achieved between Levels 2 to 6 (red line in Figure 1E) as this region of the curve offered the best linear fitting.