1. Introduction

Signals from the three cone spectral classes in the retina form the basis for color perception [1]. Their relative proportion and arrangement shapes and constrains how wavelength information from the external world is coded in the visual system to create the rich palette of hues that most humans perceive. Short-wavelength cones (S-cones) can be differentiated in histology [2], but the remarkably similar morphology and protein structure of long and middle-wavelength cones (LM-cones) has long precluded the creation of imaging techniques and antibodies to segregate them in a cellular-scale assay ex vivo. Until the advent of adaptive optics (AO) to overcome the optical aberrations of the eye [3] and resolve cone photoreceptors in living humans, there was virtually no information about the arrangement of individual LM-cones in primate retina. In fact, AO retinal densitometry to delineate the spectral types of cones by Roorda and Williams [4] was one of the very first scientific applications of AO upon its introduction to the eye in 1997.

Densitometry has long been used to study various biophysical properties in cone photoreceptors, including waveguiding efficiency of the inner and outer segments, optical density, and photosensitivity of chromophores housed in their outer segments [5–8]. Its principle rests on using the absorption of light by photopigments as a gauge to measure the aforementioned biophysical properties. For cone spectral typing, differential absorption by the three photopigments in LMS-cones is used to separate them by their spectral type. An AO fundus camera was first used for this task, enabling the visualization of cone spectral types in humans and non-human primates [4,9,10]. In 2015, the feasibility of using an AO scanning laser ophthalmoscope (AOSLO) to classify cone spectral types was tested [11]. In addition to the advantages of increased axial and lateral resolution, contrast and signal-to-noise ratio, AOSLO densitometry allowed imaging the dynamic changes in photoreceptors upon light capture. Together, this enabled cone spectral classification with high efficiency with the time taken per subject ranging between 3 to 9 hours, and requiring an average of 8–15 repeat measurements per bleaching condition to boost the signal-to-noise ratio. In contrast, the AO fundus camera needed averages of up to 50 repeat measurements across 5 days in each bleach condition [4].

Recently, Zhang et al. 2019, 2021 [12,13] showed improved precision and efficiency for mapping cone spectral types using an AO point-scan spectral-domain optical coherence tomography (OCT) instrument in phase-resolved configuration. The principle of the measurement rests on imaging light-evoked changes in the backscattered optical phase arising in the cone outer segments in response to a bleaching stimulus, also referred to as optoretinography or ORG [14,15]. These changes in phase were converted to optical path length (ΔOPL) to serve as a quantifiable gauge for measuring light propagation after bleaching in the outer segment and the magnitude of ΔOPL scaled with increasing bleach strength. By adopting a bleaching light that differentially activated L, M and S-cones, it was demonstrated that the ΔOPL response arising from individual cones could be readily segregated into three cone sub-types. In 2020, Pandiyan et al. demonstrated a line-scan AO-OCT instrument for volumetric imaging of the retina that could record retinal cross-sections up to a maximum rate of 16000 B-scans/sec [14,16]. The increased speed is the result of parallel acquisition of a linear spatial imaging field in a 2D camera after diffraction through a transmissive grating, such that after a Fourier transform, a cross-sectional B-scan is obtained in a single acquisition. This comes at the cost of loss of confocality in one dimension compared to a traditional point-scan configuration. Despite this, it was demonstrated that the instrument could image different scattering structures within the retina at cellular resolution, including retinal ganglion cells, macrophages, rod cells, foveal cones and retinal pigment epithelial cells [17]. Using this line-scan OCT instrument, Pandiyan. et. al [16,17] demonstrated the feasibility of cone spectral classification, following the general protocol established by Zhang et al. 2019. Notably, cone spectral types were classified at an eccentricity of ∼0.3° from the foveola, enabling the first in vivo demonstration of decreasing S-cone proportion in the fovea [17].

It is instructive to compare the ORG and densitometry in terms of their mechanism of action and dynamic range. As explained below, ORG has 20x greater dynamic range for assessing photopigments across cone spectral types as compared to densitometry. This provides a fundamental benefit to the signal-to-noise ratio in ORG for classification; signal corresponds to a surrogate measure of photopigment concentration while noise is contributed by various common factors for both including low retinal reflectivity, the limit to light collection set by the numerical aperture of the eye, eye movements, reflections from retinal layers besides the cones, blood flow and more. The higher dynamic range is expected to reach similar or better accuracy in ORG with lower numbers of repeat measures, thus improving speed and efficiency.

Densitometry uses the same light source for imaging photoreceptors and bleaching photopigment, typically centered at 550 nm where the spectral sensitivity of LM-cones is high, and similar to each other. After dark-adaptation and regeneration of cone photopigment, the retina is exposed to the imaging and bleaching light. The backscattered light encodes the cone locations in the retinal image, and simultaneously carries information about absorption by the photopigment. For instance, the S-cones are relatively insensitive to a 550 nm imaging light and therefore maintain their backscattered intensity after exposure to it. LM-cone photopigments on the other hand absorb this light and appear darker than the S-cones immediately after light onset. With increasing bleach of LM-cone photopigments, their concentration, and hence absorption decreases, consequently leading to an increase in the backscattered image intensity. This relationship between image intensity and photopigment concentration for LM-cones is depicted in Fig. 1(a) (brown and black dashed curves respectively). The relative change in image intensity immediately after stimulus onset and after a full bleach is referred to as the optical density and signifies the maximum photopigment-dependent signal available for probing its characteristics. Optical density is computed as the base 10 logarithm of the ratio of image intensity, post to pre-bleach. Neitz et al. 2020 [18] measured individual cone optical density using AOSLO densitometry for a range of visible wavelengths from 496 nm to 598 nm. A maximum optical density of 0.41 log10 units (2.6 times) was observed at 543 nm, i.e., a 2.6x change in image intensity is available as the dynamic range to probe LMS-cone photopigments.

 

Fig. 1. The dynamic range of ORG is ∼20x improved over densitometry. Densitometry (Fig. 1(a)) measures the time varying change in the backscattered image intensity (brown curve, left y-axis, Fig. 1(a) from the retina as the photopigment concentration (black dashed curve, right y-axis, Fig. 1(a) is reduced in the cone outer segment with a bleaching stimulus. ORG (Fig. 1(b)) measures the time varying change in ΔOPL in the cone outer segment (green curve, left y-axis, Fig. 1(b) as the photopigment concentration (black dashed curve, right y-axis, Fig. 1(b)) is reduced with a bleaching stimulus. The dynamic range is the maximum obtainable change in the photopigment-dependent signal attributed to a bleach (represented on the left y-axis for both), equal to ∼2–2.6x for densitometry and up to at least 50x for ORG.

Download Full Size | PDF

In contrast, the stimulus-evoked ΔOPL shows a significantly greater pigment-dependent change, scaling in magnitude from 10–500 nm with increasing bleach from 0.3 - 50% [14]. Figure 1(b) shows ΔOPL vs. time for a 50% bleach using a 528 nm stimulus. Considering 10 nm as the noise floor of the measurement, the ratio of post to pre-bleach ΔOPL can be as high as 50x, indicating the dynamic range available to probe LMS-cone photopigments in ORG. The mechanism of action of ORG is has recently been proposed and its bleach dependence suggests an involvement of the phototransduction amplification cascade [14,19, 24] that ultimately benefits its dynamic range.

In this article, we sought to investigate the correspondence between densitometry and ORG for cone classification in the same subjects and compare their fidelity for the task. First, we establish our ORG protocol for this task and ascertain its test-retest reliability by measuring the same subjects and eccentricities across different days. Variations in reflectance of cones caused due to interference, disc shedding and renewal, and shadows of blood vessels negatively affect both paradigms for classification [20–22]. Next, we ascertain whether cones with abnormal reflectivity exhibit normal function in ORG sufficient to identify their type. Finally, we compare the spectral assignments between the two paradigms with respect to the specific cone types and their spatial arrangement.

2. Methods

2.1 Reflective mirror based line-scan AO-OCT

A line-scan spectral domain OCT system equipped with AO, described previously [17], was used for ORG with cellular resolution. Briefly, a wavelength band centered at 820 ± 40 nm was selected from a supercontinuum light source (SuperK EXR-9OCT, NKT Photonics, A/S, Birkerød, Denmark) for imaging the retina. A 980 ± 10 nm superluminescent diode (SLD) (IPSDD0906, Inphenix, USA) was used for wavefront sensing and AO correction. A 660 ± 5 nm light-emitting diode (LED) was used as the bleaching stimulus to differentiate the cone types because L-cones have 13 times higher sensitivity than M-cones, and S-cones have negligible sensitivity at this wavelength [23]. The wavefront sensing and the imaging beams were combined with a short pass filter and these two beams were combined with the bleaching stimulus using a long pass filter. Light backscattered from the retina was split into two beams using the same short pass filter in the detection arm. The 820 nm light scattered from retina was captured either with a line-scan camera to optimize photoreceptor image focus in an en face image stream, or with a custom-built spectrometer for OCT. The spectrometer consisted of a diffraction grating (WP-600/840-35X45, Wasatch Photonics, USA) and a fast 2D camera (FASTCAM NOVA S16, Photron, USA). The fast 2D camera was used to capture the spatial and spectral dimensions along its 2 axes simultaneously to yield a cross-sectional OCT B-scan in a single frame. An anamorphic configuration was used to optimize the spatial and spectral dimension independently [16]. The sensitivity of the system was 92 dB.

2.2 Subjects

Five color normal subjects were recruited for the study after an informed consent explaining the nature and possible consequences of the study. The research was approved by the University of Washington institutional review board and experiments were performed in accordance with the tenets of the Declaration of Helsinki. Of the 5 subjects, one was Caucasian male, two were female Caucasians and two were Asian males. Two of the color normal retinae were cone-typed earlier using densitometry, published in Sabesan et al. [11]. Subjects S3 and S4 in Sabesan et al. are S4 and S5 respectively in this paper. In this study, cone classification for S4 and S5 with ORG, at 1.5° and 2° eccentricity, respectively, was compared against data from densitometry in Sabesan et al., separated in time by 6-7 years. In all other subjects, cone classification was performed at 1.5° temporal eccentricity. For comparison of cone classification across two days, subjects S2 and S4 were imaged at 4° and 1.5° eccentricity, respectively.

2.3 Imaging protocol

Subjects were dilated and cyclopleged using 1% tropicamide and/or 2% phenylephrine. They were placed in a dental impression to stabilize head movement. Their pupil was aligned to the optical axis of the instrument using a three-axis translation stage assisted with a pupil camera. The aberrations were measured and corrected in real-time using the wavefront sensor and deformable mirror. The retinal image captured in the line-scan camera was optimized to focus the cone photoreceptors. The system was then switched to line-scan OCT mode for further experiments. For ORGs, AO-OCT volumes were recorded after subjects were dark-adapted for 1-3 minutes to allow for sufficient regeneration of cone photopigments. The 660 nm, 20 ms LED flash was delivered after the 10^th volume. We chose 660 nm because at this wavelength L-cones have greater sensitivity than M-cones by ∼13 times, and S-cones have negligible sensitivity. The estimated photon density for 660 nm stimulus was 25.1 × 10⁶ photons/µm². We prefer to report the bleach strength in units of photon density “at the retinal surface” and not in units of percent or fraction bleach per se. This is because converting the photon density to bleach strength requires accounting for factors related to the cone photoreceptor that vary between subjects and eccentricity, for example, the outer segment length, inner and outer segment diameter and Stiles-Crawford effect [24]. We optimized the stimulus duration to maximize the ΔOPL difference to separate LMS-cones. The B-scan rate was 12000 scans/sec. Each volume contained 600 B-scans, and each video contained 50 volumes. The overall recording time was 2.9 secs, wherein each volume is acquired in 58 ms (including 8 ms dead time for galvo scanner return where recording is paused).

Typically, unless otherwise stated, 10 recordings were taken for averaging and improving signal-to-noise ratio. Accounting for any rejections of acquired volumes due to motion artifacts, between 7 to 10 recordings were analyzed in each subject. The schematic for cone classification using ORG is shown in Fig. 2. The OCT volume was segmented to access the cone outer segment (Fig. 2(a)). After dark adaptation, photopigments regenerate in LMS-cones, indicated with red, green and blue shades respectively on the cone outer segment discs in Fig. 2(b). Following the bleaching stimulus, L, M, and S-cone photopigments are depleted with decreasing magnitude in turn, proportional to their sensitivity to the 660 nm wavelength (Fig. 2(c)) [23]. The timing diagram of a representative ORG cone classification experiment is shown in Fig. 2(d).

 

Fig. 2. Schematic for cone classification using ORG. (a) The registered AO line-scan OCT volume was segmented at the cone inner-outer segment junction (ISOS) and the cone outer segment tips (COST) to access the outer segment. (b) After dark adaptation, photopigments regenerate in all cones. After a 660 nm stimulus, L-cones have highest bleach and activity, followed by M-cones and then S-cones, as dictated by their decreasing spectral sensitivity to 660 nm wavelength in (c) [23]. (d) The onset of the 20 ms stimulus flash occurs at τ_d = 0.58 sec or after the 10^th volume. OCT volumes (n = 50) were recorded, reconstructed and registered. The phase difference between the ISOS and COST in each volume was calculated and treated to the steps in Section 2.4 to obtain a measure of ΔOPL in each cone, evoked by the stimulus.

Download Full Size | PDF

2.4 OCT image processing

Typical steps were followed for OCT image processing. The details of the processing and ORG extraction are detailed in Ref. [16]. Each recorded 2D spectrum $({I({x,y,\lambda } )} )$ was treated to background subtraction, k-space resampling $({I({x,y,k} )} )$, and a Fourier transform $I({x,y,z} )$ to yield an OCT volume. All volumes were registered using segmentation-based 3D registration [25]. Once the volumes were registered, each OCT volume was referenced to the mean of all the volumes that were recorded before the start of the stimulus to cancel the arbitrary phase at each pixel. Then, the mean of 3 axial pixel complex values, centered at the boundaries of the inner-outer segment junction (ISOS), ${I_{ISOS}}({x,y} )$ and the cone outer segment tips (COST), ${I_{COST}}({x,y} )$ was calculated. The phase difference ${I_{ISOS/COST}}({x,y} )$ between these two layers was obtained by multiplying the complex conjugate of COST with the ISOS layer. The complex numbers ${I_{ISOS/COST}}({x,y} )$ were averaged over the collection aperture (5 × 5 pixels, ∼3.7 × 3.7 µm) of a cone photoreceptor to yield one complex number per cone in a volume. The same was repeated for the 50 volumes in the 17.2 Hz time series to give ${I_{ISOS/COST}}(t )$.

2.5 Complex average vs. angle average

We evaluated two approaches for averaging across repeat measurements to improve signal-to-noise. The first approach involved taking the complex average of individual ΔOPL vs. time series as below:

In the second approach, the argument of each individual complex valued time series was obtained to yield the phase, unwrapped if needed and converted to ΔOPL using Eq. (2). This approach represents the angle average. The mean $\overline {\varDelta OPL(t )} $ of all recordings was calculated as

2.6 Clustering analysis

The ΔOPL vs. time trace for each cone was obtained by averaging the individual repeat trials. For segregating L, M and S cones, the ΔOPL averaged over the time after stimulus onset between 0.7–1 sec was calculated for every cone that was visible in AO-OCT structural images. A histogram of the time averaged ΔOPL (maximum ΔOPL elsewhere in this article) was fit with a sum of three 1-dimensional Gaussians. The clustering method followed here is analogous to a 3-component, 1-dimensional Gaussian mixture model [27]. The intersection of the component Gaussians separated the clusters into the three cone types. Each classified cone within a cluster is associated with an error percent or uncertainty, and a probability or likelihood of its assignment. Error is a group measure for S vs. M, and M vs. L-cones, defined as the area of overlap between their respective component Gaussians. Probability is a measure specific to each cone within a cluster, defined as the ratio of its component Gaussian value to the sum of all Gaussians at a given x-axis location (ΔOPL). High or low probability is observed when a given cone’s ΔOPL lies near or distant from the fitted Gaussian respectively.

2.7 Comparison of cone spectral classification across different days to assess test-retest reliability

To assess test-retest reliability, cone classification was conducted on two different days for two subjects, at one eccentricity each. Subjects S2 and S4 were classified at 1.5° and 4° eccentricity respectively. The recorded OCT volumes on both days were registered together to select the same cones for classification. Cones so-classified were compared for cone-type specific mismatches between days. Bland-Altman analysis [28] was performed to compare the maximum ΔOPL across both days to evaluate their repeatability.

2.8 Cone-by-cone alignment for comparing densitometry and ORG

To compare outcomes of the two classification methods, we were first required to identify matching cones between the datasets obtained in AOSLO densitometry and AO-OCT based ORG. However, the two datasets obtained from these different instruments feature significantly different scales, orientations, cone appearances and motion-induced distortions. To address this problem, we used a specially designed Retina Map Alignment Tool. This tool aligns two cone-labeled retina images (‘retina maps’) via two-stages: automatic coarse image registration followed by human-assisted refinement and verification. This tool allows us to match hundreds to thousands of cones per subject with confidence in a fraction of the time it would take manually. See Fig. 3 for a detailed step-by-step overview of this pipeline with accompanying visualizations. The remainder of this subsection provides additional technical detail of the tool’s implementation.

 

Fig. 3. Description of the Retina Map Alignment Tool to find matching cones between an SLO image and the OCT en face image at the COST. In this example, we refer to a cone’s coordinate and LMS classification as its ‘label,’ and show each label as a false-color red, green or blue spot. We use the tool to warp the labeled SLO image such that the cones’ coordinates align with the OCT image. This warping then provides the desired cone pairing between the two datasets. The steps involved are shown left-to-right going from inputs to outputs. The bottom row shows a zoom-in of the SLO labels as they are gradually warped to align with the cones visible in the OCT image. From left to right: 1a) The inputs are the cone-labeled SLO image and the destination OCT image. 1b) Naively copying the SLO label positions directly onto the OCT image fails to pair the labels with the visible cones, due to their relative image distortion. 2a) We first rigidly match labeled patches of the SLO image to get a rough alignment with the OCT image. 2b) Notice how many SLO labels fall in the gaps between the cones, making pairings unclear. 3a) The center of each patch is converted into a vertex in a flexible triangle mesh. This mesh defines a piecewise linear continuous warping of the labeled SLO image. The user adjusts this mesh interactively by adding, deleting, or moving vertexes of the mesh. Shown in cyan are the vertexes and triangle mesh after user adjustment. 3b) Each label now overlaps with its paired cone in the OCT image. However, notice how many of the labels are still not centered on the cone cells. This makes visual verification of the pairing unnecessarily difficult. 4a) To fix this, we automatically move each label to the center of its paired cone cell in the OCT image. This is done by attaching a new mesh vertex to each label and moving this vertex to the nearest cone center in the OCT image. Shown in cyan are the resulting vertexes and triangle mesh. 4b) Notice that every label is now centered with each cone. 5a + b) The tool outputs the final SLO labels aligned to the OCT image, along with the warped SLO image. The user can now visually verify the pairing by checking cone-alignment between the warped SLO labels and the destination OCT image.

Download Full Size | PDF

The automatic stage of the Retina Map Alignment Tool finds a piecewise linear warping between the SLO image and the OCT en face image at the COST [29]. First, we find a rough translation, rotation, and scale between the images via the method by Reddy and Chatterji [30]. Then, the SLO image is split into 0.1° wide patches and separately matched against the OCT image via normalized cross-correlation [31]. Finally, Random sample consensus (RANSAC) is used to filter outlier matches [32]. The center of each patch gives a rough correspondence point between the two images. To obtain a continuous warping, a triangle mesh is constructed with a vertex at each correspondence point. The tool then computes per-triangle transforms and applies texture mapping techniques to warp the SLO image onto the OCT image.

The human-assisted stage enables the user to adjust the image warp interactively, and visually verify the accuracy of the final cone matches. The user inspects the initial image warp alignment by flipping between the two registered images and checking that each cone visually aligns with its match. If further adjustment is needed, the tool allows the user to add, delete, or move correspondence points and see the updated warp in real time through the tool’s graphical user interface. Once the user verifies SLO image warp, cone matches are identified by pairing each cone with its nearest neighbor across the images. Individual matches can also be manually adjusted. Finally, the tool outputs a file with the cone coordinates and identities relative to the new image.

The human assisted stage was conducted by two graders. Cone matches that agreed between the graders were included in analysis. Disagreement primarily occurred at the edges of the regions-of-interest, where there is less information to ascertain the match. Some cones matched were clearly visualized in one image but aligned to a dark or abnormally reflective cone in the corresponding image and surrounded by clearly matched cones. These cones were dubbed “dysflective cones” (see section 3.3 below), and likewise agreed upon between the graders. Overall, the disagreement between the graders led to the exclusion of 12 of 854 cones for subject 4 and 2 of 567 cones for subject 5 from the analysis.

3. Results

3.1 Complex average vs. angle average

In this section, we compare the two methods of averaging – complex and angle average – across repeat measurements. Figure 4(a) shows an example of a cone where the individual trial ORG traces have low inter-trial variability, low noise and no unwrapping artifacts. In this case, the angle average and the complex average traces align well with each other. Figure 4(b) shows an example of a cone where repeat measurements are noisy and have unwrapping artifacts. These factors corrupt the angle average trace, and in this specific case, lead to a lower overall ΔOPL change compared to the complex average. This comparison was repeated for the same population of cones (n = 512) in Figs. 4(c) and 4(d), where there is expected to be a mixture of cones with high and low inter-trial variability. The complex average of the constituent trials shows a clear separation between the three cone types in Fig. 4(c). In comparison, the angle average contains more noise in the individual cone ΔOPL versus time traces leading to an ambiguity in the cone assignments. We follow the complex average for all data presented in this article. Phase variance is directly related to the inverse of signal-to-noise – higher the SNR, lower the phase noise and vice versa. The complex average effectively weights the phase information with the signal amplitude and is more robust to noise artifacts including phase wrapping.

 

Fig. 4. Comparison between complex and angle average for calculating the mean ΔOPL vs. time. Figure 4(a) and 4(b) show examples for 10 repeat trials for 2 different cones in cases where the individual trials have low (a) and high (b) inter-trial variability respectively. When the two approaches are applied to a population of the same 512 cones, the complex average (Fig. 4(c)) readily segregates the cones into its respective spectral types. Angle average (Fig. 4(d)), on the other hand, is more susceptible to inter-trial variability and leads to ambiguity in the cone assignments.

Download Full Size | PDF

3.2 Separating LMS cones

Figure 5 shows a representative example of the procedure for cone classification under a 660 nm stimulus bleach at 1.5° temporal eccentricity in subject S4. The three cone types are well segregated owing to the variable bleach in each type. The absolute values of the registered AO-OCT complex volume segmented at the ISOS and COST (${I_{ISOS}}({x,y} )$ and ${I_{COST}}({x,y} )$) provides en face images revealing cone photoreceptors at the two layers (Figs. 5(a) and 5(b)). The mean ΔOPL between the ISOS and COST at time t = 0.58 sec after stimulus onset is shown in Fig. 5(c) for each cone (see Visualization 1 for time evolution of ΔOPL in individual cones). The magnitude of ΔOPL for each cone in the image is color coded according to the color bar, and shows the variability and distribution between cones. Three colors – red, yellow and blue on the color bar are dominant while the intervening colors are sparse, suggesting the correspondence to the three cone types. The time series of ΔOPL for individual cones is shown in Fig. 5(d). The maximum ΔOPL for every cone is plotted as a histogram in Fig. 5(e). The maximum ΔOPL for each cone was obtained after complex averaging the repeat trials. The histogram is separated into the corresponding clusters of L, M and S cones by the intersection of the three component Gaussians. The individual cone ΔOPL traces are shaded red, green and blue in Fig. 5(f) according to their cluster assignment in Fig. 5(e). The mean of the maximum ΔOPL for LMS-cones was 585 ± 5, 157 ± 20, 9 ± 18 nm respectively. After clustering assignment, cones that have noisy ΔOPL traces and wrapping artifacts are filtered manually by comparing their mean ΔOPL trace across all repeats, to the mean of all cones within that cluster. Within the filtered cones, individual repeat trials with excess eye movements or noise are removed, and the complex average is re-computed with the remaining trials. After these steps, the cones that exhibit normal ORG characteristics, as defined by the mean of the three clusters, are included in the final analysis. The remaining are flagged as having insufficient fidelity to unambiguously assign their type. These cones are excluded from the classification and labeled as ‘not classified (NC)’. The number of NC cones is marked in Table 1.

 

Fig. 5. Representative example for cone classification for Subject S4. Registered AO-OCT volumes segmented at ISOS and COST yields the en face cone photoreceptor images at these layers in (a) and (b) respectively. The mean ΔOPL of seven recordings in cone outer segments at time t = 0.5 is color coded according to the color bar in (c) and overlaid on the cone mosaic in the COST image. Individual cone time series of ΔOPL for all the cones identified in the AO-OCT structural image is plotted in Fig. 5(d). A histogram of the maximum ΔOPL averaged over the time after stimulus onset between 0.7–1 sec is shown in Fig. 5(e). Gaussian mixture model clustering analysis is used to segregate the cone types based on the maximum ΔOPL in the histogram. Figure 5(f) shows the same traces from Fig. 5(d), now color coded as red, green and blue to represent LMS-cones respectively. Scale bar: 10 arcmin

Download Full Size | PDF

Table 1. Cone Classification Summary for 3 color normal subjects at 1.5° eccentricity

View Table | View all tables in this article

Figure 6 shows the cone classification for three subjects (S1, S2, and S3, from left to right columns) following the procedure in Fig. 5. The x,y cone locations in arc-min, and spectral types are provided in Data File 1[33]. Figures 6(a,b,c) shows the gray scale image and (d,e,f) shows the histogram of ΔOPL and cluster separations into the three cone types. Figures 6(g,h,i) shows the corresponding LMS-cone mosaics labeled in false-color according to their segregation obtained from Gaussian mixture model clustering analysis. Table 1 shows the total number of cones classified, percentage of S-cones, L:M ratio, the variation of saturated ΔOPL in each cone subtype, the probability and error of assignment for the 3 subjects. The remaining two subjects appear below in section 3.5. The average maximum ΔOPL between all 5 subjects (3 subjects from Table 1 and 2 subjects from Table 2) for L, M, S-cones was 648 ± 67, 170 ± 28, 15 ± 15 nm respectively. In the 5 color normal subjects, the L:M cone ratio ranged between 1.5 - 2.4, and S-cone percent between 4.8–8.1%, and both are comparable to the literature [9,34].

 

Fig. 6. Cone spectral classification in three color normal subjects using ORG. Subjects S1, S2 and S3 are shown from left to right columns. (a-c) Gray scale en face images obtained by segmenting the registered and averaged AO-OCT volume at the COST (d-f) Histograms of maximum ΔOPL calculated in the time window 0.7–1 sec after stimulus onset overlaid with the component Gaussian fits obtained from Gaussian mixture model clustering analysis. (g-i) LMS-cone mosaics false-colored as red, green and blue to represent the three cone types. Scale bar: 10 arcmin

Download Full Size | PDF

3.3 Dysflective cones

Cone photoreceptor reflections are variable, potentially due to a combination of interference of multiple reflections from the outer segment, disc shedding and renewal, abnormal waveguiding and perhaps other unknown factors [7,8,20,21,35,36]. In addition, it is also known that a visible stimulus causes reflectivity variations that can be analyzed to infer their function, also termed intensity ORG, or iORG [37–39]. Support that these variations are caused due to interference from multiple reflections from the outer segment comes from Cooper et al.2020 [38]. They used a long coherence light source greater than the length of the outer segment that reduced the variability in cone reflectivity considerably. In the extreme case, some cones have abnormally low reflectivity, dubbed ‘dysflective’ [40,41]. It has been noted that normal eyes with no known disease have patches of dysflective cones that sometimes regain their reflectivity over time [42]. Also, using retinal tracking and targeted stimulus delivery, it has been shown cones with low reflectivity exhibit normal perceptual sensitivity [43]. In this section, we sought to investigate whether any potential functional deficit in ORG can be observed in dysflective cones.

Image formation in OCT relies on optical interference, and it may be plausible that reflections from some cones undergo destructive interference to exhibit abnormally low reflections. In the AOSLO, cones lying under the shadows of blood vessels have lower apparent reflectivity, especially when imaged with a visible wavelength as in retinal densitometry. Examples of dysflective cones are shown in Fig. 7 in subject S4 (a) and S5 (b) in their respective AO-OCT en face images of the COST. With AO-OCT alone, it was not possible to ascertain whether cones were indeed present in these dark areas indicated by the yellow circles. However, with the ‘Retina Map Alignment Tool’ described in section 2.8, we observed cone reflections present in the corresponding AOSLO image from the same subjects at these locations (S4: n = 2 cones, S5: n = 3). Figure 7(c) and 7(d) show the ΔOPL versus time for the dysflective cones indicated in Fig. 7(a) and 7(b) after a 660 nm stimulus bleach. For these dysflective cones, the phase difference was calculated at the same axial locations as where ISOS or COST would have been had there been a reflection from them. This axial location was determined based on the neighboring cones’ reflections. This approach assumes that the axial reflection locations do not change substantially in a local area. These traces indicate normal ORG characteristics consistent with their classification into either L-cones or M-cones, based on their high and low maximum ΔOPL respectively. This indicates that, in addition to having normal perceptual sensitivity, dysflective cones exhibit normal functional responses in the ORG and transduce light with similar efficacy as cones with normal reflectivity.

 

Fig. 7. Optoretinography in dysflective cones. Figure 7(a) and (b) show AO-OCT grayscale images for subject S4 (a) and S5 (b) at the COST layer where yellow circles indicate cones whose reflections were very low or absent. That cones were present in these ‘dark’ spaces was confirmed by finding their locations in the corresponding AOSLO image in the same subjects using the ‘Retinal Map Alignment Tool’ (see section 2.8). Dark shadow in 7b represents blood vessels overlying the cones. Figure 7(c) and (d) show the corresponding ΔOPL vs. time traces showing normal characteristics consistent with their classification into L-cones (Fig. 7(c)) and LM-cones (Fig. 7(d)). Scale bar: 5 arc-min.

Download Full Size | PDF

3.4 Test-retest reliability in ORG cone classification

Figure 8 evaluates the test-retest reliability of classification across two different days. Figures 8(a) and (c) show the maximum ΔOPL across the two days, labeled ‘day 1’ and ‘day 2’, at 1.5° and 4° eccentricity in Subjects S4 and S2 respectively. The cones were unambiguously selected on the images at 4° eccentricity without intrusion from rods (thanks to the OCT axial resolution) for the ORG measurements. The maximum ΔOPLs obtained in each cone for the two days are plotted against one another; each data point represents an individual cone. Three distinct clusters with increasing maximum ΔOPL are observable corresponding to the LMS-cones. Cones were classified with high fidelity using Gaussian mixture model analysis based on each day’s measurements alone, with errors less than 0.02% and mean probability of 0.99. The distribution of data points in both x and y dimensions represents the variability in ΔOPL for each day, increasing in magnitude with the maximum ΔOPL. A few data points in each case lie outside the clusters and beyond the equality line. These correspond to the cones that did not match across days in their magnitude of maximum ΔOPL. Phase wrapping errors that lead to 2π phase ambiguities can partly explain these mismatches where a few cones exhibit high maximum ΔOPL on one day and low on the other. The mean Euclidean distance from the origin for the mismatched cones was 340 nm; in comparison, for the OCT center wavelength of 820 nm, a 2π ambiguity would result in a 410 nm error in ΔOPL.

 

Fig. 8. Comparison of maximum ΔOPL across two days, labeled day 1 and day 2, in the same cones. Figure 8(a,b) and (c,d) show data at 1.5° and 4° eccentricity for S4 and S2 respectively. Note the different y-axis scale at the two eccentricities. In a) and c), the magnitude of maximum ΔOPL for days 1 and 2 is plotted against each other, b) and d) show Bland-Altman plots for difference vs. mean of the maximum ΔOPL across the days. Horizontal orange line indicates ±95% confidence intervals in the difference between days, and the horizontal black line indicates the mean difference or bias.

Download Full Size | PDF

To analyze the agreement in maximum ΔOPL as a whole between days, Figs. 8(b) and (d) show a Bland-Altman plot for 1.5° and 4° eccentricity respectively. Each data point is the measurement from the same cone across days. The three clusters corresponding to the cone types and the greater magnitudes of variability with increasing maximum ΔOPL are also observed in these plots (also see the standard deviation for S, M and L-cones in Tables 1 and 2). Typically, Bland-Altman plots represent the relationship between the difference (y-axis) versus the mean (x-axis) of two measurements to gauge their agreement. At 1.5° and 4° eccentricity, 65 out of a total of 1397 cones and 13 out of a total of 305 cones lie beyond the ±95% confidence interval (orange line) respectively. The bias (black line), equal to the mean difference across all cones for both days, is meant to indicate any potential systematic difference between the days. The values of bias were 8.06 and 2.45 nm for 1.5° and 4° respectively, and is close to the noise floor of the measurement for the ORG paradigm in the line-scan OCT [16]. At 1.5° eccentricity, of the 1397 cones classified across both days based on Gaussian mixture model analysis on maximum ΔOPL, 4.2% (59 cones – 48 L/M, 9 L/S, 2 M/S) were inconsistent in their assignment between the days. At 4° eccentricity, 305 cones were classified across both days, of which 2.6% (8 cones – 3 L/M, 0 L/S, 5 M/S) were inconsistent in their assignment between days. Overall, a high agreement was noted between assignments across different days. The slightly larger inconsistency at 1.5° could be related to phase wrapping errors. It is likely these were not sufficiently eliminated due to inadequate repeat measurements - 6 repeats on day 1 and 5 repeats on day 2 at 1.5° - compared to the 10 repeats on each day at 4°. Subject S4 has good fixation and typically has a lower probability of individual trial rejection due to large eye movements; hence we recorded 6 and 5 repeats on day 1 and day 2. No repeat trials were removed from analysis. Another potential reason could be the higher cone density and smaller cone sizes at the closer eccentricity that may pose greater noise in the measurements across days.

3.5 Densitometry vs. ORG classification comparison

We observed high classification agreement using densitometry and ORG ∼6-7 years apart in Subjects S4 and S5 in Figs. 9 and 10 (S4: 771/842 cones, 91.6%; S5: 511/565 cones, 90.4%). Figure 9 shows the false-color LMS-cone mosaics from S4 and S5, classified at the same retinal eccentricity in densitometry and ORG. The x,y cone locations in arc-min, and spectral types obtained via ORG and densitometry are provided in Data File 1, [33] for both subjects. After cone matching and visual inspection by 2 separate graders, a few cones (S4: 12; S5: 4) present in Fig. 9 were misaligned to a degree not amenable for comparison, and these were removed from further comparative analysis. Regions classified with ORG fully encapsulated regions classified with densitometry, and all cones classified with densitometry (S4: n = 830, S5: n = 559) were classified with ORG. If the comparison were restricted only to cones that were successfully classified in densitometry, the proportion agreement in S4 and S5 was 92.9% (771/830 cones) and 91.4% (511/559).

 

Fig. 9. Cone spectral classification with densitometry and ORG. Figure 9(a) and (b) show the cone spectral assignments in false color in S4 for densitometry and ORG, respectively overlaid on a grayscale image. Figure 9(c) and (d) show the same for S5. Scale bar: 5 arc-min.

Download Full Size | PDF

 

Fig. 10. Classification mismatch between densitometry and ORG. Figure 10(a) and (b) show ORG-classified mosaics highlighting densitometry-mismatches. Cone color indicates ORG LMS-classification. Cones with mismatched classification are plotted larger with a halo that indicates their densitometry classification. Figure 10(c) and (d) show classification confusion matrices which quantify the number of cone-type specific matches in ORG and densitometry. Colors of each square and the corresponding color bars indicate the proportion matched between both modalities. Figures 10(a,c) are S4 classification mismatch and confusion matrix. Figures 10(b,d) are S5 classification mismatch and confusion matrix.

Download Full Size | PDF

A small number of cones were not visible in the AO-OCT image (dubbed dysflective cones as section 3.3 above), but were nonetheless classified thanks to their transferred location from the densitometry image (see ‘Retina Map Alignment Tool’ in section 2.8 and Fig. 7) (S4: n = 2, S5: n = 3). Finally, several cones which were not classified (NC) with densitometry owing to low reflections, in part caused by their lying under shadows of blood vessels, were classified with ORG (S4: n = 12, S5: n = 6).

A cone-by-cone summary of classification agreement across the two paradigms is summarized in Figs. 10(a-d). Figures 10(a) and (b) shows the distribution of cones across the regions of interest and how their cone assignment is consistent or different between the two modalities in S4 and S5. The red, green and blue dots in the foreground indicate the cone assignments obtained from ORG. Comparing the assignments as a whole in this manner helps appreciate potential factors related to the spatial arrangement or specific cone-types that may lead to inconsistencies. For example, it has been suggested that optical blur in densitometry might lead to an error in classification where cones of like-type are clumped together [4,9]. Visual inspection shows few examples in both subjects where patches of 3-6 contiguous cones are mismatched. Two additional features are worth noting. Cones lying near the edge of the region-of-interest have a higher likelihood for mismatch than the center. One potential reason for this is eye motion that reduces the number of samples averaged at the borders compared to the center. Second, the bottom left quadrant of the region-of-interest in subject S5 has higher degree of mismatch, a location that has blood vessels overlying the cones casting a shadow, as shown in the gray scale image in Fig. 7(b). This causes decreased back reflectivity of cones lying in the shadow, increased absorption by the blood vessels for the bleaching light and increased forward and backward scatter of the imaging light in both modalities.

Confusion matrices in Figs. 10(c) and 10(d) summarize the mismatch by cone type in both modalities. In both subjects, classification agreement (relative to ORG assignment) was highest for L-cones (S4: 92%, S5: 97%), decreasing for M-cones (S4: 90%, S5: 88%) and S-cones (S4: 93%, S5 = 59%). Table 2 summarizes the classification agreement, L:M ratio and % S-cones measured from each method. Mismatches had a slight effect on the L:M ratio measured for both subjects (S4: 2.0 vs. 2.2, S5: 2.0 vs. 1.6), negligible effect on the percent S-cones measured for Subject S4 (4.6 v 4.8%), and a more substantial effect on percent S-cones measured for Subject S5 (5.8 vs. 8.1%). For Subject S4, classification mismatches were dominated by L/M mismatches (of 71 mismatches, L/M = 76.1%, S/(LM) = 7.0%, NC = 16.9%), whereas a higher proportion of mismatches for subject S5 involve S-cones (of 54 mismatches, L/M = 44.4%, S/(LM) = 44.4%, NC = 11.1%). In particular, for Subject S5, a disproportionate number of cones assigned S in ORG were assigned L with densitometry (17/54 mismatches). Of these 17 L vs. S-cone mismatches, 5 occur at the borders of the region-of-interest, and 7 occur in the area overshadowed by the blood vessel (bottom left quadrant of Fig. 10(b), corresponding AO-OCT structural image at COST layer with blood vessel shadow is shown in Fig. 7(b)).

Table 2. Summary of comparison between densitometry and ORG for cone classification

View Table | View all tables in this article

4. Discussion

This article presents the optimization and characterization of cone spectral classification using ORG and its comparison against densitometry. The cone mosaics of five color normal subjects were classified with ORG, showing high probability (∼0.99), low error (<0.22%), high test-retest reliability (∼97%) and short imaging durations (< 1 hour). The values for error are comparable to those reported previously for cone classification with point-scan AO-OCT [13]. The L:M cone ratio, and percent S-cones measured in this cohort are consistent with literature [9,34]. In comparison, slightly lower probability (∼0.95) and higher error (∼ 5%) was observed in densitometry for classifying cones [11].

Some caution is suggested in interpreting error rates and probability of cone assignments obtained from clustering analysis shown here and in previous reports. An approximate analogy could be made to the concepts of ‘precision’ and ‘accuracy’ in this regard. The error rates and probabilities as derived from the clustering analyses are very useful in indicating the precision by which the segregation of cone types can be obtained by a particular technique, when the task is faced with variability related to biological and methodological factors. However, accuracy, loosely defined in the current perspective, relates to how closely the cone types so-obtained relate to the “ground truth”. When the ground truth is unknown as is the case here, repeat measurements are carried out to assess the degree to which the results remain resilient amidst the myriad sources of variability. To this end, the test-retest analysis shown here evaluates the degree to which cone assignments remain uniform across days and can be considered the upper bound error estimate of the accuracy of classification. Such an analysis has not been reported for densitometry.

The laboratory time taken for ORG data collection ranges from 30-50 min per one region-of-interest. This includes the optimization of image focus in line-scan ophthalmoscope and OCT en face images (15-20 min) and 10 repeat measurements (15-30 min). This is significantly more efficient compared to densitometry that could take 3-9 hours of laboratory time per region-of-interest [11]. Further, it was suggested that the imaging sessions be distributed across a few days to account for reflectivity variations in cones, independent of photopigment absorption, leading generally to a greater time needed per subject [9]. Averaging is essential for both modalities given how ORG cone assignments can also be affected by phase wrapping errors when there are inadequate repeat measurements. The data processing in MATLAB for ORG takes about 4-5 hours per subject where ∼50% of this time is consumed for reading and writing the data between the solid-state hard drive and MATLAB workspace. This data processing time is not optimized, and it can be reduced drastically by using C++ with CUDA.

The overall higher fidelity of ORG is the result of greater than 20-fold dynamic range available in ORG compared to densitometry for capturing the perturbations caused by bleaching in the backscattered reflection from cones. Since the OCT layer reflections encasing the outer segment are used for computing the light-evoked ΔOPL, the optical signal is largely confined to the cone outer segment and not corrupted by reflections from other layers. Optical phase difference used to compute the ΔOPL has high sensitivity of about a few nanometers when phase is referenced between the ISOS and COST to remove common-mode noise. The bleach dependence of the ΔOPL is non-linear, scaling in magnitude with a gain greater than 1 suggesting a relationship to the amplification cascade of phototransduction [14]. The mechanism of the change in outer segment OPL has been attributed to the influx of water to maintain osmotic balance during phototransduction [24]. This is suggested to be in response to the disturbance of osmotic balance caused by the intermediary osmolytic products of the cascade that are amplified in their concentration compared to photopigments. This is critical to recognize when using ΔOPL as a biomarker for retinal diseases where mutations in specific proteins of phototransduction are implicated. In such cases, ΔOPL will not directly provide a gauge for photopigment concentration. On the other hand, the ORG early response, that manifests as a reduction in ΔOPL immediately after stimulus onset and occurs prior to the onset of the amplification cascade, can provide a direct measure of photopigment activation [14,44]. A few examples of early response maximum amplitude have been observed (unpublished data) that vary in a manner consistent with their classification into the 3-cone types. But, the dynamic range from ∼5 - 50 nm and the associated noise at these low OPL values are not sufficient to reliably differentiate the three cone types based on the early response alone. Furthermore, the early response peak amplitude occurs between 3 ms to 10 ms depending on bleach strength. To resolve this temporally, volume rates of at least 333 Hz would be required, corresponding to approximately 133k −200k B-scans/sec, in order to achieve the 400-600 B-scans per volume acquired here. Currently, this speed and SNR is beyond the reach of the instrument.

These aforementioned benefits of ORG do not apply to densitometry. Since absorption imaging is not depth-resolved to cone outer segments, it is confounded by the contribution of stray light due to scattering from other retinal layers, blood flow and anterior optics. With the confocality provided by AOSLO, the contribution of stray light can be reduced substantially. Still, maximum measurable change in absorption is limited to 0.41 log units (2.6-fold) at 550 nm [18]. Changes in optical intensity due to absorption cannot capitalize on the benefits of phototransduction amplification. The variations in cone reflectivity following a stimulus, as observable in iORG, have been used to separate S-cones [38] and can in principle be used to infer all three types.

Given that the mechanism of densitometry is well established and attributed directly to photopigment concentration, it was essential to compare the cone assignments obtained by the ORG method against it. A high agreement (mean: 91%) was observed between the two modalities in 2 subjects, with measurements conducted 6-7 years apart. Cones lying in the shadows of blood vessels in AOSLO were likely to be inconsistent in their assignments. Contiguous clusters of 3-6 mismatched cones were observed, attributable to optical blur and scatter. A greater degree of mismatch in S-cones was observed, likely in part due to the fact that S-cones are identified by a lack of their response to the bleaching stimulus in both modalities. Instead, a bleach that activates S-cones could be used in the future. Dysflective cones were found in both modalities. Cones dysflective in ORG were located in their counterpart image, and analyzed for their spectral type using ORG. Further work is needed to ascertain how the reflectivity in these cones changes with time across shorter timespans, i.e. repeat trials within a day, and longer timespans, i.e. different days. Also, it is of interest to gauge if or whether the ORG characteristics follow the changes in backscattered intensity. Such cones with abnormal reflections have been observed in healthy and diseased eyes [41,42,45]. They have previously been shown to have normal visual function also [43,45]. That they exhibit normal functional responses in the ORG has important clinical implications. A preponderance of ‘dark’ or dysflective cones in pathological eyes may not provide the complete picture of the dysfunctional state, and a functional measure of phototransduction such as ORG may be a more sensitive biomarker.

In summary, the protocol for cone classification was characterized using ORG and compared against densitometry. The benefits obtained in ORG from its higher sensitivity and dynamic range led to high efficiency and fidelity for cone spectral classification.