1. Introduction

OCT is an optical imaging method relying on the interference of a reference light with light back-scattered from the sample allowing depth-resolved measurements of back-scattering profiles (A-scans). Acquisition of many A-scans mapped over imaged region allows three-dimensional (3D) visualization of semitransparent tissue [1–3]. Over the last decade, OCT became a gold standard in the non-invasive evaluation of the volumetric structure of the eye. With the increasing need for improved sensitivity of ocular diagnostic, the simple measurements of tissue morphology, based on light scattering intensity alone, are often insufficient. Several functional extensions of OCT have been developed, over the years, that offer different contrast, greatly improving sensitivity for ocular diagnosis. These include Doppler OCT [4–7], OCTA [8], STdOCT [9], OCE [10,11], polarization-sensitive OCT [12], and spectroscopic OCT [13]. However, all these methods do not directly measure tissue functions but rather higher-order optical properties of tissue beyond simple scattering potential (blood flow, birefringence, absorption).

There are several functional objective tests of ocular tissue that extract light-evoked responses of the retina. These measurements include electroretinograms (ERGs) that measure light-evoked changes in the electrical activity of the retina. Following the same paradigm, there have been many successful attempts to observe light-evoked changes in the retina using optical imaging only [14–26]. This paper uses the term “optoretinogram” (ORG) to describe the light-evoked photoreceptor function measured by OCT. The term ORG was initially coined by Don MacLeod [27] and recently adopted by other investigators [28–32] for the paradigm of using NIR OCT to measure bleaching-induced back-scattering and/or elongation changes in the eye in vivo noninvasively. The name ORG draws an instructive parallel to the ERG, which has long been used to assess retinal function in vivo. Specifically, the ORG, like the ERG, comprises multiple components arising from distinct cells and mechanisms. Moreover, both ERG and ORG require an explanation of the underlying cellular and molecular mechanisms to achieve their full scientific and clinical utility. The most recent measurements of light-evoked intrinsic signals in photoreceptors can be divided into two categories, one extracting changes in retinal/photoreceptors reflectivity [33–36] and the other extracting changes in photoreceptors morphology. The latter splits into sub-categories: phase-based methods [17,26,28–31,37–41] and intensity-based methods [32,42,43]. Application of the phase-based ORG has been implemented mostly in cellular resolution AO-OCT systems that can visualize the photoreceptor mosaic allowing phase measurements of path length changes in outer segments of individual photoreceptors. All these methods required dedicated image acquisition systems (cellular resolution (3–4 um) or ultrahigh acquisition speed (MHz)), often not available in the current generation of clinical retinal imaging systems. However, the ultimate goal for this work should result in the implementation of light-evoked functional retinal tests using commercial OCT setups. In the OCT systems with standard clinical resolution (10–15 µm) and acquisition speed (∼100 kHz), the phase-based ORG is much more demanding to implement and explore because the systems cannot resolve individual photoreceptors and show instead speckle pattern corresponding to groups of cells [44], with no meaningful phase information due to temporal speckle decorrelation. However, a framework implemented to study phase-based light-evoked photoreceptor activity (ORG) with full-field swept-source optical coherence tomography (FF-SS-OCT) [29] proves that phase comparison between decorrelated speckle patterns is still possible by utilizing a series of images acquired before full decorrelation.

The main problem in implementing such a framework to raster scanning systems with standard acquisition speed is that the signal decorrelates very fast (tens of ms) due to saccadic and axial eye movements [41], which results in loss of phase-based ORG signals. Many groups have made a substantial effort to develop eye tracking and eye motion correction to reduce or remove tissue motion, allowing phase-sensitive measurements over seconds [45]. In intensity-based ORG measurements on mice, we observed that the light-evoked changes of the retina morphology continue over several seconds to minutes (an order of magnitude longer than in human studies). Therefore, phase-based analysis in mice faces challenges such as tissue movements caused by breathing, heartbeat, and eye drift, which are inevitable even after anesthesia.

In this paper, we describe the application of the modified Knox-Thompson method [46], where the cross-correlation of the phase difference between the measured layers of the object allows identification of the areas in which phase decorrelation occurs. In general, the method assembles the ORG signal from the most correlated elements of the inter-layer differential phase signal measured over time [47]. This work presents an evaluation of noise in phase-based length measurements in the sample consisting of a microscope coverslip and in a mouse retina in vivo. We also show that the ORG measurements, despite the relatively slow acquisition speed (equal to 100 kHz for A-scans and 100 Hz for B-scans), combined with the modified Knox-Thompson method, allow reconstruction of the phase-based ORG signal of the retina in an anesthetized mouse.

Several anatomical differences between mouse and human photoreceptors (both rods and cones) result in different appearances of the photoreceptor bands as measured by OCT. In human retinal imaging, cones reflectance provides a dominant signal and thus is naturally investigated for ORGs. On the other hand, human rods are much more challenging to image and analyze due to their lower reflectivity and size, but there are efforts to measure ORGs in rods [38]. There is no distinct difference between rods and cones reflectivity in mice, making the separation of two classes of photoreceptors very challenging. Additionally, due to the composition of the mouse retina (rod dominant), the measured ORG signal is believed to come mostly from rods. There are, however, efforts to measure cone ORGs in mice, and proposed phase-based analysis might provide access to the cone-based ORGs in mice [32].

The phase-sensitive ORG signal analysis presented here could be in principle extended to clinical retinal raster scanning OCT systems, potentially opening doors for clinically friendly ORG probing.

2. Methods

This paragraph will briefly describe details of our experimental system, including ORG data acquisition procedure and presentation of the main components of phased-based ORG processing.

2.1 Experiments in the living mice

All mouse husbandry and handling were in accord with protocols approved by the University of California Animal Care and Use Committee, which strictly adheres to all National Institutes of Health guidelines. Mice (B6-albinos (B6(Cg)-Tyrc-2J/J)) were obtained from Jackson Labs. A contact lens and gel (GelTeal Tears, Alcon) were used to maintain corneal transparency. All ORG experiments were performed in a darkroom. During the experiment’s mice were anesthetized with 2% isoflurane delivered in O₂, and maintained on a heated, adjustable platform adjacent to the imaging apparatus, as described previously [48].

2.2 Experimental system

For in vivo imaging of the mouse retina, we used a custom multimodal system with OCT and SLO sharing the same scanning system for simultaneous imaging of the retina. In the experiments here, only horizontal scanning is on without vertical scanning. Each horizontal scan collects 1024 pixels for SLO and 1024 A-scans for OCT at a 100 kHz rate. The SLO subsystem is equipped with a 488-nm laser (Coherent, 488-30FP) with an external trigger that enables precise control of light delivery during scanning. The SLO was used to deliver precisely calibrated 488-nm bleaching exposures. The SLO laser power was set to 4.7 µW at the mouse pupil. A single SLO scan (50 degrees) at this power bleaches ∼10% of the rhodopsin. The NIR Fourier domain-OCT system has a superluminescent diode light source (Superlum, T-860-HP) with an 82-nm bandwidth (effective fullwidth at half maximum [FWHM] after Hann windowing) centered at 860 nm and delivers 600 µW at the mouse pupil. The total 860 nm light exposure during an imaging session produced negligible bleaching or activation of rhodopsin. The lateral resolution of the system was characterized as ∼3 µm for SLO [49] and ∼5 µm for OCT, respectively.

The timing protocol for the experiment described in this paper is different from our standard protocol for in vivo imaging of the rod ORG, consisting of a series of scanning cycles (each taking about ∼0.92 seconds) with 4 B-scans collected in each cycle. Standard ORG protocol, used in our previous work relying on intensity-based light-evoked signal extraction, includes repeated raster scanning of the retina over field-of-view (FOV) of 50 × 25 degrees and a pixel sampling of 512 × 128 pixels at 100 kHz A-scan rate; the OCT signal is collected only from 4 centered adjacent B-scans in each raster scan. Considering the scanner fly-back time and the software resetting time, the time interval between two raster scans is about 0.92 s [32] which is ∼100 times slower than what we present here. Instead, the data is collected in a BM-scan mode, where consecutive B-scans are acquired over the same spatial location at a 100 Hz frame acquisition rate. For the “stimulus ON” experiment, the 488 nm laser was triggered to deliver the bleaching exposure for ∼1.4 s between B-scans 129 to 270. The scanning continued for an additional ∼40 s (4,000 B-scans) to record the ORG response. The “stimulus ON” scanning protocol was designed to analyze the phase stability between B-scans and ensure all the rhodopsin is bleached during light stimulation.

2.3 Initial OCT signal processing

The OCT interferometer acquires back-scattered light from the sample arm and allows its interference with the reference arm light. A line camera detects the spectrally modulated spectrum of light that undergoes the typical A-scan reconstruction steps: fixed pattern noise subtraction followed by wavenumber linearization, dispersion mismatch compensation, and spectral shaping. Finally, the Fourier transformation produces a complex-valued A-scan [50]. Several such complex A-scans acquired by moving the imaging beam with a fast x-scanner are combined to create a two-dimensional complex OCT B-scan:

2.4 Intensity-based processing of ORG signals in mice

The intensity-based processing of ORG signals in mice has been described in detail previously [32]. Here we just highlight its major components: the center 60%−80% (600–800 A-scans) portion of the BM-Stack was manually selected, then it was divided into sub-stacks, each with 20 A-scans widths, they were cross-correlated, axial position corrected and averaged, then, the averaged A-scans were interpolated axially to increase digital sampling to 0.1 µm/pixel and further cross-correlated, and axial potion corrected in the time direction to form a kymography flatten to BrM. Finally, the peak positions and peak intensities of BrM, ELM, ISOS, and other layers of interest were calculated using gaussian fit of individual layers to allow comparison with ORG extracted using phase-based processing. Note that the axial peak position estimation method is similar to the centroiding methods used in photoactivated localization microscopy (PALM), one of the super-resolution microscopy techniques [51]. Therefore, we can achieve sub-micron sensitivity of peak position shifts between individual OCT layers.

2.5 Phase-based processing of ORG signals

The phase-based processing of ORG signals in mice has been implemented to improve the sensitivity or ORG measurements (change in thickness of photoreceptor structures). In general, we expect phase-based position detection methods to be more sensitive in OCT processing than intensity-based position detection methods. Figure 1. shows an outline of phase-based processing of ORG signals with visualization of intermediate results. In the first step to extract phase information from layers of interest, an intensity-based OCT segmentation is performed by semi-automatic procedure derived from the work of Yin [52]. Initially, manual segmentation is performed on one specific Intensity B-scan to define the layer of interest. For all other B-scans, segmentation is performed automatically by an iterative Monte-Carlo method minimizing a cost function based on intensity changes between the line and its envelope (so-called double-line algorithm). After this intensity segmentation, we obtain reference positions for the phase processing. We refer to them as lay1 and lay2 for two retinal layers, respectively. To observe phase changes between two segmented layers in time t, we adopted the method of cross-spectra analysis first proposed by Spahr et al. [47]. Cross-spectrum is a two-dimensional representation of time-dependent phase-correlation between complex numbers. We calculate the cross-spectrum of the time evolution of a point with coordinates x and z, as:

 

Fig. 1. Outline of the processing steps for phase-based ORG signal extraction, with visualization of intermediate results. (a) OCT BM-Scans. (b) Time evolution of the cross-spectrum phase from single point in the retina (lay), representative position of point is marked by orange rectangle in (a). (c) Time evolution of a phase difference of the cross spectrum between two spatially separated points from two layers (lay1 & lay2), representative position of points is marked by two green rectangles in (a). (d) Time evolution of a phase difference of the averaged cross spectrum, over window od size M, between two spatially separated layers (lay1 & lay2), representative position of windows is marked by two blue rectangles in (a). (e) Weights w showing phase correlation between points within the window M. (f) Binary-valued weights b, calculated to extract well-defined regions of phase cross-spectrum. (g) Reduced weighted time evolution of a phase difference of the averaged cross spectrum from (d). (h) Visualization of phase extraction paths from g, brown arrow (along the first (left) vertical line) and green arrows (using the Knox-Thompson paths). (i) Time-dependent phase difference between two layers extracted using single path. (j) Time-dependent phase difference between two layers extracted using the Knox-Thompson paths. See Eq. (1–9) in text for details regarding mathematical operations.

Download Full Size | PDF

Next, a phase difference between two spatially separated points (in our case, two layers of interest separated by a distance along the z-axis on an A-scan) is calculated by multiplying the cross-spectrum of the first layer with conjugation of cross-spectrum of the second layer, see Fig. 1(c):

The difference signal between two layers U_lay1/lay2 is averaged in space (over multiple A-scans) over window od size M, to improve the accuracy of phase difference estimation for ORG, see Fig. 1(d):

The weights w can be interpreted as the phase correlation between points within the window M. The weights w has the value of 0 for totally random phase shifts and 1 for fully coherent phase shifts. Next, we calculate w² for easier visualization of phase correlations, and we use it for thresholding of the $U_{lay1/lay2}^M$. The threshold was selected experimentally on the level th = 0.02. In our experiments, this thresholding removed ∼ 55% of the most decorrelated pairs of B-scans (separated by Δt) in mouse data and ∼ 14% of the final phase ORG signal. A real-valued weights b, shown in Fig. 1(f), used to extract well-defined regions of phase cross-spectrum are defined below:

Providing the reduced weighted cross spectra between layers 1 and 2 over window M, see examples in Fig. 1(g), Fig. 7(c), Fig. 9(c), Fig. 10(c), and Fig. 11(a,d,g,j):

The regions where $U_{lay1/lay2}^{bM}$ was reduced to zero (created by zeros in the weights b (Eq. (6)) were interpolated using an averaged value from ten nearest points (+/- 5 $\Delta t$), resulting in weighted corrected cross-spectra $U_{lay1/lay2}^{cM}$. This worked well for the $U_{lay1/lay2}^{cM}$ values near the base ($\Delta t$ = 0), which is sufficient for Knox-Thompson reconstruction.

To extract time-dependent phase difference between two layers using weighted corrected cross spectra, which is needed for ORG signal extraction, one needs to define the reference time t₀. Once this is done, the ORG signal can be reconstructed using:

In order to improve ORG signal fidelity over larger time intervals, the final phase signal $\Delta \varphi _{lay1/lay2}^{KT}$ was computed using the Knox-Thompson paths [46,47]. The number of paths has been selected in such a way that allowed phase wrapping and decorrelated region (Fig. 1(h)). Signal reconstruction was done by assembling Knox-Thompson paths of length p between B-scans indicated with arrows. The total number of Knox-Thompson paths Π is given by: $\Pi = {{\rm Z}_{quotient}}({\Delta {t_{\max }}/(p \cdot \Delta {t_{\min }})} )$, where Δt_max is the total measurement time. Current Knox-Thompson path number δ for time difference Δt is defined as $\delta = {{\rm Z}_{quotient}}({\Delta t/(p \cdot \Delta {t_{\min }})} )$ where Δt_min is the time step between two consecutive B-scans. Then we define $\Delta \varphi _{lay1/lay2}^{KT}$ as a phase difference at the end of δ-th of Knox-Thompson paths, see Fig. 1(j):

Binarization function b was computed the similar way as $\Delta \varphi _{lay1/lay2}^{KT}$ by path collection from Π Knox-Thompson paths:

The phase differences with weights are then defined as:

Finally, the time-depended layer displacement can be calculated from:

2.6 Phase-based distance measurements error

The theoretical sensitivity of phase difference measurements depends on two components. The first component is defined by the signal-to-noise ratio of OCT data [53,54]. In the case of the ORG measurements, since we compare phases between two layers, SNR of each layer needs to be taken into account:

The second component is defined by a reduction in coherence between two phase measurements. In classical Doppler measurements, the reduction of coherence is simulated by phase error caused by lateral displacement, Δx between two consecutive phase measurements acquired with the beam of diameter d [53,54]:

This results in the total phase sensitivity estimation of:

Note that from Eq. (14), the maximum value of Doppler phase error calculated between two totally decorrelated measurements (d >> Δx) should be equal to $\sqrt {4\pi /3}$ rad. Considering Eq. (12) and our central wavelength of 860 nm, this translates to σ_Δx,max, ∼ 105 nm. In the ORG measurements, the phase difference is calculated between two points separated in space, not in time. This means that the phase variance of the ORG signal calculated between two positions should depend on how correlated the two fields are. Resulting in the experimental values consistently higher than theoretical phase variance estimation σ_{Δph, ORG} for correlated fields:

In the phase estimation presented here, two random independent speckle fields are used to calculate phase difference. Thus, we should expect to have the phase error σ_incoherence in the range of σ_Δx,max.

Another consequence of this phase decorrelation between two layers is that ORG signal averaging between different A-scans will not follow simple SNR-dependent improvement of phase sensitivity (with the square root of the number of averaged signals M). Therefore, we decided to use a more elaborate Kernel to investigate the effect of averaging on ORG signal variance. We found experimentally that the following equation describes the improvement of differential phase error, due to A-scan averaging, very accurately (see Fig. 3):

3. Results

3.1 System stability characterization

The system was characterized by signal-to-noise ratio (SNR) and temporal phase stability. A measure of phase sensitivity or stability is given by the standard deviation of the measured phase differences. Figure 2 showcases the effect of the galvo-scanner jitter quantified by measurements of phase sensitivity under three conditions: (cond. I) galvo-scanner powered off, 4000 BM-scans of a coverslip glass plate, (cond. II) galvo-scanner powered on, and linear scan of a coverslip glass plate, 4000 BM-scans, and (cond. III) a galvo-scanner powered on and linear scan of the living mouse retina. The first column of the Fig. 2 shows an example of OCT intensity B-scans for three experiments. The following two columns show the time-varying phase extracted from the OCT data over the whole duration of the experiment (40 s) and the zoomed-in portion (2 s). The fourth column shows temporal Phase Power Spectrum (FFT of phase signals from layers 1–2 and the phase difference between layers), and the fifth column shows a phase histogram for corresponding data sets. For absolute phase variation, the phase on the front surface of coverslip over time was obtained for conditions (cond. I) and (cond. II) in a BM-scan. For (cond. III), the phase in the same spatial location on a layer in mouse retina in each BM-scan was obtained along the time. Note, that cond. I and cond. II represent static sample, while cond. III represent moving sample (living mouse retina). The theoretical sensitivity σ_the (calculated using SNR) of the raw phase for one A-scan was equal to 0.04, 0.19, and 0.24 radians respectively for the three conditions (cond. I), (cond. II) and (cond. III). The standard deviation of the measured raw unreferenced phase over time was equal to 0.20, 1.03, and 1.30 radians, respectively for the three conditions (cond. I), (cond. II), and (cond. III), confirming deviation from shot-noise limited detection of the phase.

 

Fig. 2. Signal analysis of the temporal phase stability. The first column presents an example of one OCT B-scan, from series of BM-scans, for glass plate (a, f) and mouse retina (k). The second column presents phase stability for edges of the coverslip glass in the sample channel with the scanners off (b) and on (g), and for layers in a scanned mouse retina (l). Red and blue colors correspond to phase calculated from layers 1 & 2 marked by rectangles with corresponding colors in (a, f, k). The black line shows the phase difference between the signal from layers 1 and 2. The third column presents the magnified phase results from (b, g, l) - seconds 10 to 12 marked by a dashed green line. The fourth column shows temporal Phase Power Spectrum (FFT of phase signals from layers 1–2 and the phase difference between layers). The fifth column shows a phase histogram for corresponding data sets.

Download Full Size | PDF

In order to showcase that our phase stability was not negatively affected by the time interval between phase measurements and to show that our system is phase stable, we compared our previously measured temporal phase stability between B-scans to temporal phase stability between A-scans. This resulted in two orders of magnitude reduction in the time interval between phase measurements. Figure 3 shows an analysis of differential phase stability in this case for the microscope coverslip and both galvo scanners powered off (cond. I). Note that Fig. 3(c) shows virtually the same differential phase histogram as shown previously for longer time intervals in Fig. 2(e,j) for conditions (cond. I) and (cond. II).

 

Fig. 3. Signal analysis of temporal phase stability for 100 kHz A-scan acquisition (cond. I). (a) Phase stability for edges of the coverslip glass in the sample channel with the scanners off. Red and blue colors correspond to phase calculated from layers 1 & 2 marked by rectangles with corresponding colors in Fig. 1(a). The black line shows the phase difference between the signal from layers 1 and 2. (b) Temporal Phase Power Spectrum (FFT of phase signals from layers 1–2 and for the phase difference between layers). c) Histogram for phase difference.

Download Full Size | PDF

3.2 Effect of averaging on phase-based processing

To estimate the effect of averaging on differential phase measurement error, we evaluated data acquired with cond. I) galvo-scanner powered off, acquired at coverslip glass plate, and cond. III) a galvo powered and linear galvo scan of the mouse retina. The changes in phase stability as a function of number averaged A-scans were measured and fitted using Eq. (17). Figure 4 summarizes these results for both theoretical (SNR based: σ_{Δph, ORG}) and experimental (direct: σ_{phase, ORG}) differential phase error estimation.

 

Fig. 4. Estimation of phase noise theoretically – SNR-based (red line) and experimentally (black line) for a glass plate with scan off (a) and mouse retina with scan on (b). Coefficient of determination for (a) and (b) is equal R² = 0.99.

Download Full Size | PDF

The results shown in Fig. 4 allow us to estimate the phase error σ_incoherence using Eq. (16), due to decorrelation between two layers. This error is equal to ∼9.6 nm and ∼102 nm for glass plate and mouse retina differential phase measurements, respectively. This result is expected since reflections of two surfaces of a glass plate should be well correlated while scattering between two retinal layers should be decorrelated. Another interesting observation is that averaging of A-scans allows a very efficient reduction of the differential phase error bringing the theoretical limit of improvement by averaging to the same range (see from (Eq. (17)) c_exp and c_the for both measurements). This also showcases the need for averaging of differential phase measurements between A-scan when extracting light-evoked changes in layer thickness. The similar values of fitting parameters a_the/exp, and b_the/exp (both in SNR-based estimation and direct phase measurements) suggest that degree of differential phase decorrelation does not affect the averaging process.

3.3 Measurements of phase-based ORG signal

In our ORG measurements, the light-evoked changes in retinal activity are measured as the relative swelling of the outer retina, which is extracted by measuring a change in relative distance between the BrM and ELM directly proportional to the phase difference (Eq. (12)). However, a simple calculation of the phase change between the first and subsequent B-scans is not feasible due to the OCT signal decorrelation (within the same layer) over time. This signal decorrelation (due to the eye movements and tissue internal dynamics) causes the measured phase differential signal to disappear, as shown in Visualization 1.avi and Visualization 2.avi. Selected frames are shown in Fig. 5 and Fig. 8, where we can see that for a time difference of Δt = 0.01 s, the phase difference between two B-scans is well defined. After 11 seconds, the common phase information between the same structures in the 1^st and 1060^th B-scan has decreased, manifested by a low phase difference signal (phase difference becomes random). This is why the use of periodic signal self-referencing is necessary to allow the comparison of successive phase differences over a longer time.

3.3.1 Control experiment (no light stimulation)

ORG control experiment measured changes in the relative swelling of the outer retina (distance between the BrM and ELM) without any visible light stimulation. The only source of bleaching is the imaging OCT light itself. Figure 5 showcases the phases stability during the experiment performed over 40 s acquisition window.

The problem of time-dependent signal decorrelation observed in Fig. 5. can be solved using phase analysis employing Knox-Thompson paths [46]. The light-induced retinal activity in the mouse retina was mapped as the time-dependent change in outer retina thickness (between ELM and BrM). These reference layers were presented in Fig. 2, Fig. 5, and Fig. 8. For comparison, we measured the light-evoked activity using our previously developed intensity-based methods and newly developed phase-based method, as described in section §2.7. In order to allow phase-based ORG signal extraction, we calculated cross-spectra differences $U_{lay1/lay2}^{M = 150}$ from 150 A-scans (from 350^th A-scan to 500^th A-scan). The number of the A-scans that need to be averaged will depend on the OCT signal intensity and phase stability of the optical system. Figure 4 shows the change of the phase noise level depending on the number of averaged A-scans. According to Eq. (17), we can estimate the theoretical phase error depending on the number of averaged A-scans, corresponding to the accuracy of the phase-based ORG signal. We have chosen the number of A-scans to be averaged based on the difference between the experimental standard deviation of phase and the theoretical phase error levels. We found that the number of A-scans around 150 (see c_exp and c_the in Fig. 4(b)) provide the smallest difference. Figure 6(a) and Fig. 9(a) show examples of the cross-spectrum before averaging. Based on this, we conclude that solid ORG responses could be traced with a low number of averaged A-scans, while week ORG responses would require a large number of averaged A-scans.

 

Fig. 5. Phase stability during measurement of a mouse retina. (a & e) The OCT intensity in log₁₀ scale of B-scan #2 and #1060 was collected in time 0.01 s and 10.59 s from the first reference B-scan. The red and blue rectangles indicate ELM and BrM layers accordingly. The yellow dashed lines presented a range of A-scans that were taken to analyze with Knox-Thompson paths. (b & f) Reference phase of first B-scan in radians. The color scale bars (b-d & f-h) indicate values of phase from –π to π from blue to red, respectively. (c & g) The phase of B-scans collected in time 0.01 s and 10.59 s from the first reference B-scan, respectively. (d & h) The phase difference between phases of first reference B-scan and B-scan collected after time interval Δt. The phase stability over time is presented in Visualization 1.avi.

Download Full Size | PDF

 

Fig. 6. The results of the control experiment (no visible light stimulation) in the mouse retina with a distance between ELM and BrM layers followed over 40 s. (a) The phase of $U_{lay1/lay2}^{M = 1}$ cross-spectra layer differences for one A-scan. (b) Result of averaging M = 150 cross-spectra differences $U_{lay1/lay2}^{M = 150}$ for ELM and BrM layers. A-scan range taken to average is shown in Fig. 4(a,c). (c) The binary map b was used to threshold and remove less correlated phasors. The white region indicates a well-defined phase signal that could be taken to analyze. (d) Comparison of Knox-Thompson signal without thresholding (red line) and with weights (black line). (e) Comparison ORG phase-based signal with ORG intensity-based signal. The ORG phase signal was reconstructed from Π=5 Knox-Thompson paths. (f) Power spectrum (FFT of phase signals from ELM and BrM layers).

Download Full Size | PDF

The range of A-scans used to extract the ORG signal is indicated in Fig. 5(a). As a control experiment, the results of phase-based ORG signal extracted from 40 sec. acquisition without light stimulation are presented in Fig. 6. The phase difference signal from a single A-scan was characterized by a standard deviation σ_lin. The σ_lin varied from A-scan to A-scan in the range of 1.0–2.1 rad (51–107 nm). This variation in σ_lin is probably due to the local tissue properties of ELM and BrM and represents a varying contribution of back-scattering and back-reflection from these structures. Low σ_lin corresponds to reflection dominant signal, and high σ_lin corresponds to scattering dominant signal (see a change in ELM and BrM brightens depending on the lateral position of A-scans in Fig. 5(a,e). This is why spatial signal averaging was needed to reduce the phase variance, as presented in Fig. 6(b). The spatial averaging over M = 150 A-scans reduced the standard deviation of phase difference to σ_lin=0.19–0.39 rad (10–20 nm), depending on the specific A-scans used for averaging. The initially constant phase difference (before light stimulation) is associated with the baseline “dark adapted” state of the retina and can be visible in the averaged image Fig. 6(d), better depicted in the portion of the experiment acquired during the first 10 s Fig. 7(d).

 

Fig. 7. Magnified results of the distance between ELM and BrM layers in the first 10 s without light stimulus. (a) The phase of $U_{lay1/lay2}^{M = 1}$ cross-spectra layer differences for one A-scan. (b) Result of averaging M = 150 cross-spectra differences $U_{lay1/lay2}^{M = 150}$ for ELM and BrM layers. A-scan range taken to average is shown in Fig. 4(a). (c) The averaged $U_{lay1/lay2}^{b,M = 150}$ signal after thresholding to remove less correlated phasors. The yellow rectangles in (b-c) indicate a region where we can see decorrelated signal reduction using weights. (d) Comparison of Knox-Thompson signal without thresholding (red line) and with weights (black line). (e) Comparison ORG phase-based signal with ORG intensity-based signal. The ORG phase signal was reconstructed from Π=10 Knox-Thompson paths. (f) The power spectrum for phase difference signal (FFT of time-varying phase signals from ELM and BrM layers).

Download Full Size | PDF

The phase ORG signal in Fig. 6 was reconstructed using Π = 5 Knox-Thompson paths corresponding to Δt_KT= 8 s. The phase ORG signal in Fig. 7 was reconstructed using Π=10 Knox-Thompson paths corresponding to Δt_KT= 4 s. The thresholding on $\Delta \varphi _{lay1/lay2}^{KT}$ signal for Δt_KT= 8 s removed 3% of the least correlated phase and 26% of data from a total phase of $U_{lay1/lay2}^{M = 150}$. The thresholding on zoomed $\Delta \varphi _{lay1/lay2}^{KT}$ signal for Δt_KT=4 s removed 2% of the least correlated phase and 5% of data from a total phase of $U_{lay1/lay2}^{M = 150}$. The phase-based ORG results are compared with averaged intensity-based ORG calculations for the same data set in Fig. 6(d,e) and Fig. 7(d,e) for shorter time evaluation. For comparison, the standard deviation for the phase-based ORG method, the variance of position difference between ELM and BrM intensity maxima, was σ_lin ∼ 12 nm (Fig. 7(e)). The standard deviation for the intensity-based ORG method, the variance of position difference between ELM and BrM intensity maxima, was σ_lin=45 nm (∼4 times higher than for phase-based ORG method).

The increase in distance between ELM and BrM layers in this experiment represents the unwanted activation of photoreceptors by the NIR OCT imaging beam bleaching 0.01% of rhodopsin during the whole imaging session (∼12000 B-scans (∼2 min) in total, including the alignment, the pre-imaging scans, and the imaging acquisition). The method to calculate the bleach levels in mice was described in detail previously [55], and includes formulas for pigmented and albino mice. Briefly, the spectral power distribution function (W/nm) of the OCT light at the pupil was determined from the spectrometer, and the total power was measured by Thorlabs S120C photodiode with its wavelength setting at 860 nm. For calculation of bleaching, the spectral energy density was converted into rhodopsin-equivalent units, using a Lamb–Govardovskii [56,57], template for mouse rhodopsin, assuming the axial density of rhodopsin in mouse rods to be 0.35. Note that total elongation of ELM and BrM due to this continuous stimulation equaled 100 nm over 40 s and that it can be only seen on phase-based ORG analysis. Also, note that in this paper, we use serial B-scans for fast ORG signal measurement rather than serial volumes used in our standard intensity-based ORG signal measurements [32,43], when local bleaching by OCT light is 100x lower, and thus its effects are negligible. The periodic pattern seen on magnified, and A-scan averaged cross-spectra in Fig. 7(c) is depicted by the power frequency analysis seen on Fig. 7(f) and corresponds to the modulations of 1.3 Hz and 9.5 Hz. The potential origins of these modulations are breathing and heartbeat, respectively. Their effect on phase-based ORG signal variance will be briefly discussed in the discussion and conclusions section.

 

Fig. 8. Phase stability during measurement of the light-evoked responses of the mouse retina. (a & e) The OCT intensity in log₁₀ scale of B-scan #2 and #1060 was collected in time 0.01 s and 10.59 s from the first reference B-scan. The red and blue rectangles indicate ELM and BrM layers accordingly. The yellow dashed lines presented a range of A-scans that were taken to analyze with Knox-Thompson paths. (b & f) Reference phase of first B-scan in radians. The color scale bars (b-d & f-h) indicate values of phase from –π to π from blue to red, respectively. (c & g) The phase of B-scans collected in time 0.01 s and 10.59 s from the first reference B-scan, respectively. (d & h) The phase difference between phases of first reference B-scan and B-scan collected after time interval Δt. The phase stability over time is presented in Visualization 2.avi.

Download Full Size | PDF

3.3.2 Light evoked ORG experiment

After our baseline experiment (no visible light stimulation) present in the previous section, here we want to present our phase-based ORG analysis of ELM-BrM distance changes due to visible light stimulation. The bleaching light was introduced after 1.29 s from the start of the OCT acquisition. The stimulus continued for ∼1.41 s. Similar to previously observed time-dependent signal decorrelation depicted in Fig. 5 in our baseline experiment, here as well the same effect is evident as showcased in Fig. 8. To resolve this problem, we used the same phase analysis employing Knox-Thompson paths and mapping light-induced retinal activity as the time-dependent change in outer retina thickness (between ELM and BrM). Here as well the phase-based extraction of ORG signals was compared with the intensity-based method. To allow phase-based ORG signal extraction, we calculated cross-spectra differences $U_{lay1/lay2}^{M = 150}$ from 150 A-scans (from 350^th A-scan to 500^th A-scan). The range of A-scans used for our analysis is indicated in Fig. 8(a).

 

Fig. 9. The light-evoked activity of the mouse retina for ELM and BrM layers. (a) The phase of $U_{lay1/lay2}^{M = 1}$ cross-spectra layer differences for one A-scan with light stimulus. (b) Result of averaging M = 150 cross-spectra differences $U_{lay1/lay2}^{M = 150}$ for ELM and BrM layers. A-scan range taken to average is shown in Fig. 8(a). (c) The averaged $U_{lay1/lay2}^{b,M = 150}$ signal after thresholding to remove less correlated phasors. (d) Comparison of Knox-Thompson signal without thresholding (red line) and with weights (black line). (e) Comparison ORG phase-based signal with ORG intensity-based signal. The ORG phase signal was reconstructed from Π=20 Knox-Thompson paths. (f) The power spectrum for phase difference signal (FFT of phase signals from ELM and BrM layers).

Download Full Size | PDF

The results of the phase-based ORG signal extracted in the experiment with light stimulation are presented in Fig. 9 for the whole duration of the experiment. Here the ORG was reconstructed using Π = 20 Knox-Thompson paths corresponding to Δt_KT= 2 s. The phase ORG signal in Fig. 10 was reconstructed using Π = 10 Knox-Thompson paths corresponding to Δt_KT= 4 s. The thresholding on $\Delta \varphi _{lay1/lay2}^{KT}$ signal Δt_KT= 2 s removed 15% of the least correlated phase and 58% of data from a total phase of $U_{lay1/lay2}^{M = 150}$. The thresholding on zoomed $\Delta \varphi _{lay1/lay2}^{KT}$ signal for Δt_KT= 4 s removed 10% of the least correlated phase and 29% of data from a total phase of $U_{lay1/lay2}^{M = 150}$. Figure 9(e) and Fig. 10(e) show the comparison of the ORG signals extracted with phase-based and intensity-based calculations. Here as well, the phase difference signal from a single A-scan was characterized by a standard deviation σ_lin. The σ_lin varied from A-scan to A-scan in the range of 1.0-2.4 rad (51–123 nm). After spatial signal averaging, the phase variance was reduced, as presented in Fig. 9(b) and Fig. 10(b). The spatial averaging over M = 150 A-scans reduced the standard deviation of phase difference to σ_lin = 0.29–0.67 rad (15–34 nm), depending on the specific A-scans used for averaging. Note that the constant slope was removed to estimate the standard deviation of phase difference during light stimulation. The initial brief constant phase difference (before visible light stimulation) can be visible in the averaged image Fig. 9(d,e), better depicted in the magnified portion of the experiment acquired during the first 10 s Fig. 10(d,e). The standard deviation for phase-based ORG method, the variance of position difference between ELM and BrM intensity maxima, was σ_lin ∼ 15 nm (Fig. 10(e)). The standard deviation for the intensity-based ORG method, the variance of position difference between ELM and BrM intensity maxima, was σ_lin = 68 nm, ∼ 4 times higher than for the phase-based ORG method.

 

Fig. 10. Magnified results of the light-evoked activity of the mouse retina for ELM and BrM layers in the first 10 s in response to 1.41 s of light stimulation. (a) The phase of $U_{lay1/lay2}^{M = 1}$ cross-spectra layer differences for one A-scan with light stimulus. (b) Result of averaging M = 150 cross-spectra differences $U_{lay1/lay2}^{M = 150}$ for ELM and BrM layers. A-scan range taken to average is shown in Fig. 7(a). (c) The averaged $U_{lay1/lay2}^{b,M = 150}$ signal after thresholding to remove less correlated phasors. The yellow rectangles in (b-c) indicate a region where we can see decorrelated signal reduction using weights. (d) Comparison of Knox-Thompson signal without thresholding and (red line) with weights (black line). (e) Comparison ORG phase-based signal with ORG intensity-based signal. The ORG phase signal was reconstructed from Π=10 Knox-Thompson paths. The blue rectangle indicates light exposure time. (f) The power spectrum for phase difference signal (FFT of phase signals from ELM and BrM layers).

Download Full Size | PDF

The increase in distance between ELM and BrM layers in this experiment is much more evident than in our baseline studies due to the estimated total bleaching (100%) of photoreceptors rhodopsin during 1.41 s visible light exposure. Total elongation of ELM and BrM due to this stimulation equaled 1.5 µm over 40 s (15 times more than in our baseline experiment), and that it can be easily seen on both phase-based and intensity-based ORG analysis. Note that here as well, the periodic pattern can be seen on zoomed and A-scan averaged cross spectra in Fig. 10(c), and that it is clearly depicted by the power frequency analysis seen on Fig. 10(f) and corresponds to the modulations of 1.3 Hz and 10.4 Hz, similar to the values observed in control experiment presented in section §3.3.1 and assigned to breathing and heartbeat, respectively.

3.3.3 Repeatability of phase-based ORG signal extraction

To showcase the repeatability of the phase-base ORG signal extraction, Fig. 11 shows the results of four more representative experiments (two without light stimulation and two with light stimulation).

 

Fig. 11. The results of additional control experiment (Experiment 1 & 2) and light-evoked activity experiment (Experiment 3 & 4) of the mouse retina showing the distance between ELM and BrM layers in the first 10 s in response to 0.85 s (Experiment 3) and 0.62 s (Experiment 4) light stimulation. (a, d, g, j) The averaged $U_{lay1/lay2}^{b,M = 150}$ signals after thresholding to remove less correlated phasors are shown for each experiment. (b, e, h, k) Comparison of Knox-Thompson signal without thresholding (red line) and with weights (black line) for each experiment, respectively. (c, f, i, l) Comparison between ORG phase-based signal with ORG intensity-based signal. The ORG phase-based signal was reconstructed from Π=4 Knox-Thompson paths and threshold b = 0.03. The blue rectangle indicates blue light exposure time.

Download Full Size | PDF

These results confirm the improved accuracy of phased-based ORG signal extraction compared to our standard intensity-based ORG analysis.

Note that decorrelation between B-scans can be estimated by evaluating the averaged $U_{lay1/lay2}^{b,M = 150}$ signals after thresholding. In all experiments, phase variations corresponding to breathing and heartbeat can be easily identified, further confirming similar phase sensitivity with previous data. As an example, “Experiment 2” data shows a breathing artifacts-free dataset. In all experiments, phase-based ORG signal extraction resulted in much less noisy estimation of the ORG signals, see Fig. 11(c, f, i, l). Application of the Knox-Thompson paths allowed reliable extraction of the ORG signal even in Experiment 4 (Fig. 11(j- l)), where the OCT signal becomes decorrelated after about 5 seconds from the start of the acquisition.

4. Discussion and conclusions

OCT is a coherent detection technique that carries not only information about scattering intensity but also about the axial position of the scattering object encoded in both the peak position of the signal and in a relative phase of its interferometric fringe. The FFT amplitude of the fringes allows for obtaining information about the structure of the object. Evaluation of the FFT phase of spectral fringes allows access to information about axial shifts smaller than the axial resolution of the OCT system. The work presented here described the practical implementation of the framework to use phase-based OCT signal analysis to measure axial displacement between two randomly scattering structures generating time-varying, uncorrelated speckle fields.

Lack of melanin in albino mice allows clear visualization of OCT signal back-scattered from ROS tips (ROST) and from BrM, therefore only B6-albinos (B6(Cg)-Tyrc-2J/J) mice, which are congenic with B6 mice (C57BL/6(J)) save for the null mutation in the tyrosinase gene that eliminates melanin, were used for this study. However, the methodology developed here could be applied to pigmented mice as well. We intend to work in the future on extending our phase-based ORG method analysis to mouse models with varying pigmentation levels [58].

Our previous work [43] has shown that intensity-based processing can be applied to measure the light-evoked response of the retina in an anesthetized mouse. These methods allow the extraction of time-dependent changes in the thickness of the retinal layers using simple and fast data processing. Additionally, intensity-based processing was relatively not sensitive to the eye motion artifacts. Alternatively, in OCT data, changes in the relative position of retinal layers can be measured using phase-based processing [39]. The phase-based position detection in interferometric methods is in general more sensitive than intensity-based peak position detection. However, phase-based measurement is more sensitive to bulk motion artifacts, which makes it particularly vulnerable to the eye motion artifacts mainly caused by saccades, heartbeat, and breath. One way to avoid these artifacts is to increase the acquisition speed. Publications [17,40,41] show that with an A-scan acquisition rate of 1 MHz and motion correction, it is possible to measure the phase shifts within individual photoreceptors caused by swelling of photoreceptors outer segments using direct phase analysis methods using flying spot AO-OCT datasets. The standard approach is to subtract consecutive A-scans from the reference A-scan, then measure the signal change over time.

However, as shown in this paper, during the measurement lasting from a few to several dozens of seconds, and in the presence of speckles, there are time varying changes in signal correlation. Here, we have visualized the mutual correlation of each B-scan by implementation of the cross-spectra analysis allowing for easy separation of highly correlated regions from weakly correlated ones. Using the Knox-Thompson method with the path analysis (folding arrows in Fig. 1(h,j)), it is possible to perform the phase analysis of the retina during light stimulation for systems with an A-scan acquisition rate in the range of 100 kHz and multiple scattering dominant (presence of speckle) of OCT signal. The threshold providing the reduced weighted cross spectra between layers 1 and 2 over window M was set so that to remove artifacts related primarily to the breath of the mice. The breath of the mice can be observed by phase jumps, as shown in Fig. 7(d) and Fig. 10(d) and Fig. 11(b, h, k). The red line in Fig. 7(d) and Fig. 10(d) show the ORG signal without thresholding compared to the final threshold signal drawn in the black line. Overestimating the threshold value could include the decorrelated signal of the mouse breathing in ORG calculations. The standard deviation of the measured ORG signal might increase, and the accuracy of the calculation decreases. On the other hand, the threshold should be chosen to retain useful phase information. Underestimating the threshold might demote the precision of the ORG measurement.

We have shown that for a 100 Hz serial B-scan acquisition, and in the case of an immobilized animal (mouse), we are able to reconstruct the phase-based ORG signal from the BrM and ELM layers using the Knox-Thompson path analysis. We achieved compliance of the phase-based ORG signal with the amplitude-based ORG signal. As expected, the phase signal shows a higher sensitivity depending on the number of averaged A-scans and signal SNR. As the sensitivity of the phase method is higher, we could observe a decrease in the ELM-BrM distance with the onset of the light stimulus (Fig. 10(d)), plausibly caused by the contraction of the photoreceptor OS. This effect was not visible in the amplitude measurement. The relatively high value of phase difference variance in the in vivo ORG experiment is due to residual phase shifts in the retina. A closer look at the power spectrum of the phase-based ORG signals shown in Fig. 6(f) and Fig. 7(f) for control and Figs. 9(f) and 10(f) for light stimulus case reveal the presence of slow (∼1.3 Hz for both experiments) and fast (∼9.5 Hz and ∼10.4 Hz for control and light stimulus experiments respectively) phase modulations. We have confirmed that this corresponds to mouse breading and heart bit frequency, respectively. This showcase that our phase-based ORG measurements are more sensitive than our current intensity-based distance analysis, offering about 10 times improvement in sensitivity of phase-based ORG estimation.

Since the distance between the ELM-BrM boundary rather than individual compartments of this band (RPE, SRS, OS, IS) were measured in the present study, we cannot directly attribute our observations with photoreceptor swelling. One alternative possibility is that ELM-BrM distance change could correspond to light-induced subretinal space (SRS) expansion caused by mitochondrial respiration [59], where SRS is defined as the extracellular fluid space between RPE cells and photoreceptors, with the growing body of work researching this topic [60–63]. We have, however, previously studied the light-evoked SRS changes [32], and did not find any SRS expansion after a visible light stimulus. One noticeable difference between these experiments is the type of anesthesia used in the study. Other groups often use a mix of ketamine (100 mg/kg) and xylazine (6 mg/kg) as injectable anesthetics, while our mice are always anesthetized with 1.5–2% isoflurane mixed with oxygen. Taking into account how different methods of anesthesia affects mouse physiology, such as blood flow [64], intraocular pressure (IOP) [65], or eye movements and electroretinogram (ERG) [66], it is possible that it could affect the light-evoked response of the retina in general and of the SRS in particular, despite some reports suggesting no difference [61]. Therefore, future ORG studies in animal models should include evaluating the different effects anesthesia could have on the extracted signals.

Over the last few years there has been rapid progress in development of the non-invasive, optical measurement of neural function in the retina. These functions have been investigated using imaging systems already used in the clinic or intended for clinical and translational use. While the most success in ORG signal extraction has been achieved with a phase-based analysis of data acquired with optical coherence tomography (OCT) with adaptive optics (AO), its broader application in clinical settings remains limited mostly due to high cost of AO-OCT systems and difficulty of its use in truly clinical settings. The work presented here contributes to developing sensitive phase-based methods to extract ORG signal from clinical-grade non-AO OCT systems. The ORG methodology has demonstrated the ability to detect responses to stimuli isomerizing less than 0.01% of photopigment. Therefore it has the potential to prove a quick, non-invasive, and objective way to measure subtle disease-related dysfunction and thus to provide an entirely new and complementary biomarker for retinal disease and recovery [67]. Using the presented cross-correlation spectrum analysis framework, the use of phase-based ORG signal extraction should be possible in the current clinical systems operating at B-scan acquisition rates of 0.1 Hz and without active eye-tracking systems or AO correction. Presented data processing was performed after the acquisition of serial B-scans, and lateral and axial displacement of tissue was waited by cross-spectra decorrelation and Knox-Thompson paths calculations. Thus, phase errors caused by large transverse or axial displacements were removed from the analysis. This development of a phase-sensitive analysis method that is insensitive to speckle signal decorrelation due to the eye movements represents an important step towards implementing non-invasive procedures of measuring the light-evoked retinal function in the human clinic.