Optoretinography of individual human cone photoreceptors

: Photoreceptors mediate the first step of vision, transducing light and passing signals to retinal neurons that ultimately send signals along the optic nerve to the brain. A functional deficiency in the photoreceptors, due to either congenital or acquired disease, can significantly affect an individual’s sight and quality of life. Methods for quantifying the health and function of photoreceptors are essential for understanding both the progression of disease and the efficacy of treatment. Given that emerging treatments such as gene, stem cell, and small molecule therapy are designed to operate at the cellular scale, it is desirable to monitor function at the commensurate resolution of individual photoreceptors. Previously, non-invasive imaging methods for visualizing photoreceptor mosaic structure have been used to infer photoreceptor health, but these methods do not assess function directly. Conversely, most functional techniques, such as ERG and conventional microperimetry, measure function by aggregating the effects of signals from many photoreceptors. We have previously shown that stimulus-evoked intrinsic changes in intensity can be measured reliably in populations of cone photoreceptors in the intact human eye, a measurement we refer to more generally as the cone optoretinogram. Here we report that we can resolve the intensity optoretinogram at the level of individual cones. Moreover, we show that the individual cone optoretinogram exhibits two key signatures expected of a functional measure. First, responses in individual cones increase systematically as a function of stimulus irradiance. Second, we can use the amplitude of the functional response to middle wavelength (545 nm) light to separate the population of short-wavelength-sensitive (S) cones from the population of middle- and long-wavelength-sensitive (L and M) cones. Our results demonstrate the promise of optoretinography as a direct diagnostic measure of individual cone function in the living human eye. densitometry. We also performed the same classification on the log 10 iORG amplitudes obtained in response to the 450 nW/degree 2 stimulus.

Our previous report of iORG signals revealed that the human cone iORG is sensitive to stimulus irradiance and that its action spectrum is well-matched to that expected from human cone spectral sensitivities [18]. These measurements were population-based, with the signals aggregated from multiple cones over a region of retina, often subtending hundreds of cones. While such aggregation enhances the robustness of the cone iORG assay, population-based measurements do not enable the assessment of function at the fine spatial scale of individual photoreceptors.
Despite prior observations of an iORG signal, quantification of the stimulus-evoked changes in an individual cone's intensity remains elusive. Here we report that we can resolve the iORG at the level of an individual photoreceptor. Moreover, we show that the individual cone iORG exhibits two key signatures expected of a functional measure: First, responses in individual cones increase systematically as a function of stimulus irradiance. Second, we can use the amplitude of the iORG in response to middle wavelength (545 nm) light to separate the population of shortwavelength-sensitive (S) cones from the population of middle-and long-wavelength-sensitive (L and M) cones. Our results demonstrate the promise of intensity optoretinography as a direct diagnostic measure of individual cone function in the living human eye.

Pre-registration
The experiments were subject to pre-registration submitted to the Open Science Foundation (https://osf.io/jtmx6/). We generally followed the protocols as written, with some adjustments made as we gained experience carrying out the experiments and analysis. The number of subjects studied was fewer than pre-registered because of limitations on human subjects studies imposed by the Covid-19 global pandemic. The pre-registration history during project development as well as key deviations from the final pre-registered protocol are listed in Table S1 in Supplement 1.

Human subjects
This research was approved by the Institutional Research Board at the University of Pennsylvania, and was conducted in accordance with the tenets of the Declaration of Helsinki. Two subjects with no known pathology were recruited for this study. Subjects provided informed consent after the nature and possible risks of the study were explained. Each subject's pupil was dilated and accommodation arrested using one drop each of tropicamide (1%) and phenylephrine (2.5%).

Intrinsic imaging of the photoreceptor mosaic
Subjects were imaged using a custom adaptive optics scanning laser ophthalmoscope ( [13]; AOSLO). To minimize imaging during disc shedding events, which have been reported to occur at their highest rate in the morning [29,30], subjects were required to have been awake for at least 4 hours before the imaging session. Each image sequence was obtained 1.5-2°from the subject's center of fixation at a rate of 17.85 frames per second, using a single-mode 785nm laser diode (Thorlabs LP785-SF20, 2nm bandwidth, 95µm coherence length) and a 1 × 1°field of view. We acquired images from three separate channels simultaneously, allowing us to obtain a confocal imaging mode alongside a split-detection mode [31]. The cone intensity optoretinogram was obtained using the confocal mode, while the split-detection mode was used as the primary sequence for the image registration process described below. This decision was made due to the significant stimulus-induced intensity changes observed in the confocal image sequence, which caused registration to fail. The intensity of structures in the split-detection mode was more stable, making it a better choice as the primary sequence for registration.
To describe the study, we adopt the following terminological hierarchy (from lowest to highest-level): • Acquisition: A single recording using the AOSLO.
• Trial: A set of 13 acquisitions, using a single stimulus condition.
• Permutation: A randomly ordered group of trials over the set of stimulus irradiances used in the particular session, where each trial had a unique stimulus irradiance.
• Session: A single visit by a subject, which consisted of 5 permutations.
Each acquisition consisted of 72 images of prestimulus recording (approximately four seconds), 18 images of recording during visible light stimulation when a 545nm stimulus (NKT Photonics, SuperK EXTREME/VARIA, 10nm bandwidth) was delivered to the entirety of the imaging field (approximately one second), and 72 images of post-stimulus recording (approximately four seconds). The subject then closed his or her eyes for approximately nine seconds to rest. In total, each acquisition lasted approximately 18 seconds, with 162 images acquired over the first nine seconds. Each trial was preceded by two minutes of dark-adaptation and consisted of exactly 13 acquisitions. We acquired five trials for each of five stimulus irradiances: 0 (control), 25, 50, 140, and 450nW/deg 2 . The control and 450nW/deg 2 stimulus trials were obtained in different imaging sessions than the 25, 50, and 140nW/deg 2 stimulus trials. In total, there were 325 acquisitions obtained per subject (13 acquisitions/trial * 5 trials/stimulus * 5 stimuli). Additional trials using higher stimulus irradiances were obtained within the permutations for the same sessions as the 0 and 450nW/deg 2 stimuli. Data from these trials are not reported here as we concluded that these higher irradiances (e.g. >900nW/deg 2 ) led to a saturated response. Exploratory analysis suggested that the signals obtained from the first three acquisitions of a trial were inconsistent with subsequent acquisitions, possibly due to a larger loss of pigment per stimulus, before pigment loss per stimulus approached a steady state ( Figure S1 in Supplement 1); thus, the first three acquisitions of each trial were excluded from the analysis. Each 1 second stimulus presentation corresponded to an average 0%, 2.0%, 3.5%, 7.0%, and 10.5% bleach of the remaining photopigment in L and M cones, at the five stimulus irradiances respectively ( Figure S1 in Supplement 1). The exclusion of the first three trials was pre-registered in advance of the data collection reported here.

Image sequence post-processing
After acquiring the data, a reference image was automatically selected from the split-detection modality within the stimulus delivery range ( [32]; selected from image indices 73-90). Next, we corrected the static intra-frame distortion caused by the sinusoidal motion of the resonant optical scanner. This was done by estimating the static spatial distortion from separately acquired images of a stationary vertical Ronchi ruling and resampling each image onto a grid of equally spaced pixels. All of the "desinusoided" images in the acquisition were then strip registered using a previously described algorithm [33] to the selected (and also desinusoided) reference image. If the chosen reference image or any image within indices 73-90 was unable to register with the rest of the sequence, or if the images were not of the intended retinal location because of inordinate eye movements, then the acquisition was removed from further consideration. Additionally, if fewer than 60% of all images or fewer than 60% of the images within the stimulus delivery range (image indices 73-90) were able to register to the reference image, then the acquisition was removed from further consideration. Out of the 325 acquisitions obtained for each subject, 73 and 9 acquisitions from subjects 11002 and 11049 respectively were discarded during this stage of processing.
Following strip registration, each acquisition was processed through a fully-automated postprocessing pipeline, which consisted of the following steps: 1) Residual distortions caused by eye motion within the reference frame were removed using a custom algorithm [32]. 2) The registered acquisition was cleaned by removing frames that had registration errors mid-image, and frames that had mean intensities two standard deviations below the mean of the entire acquisition. 3) After these frames were removed, the acquisition was cropped to regions of the acquisition that always had image data present. 4) The strip-registration algorithm uses a purely translational model to correct for eye motion and does not account for eye torsion present during imaging. If torsion is present, then the registration produces errors that are manifested as a combination of rotation and skew. To remove this residual error from our acquisition, the final step was to perform an affine registration of the strip-registered, cleaned, and cropped acquisition to the cleaned and cropped reference image.
All acquisitions were then registered to the acquisition with the largest image area using constellation matching [32]. A "super average" image was created from all post-processed acquisitions, and cone locations in this "super average" image were marked and these locations used for all following analyses. Cones underlying retinal capillaries were excluded across all registered acquisitions by highlighting the vasculature with a version of a previously described motion contrast algorithm [34], modified to apply to multiple independently acquired acquisitions. Specifically, the modifications made the algorithm more robust to imperfect video overlap, image noise characteristics, and temporal discontinuities due to image exclusion and the application of the algorithm across acquisitions. A binary vasculature mask was created from the algorithm's resultant motion contrast image by thresholding the image using MATLAB's adaptthresh function. Cones falling within the mask were excluded from analysis.

Extracting the individual cone photoreceptor intensity optoretinogram
A time-varying intensity signal (intensity vs time) was created for each analyzed cone by first projecting a 3 × 3 pixel column through the aligned acquisition at each cone location and then averaging the 9 (3 × 3) pixels in each frame. For each image in the sequence, each cone's intensity for any given time point was divided by the mean intensity of all cones at that time point, thus normalizing for any overall frame-to-frame variation in image intensity. Each cone signal was then standardized with respect to its pre-stimulus mean and standard deviation as described previously, thus taking into account cone-to-cone variability that is not related to the light stimulus ( [18]; Fig. 1(B)). Finally, if a cone did not have an intensity signal for 60% or more of the stimulus delivery time of any acquisition, that cone's signal from that acquisition was discarded. This type of discard could occur due to a frame being removed due to the post-processing described above, or due to eye motion moving the cone out of the registered area. After the discard process, cones that had fewer than 25 viable acquisitions per stimulus irradiance (out of the nominal 50) were excluded from analysis.
As our AOSLO images were obtained using NIR light, we did not expect to observe a significant mean intensity change due to photopigment bleaching [15,18]. We summarized the behavior of a single cone for a given stimulus irradiance (including the 0nW/deg 2 control condition) by calculating the moving root mean square (RMS). We used a 5 frame, or 0.28 sec window (as our stimulus is coupled to the raster:  Overview of the cone reflectance analysis process. A) A cropped section of the average image from multiple confocal acquisitions that were exposed to 450nW/deg 2 545 nm stimulus, with a single cone outlined. B) The normalized, standardized reflectance of 48 trials from the single cone outlined in (A) (blue signals). This cone is typical and responds heterogeneously to visible stimulation (grey bar). C) To create an individual cone's iORG, all acquisitions from a given cone were summarized using moving average RMS with a 0.28 second averaging window. The 95 th percentile of the iORG was used as the amplitude of the summarized signal.
of that stimulus irradiance, defined as: where N is the number of elements within R a [t], and where R a [t] is the union of all standardized intensity values r centered around frame index t across all usable acquisitions a. The mean of the control (0nW/deg 2 ) moving RMS across all cones was then subtracted from the moving RMS obtained for each cone. The result was defined as the cone's intensity optoretinogram (iORG) for that stimulus irradiance ( Fig. 1(C)). Note that the amplitude of the response at a time t will be negative in cases where a cone's unsubtracted amplitude is less than the mean unsubtracted amplitude of the control condition at the same time t.
We extracted the amplitude of the iORG by determining the 95 th percentile of the signal during and around the stimulus delivery interval (image indices 71-92, approximately equivalent to the stimulus delivery interval ± (moving RMS window duration/2). This analysis approach was used for every stimulus irradiance. We found that the variance of the iORG across the population of cones increased substantially with stimulus irradiance ( Figure S2 in Supplement 1). To stabilize the variance across irradiance, we applied a log 10 -plus-one transform to the iORG, which we refer to as log 10 cone iORG amplitude.

Classifying S versus L and M cones using adaptive optics densitometry
To classify S versus L and M cones we performed dynamic densitometry using the methods described in Sabesan et al. [35]. The purpose of this classification was to enable an analysis to validate the functional significance of the iORG amplitude, as described in Results. Briefly: Each subject was dark adapted for 5 minutes prior to every densitometry acquisition. Then, 3uW of 545nm light was used to image the cone mosaic over 72 frames (approximately 4 seconds), followed by registration and alignment of the image sequences as described above. Cone intensity signals for this analysis were again extracted by projecting a 3 × 3 pixel column through the aligned image sequence at each cone location and then averaging the 9 (3 × 3) pixels in each frame.
Because L/M cones are ∼2 log units more sensitive to 545nm light than S cones [36], the protocol described above leads to selective bleaching of the L and M cones, which in turn leads to a systematic increase in light reflected from L and M cones over the duration of stimulation, since less photopigment is present to absorb the imaging light as it passes through the L and M cones. We obtained a minimum of 15 candidate densitometry acquisitions from each subject, with data collected over multiple days. Each densitometry acquisition from each cone was normalized to the 95 th percentile of its intensity values over time, providing a relative intensity measure over time for each cone. All acquisitions from each cone were combined into one dataset and fit with the exponential bleaching model: where the parameters γ, β, and τ were obtained from a non-linear least squares fit, as described by Sabesan et al. [35]. The densitometry "amplitude" of each cone was taken as the fit value at t = 0 sec subtracted from the fit value at t = 4 sec. Histograms of densitometry amplitudes for each subject were fit with a two-component Gaussian mixture model. Cones were classified into two groups based on which Gaussian component had a greater amplitude for that cone. Cones classified under the Gaussian with the lower mean were taken as putative S cones, while all others were considered putative L/M cones. The error rate was calculated using the integral measure of overlap, which is the sum of the overlapping area of each Gaussian.

Results
We obtained intensity signals from 782 and 789 cones in subjects 11002 and 11049, respectively. The average (± standard deviation) log 10 cone iORG amplitude was 0.25 ± 0.15, 0.34 ± 0.18, 0.49 ± 0.23, 0.62 ± 0.21 for subject 11002, and 0.16 ± 0.13, 0.31 ± 0.17, 0.52 ± 0.21, 0.69 ± 0.20 for subject 11049 corresponding to 25, 50, 140, and 450nW/deg 2 respectively (Fig. 2). The top panels of Fig. 3 plot (for each subject) the log 10 iORG amplitude against stimulus irradiance for each cone. On average, the iORG increases with stimulus irradiance. To quantify this, we fit the log 10 iORG amplitude versus irradiance function for each cone with a line (least squares fit; intercept locked to 0). All individual cone log 10 iORG amplitudes increased with irradiance, as evidenced by a consistently positive slope in the fit (none of the cones for either subject had a negative slope). Histograms of the slopes for each subject are plotted in the bottom panel of Fig. 3. The average (5 th , 95 th percentiles) log slope across all cones was 0.22 (0.08, 0.34) and 0.23 (0.10, 0.32) for subjects 11002 and 11049, respectively.
The four panels of Fig. 4 show the spatial distribution of subject 11049's cone log 10 iORG amplitudes in a color-coded Voronoi diagram. On average the log 10 amplitude across panels increases (becomes more yellow across panels) with increasing stimulus irradiance, however a number of cones have a consistently low response (blue) relative to their neighbors at all irradiances (white arrows, Fig. 4).
Examination of the histograms in Fig. 3 and the spatial distributions in Fig. 4 suggest two distinct groups of cones for each subject: a small group with low amplitude versus irradiance slopes and a larger group with higher amplitude versus irradiance slopes. As the 545nm stimulus would be expected to preferentially stimulate the L and M cones while minimally stimulating the S cones, this suggests that S cones consistently exhibit a low iORG, and thus a low slope. To test this hypothesis, the slopes were classified using the Gaussian mixture model approach described in the Methods for the analysis of the densitometry data (Fig. 3, blue and orange Gaussian fits, respectively) and the iORG classification was compared to the classification of putative S versus L and M cones obtained from densitometry. We also performed the same classification on the log 10 iORG amplitudes obtained in response to the 450 nW/degree 2 stimulus. Dynamic densitometry revealed two distinct populations in both subjects (Fig. 5), ostensibly corresponding to S (low densitometry amplitude; blue Gaussian fit) and L/M (high densitometry amplitude; orange Gaussian fit) cones, corresponding to 7.2% and 4.9% of observed cones and associated error rates of 3.2% and 1.4%. Cones classified as putative S cones via dynamic densitometry were compared to our iORG results (Fig. 6). We observed a good agreement between putative S cones and L/M cones identified via dynamic densitometry and putative S cones and L/M cones identified via iORG amplitude at the 450nW/deg 2 stimulus irradiance, with an agreement of 91.1% and 96.6% and S-cone population percentages of 11.0% and 5.6% for subjects 11002 and 11049, respectively, (Fig. 6, Top left, data for subject 11049 shown). Estimated error rates for these classifications were 2.8% and 1.4%, respectively.
The correspondence was improved (93.7% and 97.7%, S-cone percentages of 10.3% and 4.7% for subjects 11002 and 11049 respectively) when comparing putative S cones from dynamic densitometry to those identified from the slope of the log 10 iORG amplitude versus irradiance function (Fig. 6, Top right, data for subject 11049 shown). Error rates also improved to 2.4% and 0.77% for the two subjects respectively. Tables of classification performance (Fig. 6, Bottom  row) revealed a tendency of the iORG technique to classify more cones as S cones in comparison to the densitometry technique in subject 11002, but have decreased error.
We did not conduct formal analyses of the topography of the putative S cones. Qualitatively, they appeared evenly distributed throughout the image patches studied (see Fig. 4). The individual cone iORG amplitudes for each subject (black lines) were fit with a fixed 0-intercept linear regression line to determine their slope. Bottom row: Histograms of the log 10 iORG amplitude versus irradiance slope for each cone and each subject. With this visualization, two distinct populations emerged. To distinguish between the two populations, each distribution was fit using a two-component Gaussian mixed model. Cones were classified by determining which Gaussian had the higher probability for any given slope.

Discussion
We show that the iORG can be used to quantify the function of individual cone photoreceptors. By acquiring multiple acquisitions of the cone mosaic while stimulating the retina, we obtained signals that captured the increasing heterogeneity of cone intensity in response to visible light stimulation. We quantified this heterogeneity and used moving RMS to obtain the amplitude of the individual cone iORG, then measured the iORG of individual cones as a function of stimulus irradiance and summarized the results using the slope of the log 10 amplitude versus irradiance function. These slopes were positive for all cones measured. Overall, these results validate the use of the individual cone iORG as a cellular scale measure of cone function, and represent a significant step toward developing the iORG as a clinically relevant tool.
Our current technique has its roots in work we reported previously [18], where heterogeneity was assessed across a population of cones on each trial, rather than within individual cones across trials. In addition to the increase in resolution afforded by aggregating signals within individual cones rather than across a population of cones, the current technique differs in a few significant ways: First, in our original report, we measured the population iORG with a 13.6µm coherence length source (superluminescent diode); here we used a light source (laser) with a significantly greater coherence length of 95µm. Practically, this causes more variability in the intensity across cones in our images ( Figure S3 in Supplement 1), but also has the effect of enhancing the stimulus-evoked variation in cone intensity over time. In addition, we increased the number of trials for a given stimulus irradiance, with the goal of obtaining 50 acquisitions per cone; this number was chosen based on a prior analysis that showed a population iORG could be measured using only 50 cones [18]. Finally, we changed our analysis approach from pointwise pooled-standard deviation to window-based moving RMS; we found that the moving RMS produced more consistent results than simply computing RMS at the individual timepoints. The use of a moving window smooths out noise.
In our histogram analysis of the iORG slope with increasing irradiance we observed a bimodal distribution consisting of a small group of cones with low slopes and a large group with higher slopes. We also observed a similar distribution and grouping in the iORG amplitudes of our highest stimulus irradiance (450nW/deg 2 ). Presuming that the low-amplitude cones represented S-cones, we classified cones as putative S versus L/M cones and compared the results to a separate classification using established densitometry methods. Overall correspondence between iORG-based classifications and densitometry was good (>91%), as would be expected if the iORG assesses individual cone function (S cones would not be strongly stimulated by our 545nm stimulus). We did not expect perfect correspondence; in addition to the previously discussed sources of error in dynamic densitometry [35] and those discussed above, our instrument has a slower image acquisition rate (17.85 Hz) than the one on which the technique was developed (30 Hz). A slower frame-rate increases the sensitivity of both techniques to common imaging difficulties such as microsaccades, and poor SNR.
In addition to the good agreement between the dynamic densitometry and iORG techniques, we found that both methods had comparable estimated error rates (1.4-3.2% vs 0.77-2.8% respectively) for classifying S versus L/M cones, implying similar sensitivity. Importantly, we note that the cones classified using only the 450nW/deg 2 irradiance amplitude suffered only a ∼0.4-0.7% increase in error rate, suggesting that obtaining a cone's amplitude versus irradiance slope may not be necessary to obtain a reasonable estimate of cone type. Using an intrinsic method for cone classification like the one described here offers superior comfort and speed relative to dynamic densitometry, at least as we instrumented the two methods. While every trial of densitometry requires fixation near an extremely bright (∼3µW) stimulus over 3 seconds, our technique only requires fixation near a much more comfortable, 450nW stimulus over 1 second. Moreover, dynamic densitometry is a significantly longer process, requiring upwards of 3 times the amount of imaging time to allow consideration for sufficient dark adaptation between acquisitions.
In this study, we did not explore the minimum number of trials required to obtain a reliable single cone iORG signal. As described, our protocol obtained 65 acquisitions per stimulus irradiance and 65 acquisitions for the control condition, with 15 of the 65 discarded from further analysis. Even with an efficient imaging team and cooperative subjects, the associated imaging session for just a single stimulus irradiance and control required roughly one hour to complete. Session length itself may contribute to decreasing image quality towards the end of the session, as tear film integrity and subject attention wanes with increasing imaging time. Future work is needed to optimize the optoretinography protocol with the aim of decreasing imaging time while maintaining signal-to-noise ratio. This will become increasingly important as we move to clinical applications of this technique.
Recent work from groups using pORGs have also shown the ability to measure individual cone function, with excellent results that assess cone spectral identity [37] and health [38]. These results together with extant models of photoreceptor reflectance [39] suggest that the stimulus-evoked changes in optical path length in cone outer segments is due to the process of phototransduction. The pORG results are complementary to the work presented here. Indeed, models of the cone photoreceptor surfaces [14,39,40] reliably describe the variability of cone intensity observed in en-face AO ophthalmoscopes by representing each cone as an interferometer. These "cone interferometer" models suggest that intensity should modulate in response to the same optical path length changes observed in AO-OCT if intra-cone interference over the length of the cone outer segments is present. The 95µm coherence length of the imaging source we use here is substantially longer than the cone outer segment, and as such should be more than sufficient to produce an interferometric effect, given model assumptions. This means that the intensity optoretinogram we observe here and the phase optoretinograms described from AO-OCT likely arise from the same physiological mechanism following phototransduction in the cone. Consistent with this conjecture, the time course of the intensity and phase cone optoretinograms is similar [18,19,24,27,37].
Th current technique is readily instrumented using any adaptive optics en face imaging device, such as an AO-fundus ophthalmoscope or AO-SLO. This feature facilitates use of the technique for clinical applications, because AO-fundus ophthalmoscopes and AO-SLOs are relatively 'simple' to design, deploy, and operate, especially in comparison to AO-OCTs. This, combined with the potential for rapid automated post-processing [41,42] makes these devices natural candidates for use in clinical research environments and clinical trials. Moreover, obtaining iORGs in an extant en face system only requires use of a visible stimulus source and the algorithms described above (our reference implementation may be found at [43]).
Overall, the results presented here represent a significant step forward in the utility of the intensity optoretinogram as a functional biomarker. The AO ophthalmoscope-based optoretinogram has proven itself a sensitive tool that is able to quantify the function of photoreceptors at the scale of an individual cell, positioning it to holistically assess the progression of retinal disease and outcome of therapies moving forward.