**Significance:** There are no label-free imaging descriptors related to physiological activity of inner retinal cells in the living human eye. A major reason is that inner retinal neurons are highly transparent and reflect little light, making them extremely difficult to visualize and quantify.

**Aim:** To measure physiologically-induced optical changes of inner retinal cells despite their challenging optical properties.

**Approach:** We developed an imaging method based on adaptive optics and optical coherence tomography (AO-OCT) and a suite of postprocessing algorithms, most notably a new temporal correlation method.

**Results:** We captured the temporal dynamics of entire inner retinal layers, of specific tissue types, and of individual cells across three different timescales from fast (seconds) to extremely slow (one year). Time correlation analysis revealed significant differences in time constant (up to 0.4 s) between the principal layers of the inner retina with the ganglion cell layer (GCL) being the most dynamic. At the cellular level, significant differences were found between individual GCL somas. The mean time constant of the GCL somas (0.69 ± 0.17 s) was ∼ 30 % smaller than that of nerve fiber bundles and inner plexiform layer synapses and processes. Across longer durations, temporal speckle contrast and time-lapse imaging revealed motion of macrophage-like cells (over minutes) and GCL neuron loss and remodeling (over one year).

**Conclusions:** Physiological activity of inner retinal cells is now measurable in the living human eye.

## 1.

## Introduction

The inner retina is composed primarily of ganglion cells (GCs) whose axons, somas, and dendrites locate to three distinct retinal layers: nerve fiber layer (NFL), ganglion cell layer (GCL), and inner plexiform layer (IPL), respectively. The central role of GCs in processing retinal images captured by photoreceptors^{1} has been extensively studied since the first observations of GCs by Cajal et al.^{2} However, much remains unknown about the GC neural circuitry and its vulnerability to aging and disease, in part because of our inability to observe the activity of these highly translucent cells in the living human eye.^{3}^{–}^{8}

Recent progress in high-resolution, high-contrast imaging has overcome the translucency barrier, enabling visualization of individual retinal neurons—most notably GCs—in living human retina.^{9}^{–}^{11} In particular, adaptive optics optical coherence tomography (AO-OCT) allows three-dimensional (3-D) imaging of the individual cells and structures that comprise the inner retina.^{10}^{,}^{12}^{,}^{13} While successful, such imaging has not revealed the physiological activity of these cells. Here we investigate a method that does, by extending AO-OCT to detect temporal cellular changes.^{14}^{,}^{15} We observe temporal dynamics in the same patch of retinal tissue at dramatically different timescales, from a fraction of a second (e.g., intracellular soma dynamics) to one year (e.g., soma loss and migration). We also observe temporal dynamics on the intermediate scale of minutes, most notably the motility of macrophage-like cells—bright, irregular star-shaped cells that sparsely cover the surface of the inner limiting membrane (ILM).

We use a new correlation analysis method to characterize the fast temporal dynamics of the inner retinal layers (NFL, GCL, and IPL) and individual retinal nerve fiber bundles (RNFBs) and GCL somas. Time-lapse imaging and a temporal speckle contrast method related (but not identical) to that used in OCT angiography allow us to characterize the intermediate dynamics of macrophage-like cells and GCL somas that occur over minutes. Finally, we use pairs of images acquired one year apart to demonstrate the slow dynamics that occur over a year (neuron loss and remodeling). The ability to measure *in vivo* a wide range of fundamentally different dynamics in the same tissue using the same AO-OCT system reflects the power of the method we have developed.

## 2.

## Methods

Description of methods is in three sections. Section 2.1 describes the Indiana AO-OCT imaging system used in this study. Section 2.2 lays out the imaging protocol, experimental procedures, and subject information. Section 2.3 describes the postprocessing methods that we developed for visualizing and quantifying the temporal dynamics of the targeted inner retinal structures and cells. Details of our postprocessing methods are presented in Appendices A–G.

## 2.1.

### Indiana AO-OCT System

The Indiana AO-OCT system used in this study is described in detail elsewhere.^{16}^{,}^{17} Importantly, the fiber-based system operated at a center wavelength of 790 nm and bandwidth of 42 nm (superluminescent diode, SM fiber output power of 20 mW, BLMD-S-HP3, Superlum, Ireland), with a theoretical axial resolution of $4.7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in tissue ($n=1.38$) and lateral resolution of $2.4\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ (beam diameter of 6.7 mm at the eye pupil). We used the system’s two-camera mode to achieve an image acquisition speed of 500K A-scans/s (Kocaoglu et al.^{17} described the available camera modes). We focused the system on the GCL to maximize the signal strength and image sharpness of this layer. Optical power delivered to the eye was below $430\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{W}$ and more than an order of magnitude below the maximum permitted by American National Standards Institute^{18} for all our imaging protocols, as described next.

## 2.2.

### Imaging Protocol

Two subjects were recruited for the study (see Table 1). Neither had a history of ocular disease. The older subject was being treated for ocular hypertension (above normal ocular pressure) but was otherwise normal. All protocols adhered to the tenets of Helsinki Declaration and were approved by the Institutional Review Board of Indiana University. We obtained written informed consent after explaining the nature of the study and possible risks. Prior to the imaging session, one drop of Tropicamide 0.5% was administered to the right eye for mydriasis and partial cycloplegia. Axial eye length was measured with an IOLMaster 500 (Zeiss, Oberkochen, Germany) and used to correct for axial length differences in scaling of the retinal images following the method of Bennett et al.^{19}

## Table 1

Subject information.

Subject | Agea | Axial eye length (mm) | Spherical equivalent power |
---|---|---|---|

S1 | 27 | 24.0 | 0 D |

S2 | 50 | 25.4 | −2.5 D |

We acquired AO-OCT volume videos 12 deg temporal to the fovea. We chose this eccentricity because GCL soma sizes at the macular edge are larger and more variable,^{10}^{,}^{20}^{,}^{21} making it easier to compare cell structure and temporal dynamics. Imaging protocols for the two experiments are summarized in Table 2. Imaging protocol A was used to track fast temporal dynamics up to 2.6 Hz (Nyquist frequency of the 5.3 Hz volume acquisition rate) over a 1-deg retinal field of view and with $1\text{-}\mu \mathrm{m}/\mathrm{A}$-scan lateral spacing. Imaging protocol B was used to characterize intermediate and slow temporal dynamics occurring over time durations of minutes (maximum of 16 min) and 1 year (352 days for subject S1 and 364 days for subject S2), trading off speed for a $2\times $ increase in imaging area ($1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}\times 1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$) while maintaining a sufficient volume rate to support effective 3-D image registration. All videos were time-stamped to confirm their acquisition times and to combine them for the three studies: fast, intermediate, and slow temporal dynamics.

## Table 2

AO-OCT acquisition parameters for retinal imaging.

Imaging protocol | Image field of view (deg) | No. of A-scans per volume | Volume acquisition rate (Hz) | No. of volumes per video | Video acquisition interval (s) | No. of videos per retinal location |
---|---|---|---|---|---|---|

A | $1\times 1$ | $300\times 300$ | 5.3 | 12 | $60\pm 46$ | 23 |

B | $1.5\times 1.5$ | $450\times 450$ | 2.4 | 11 or 12 | $45\pm 25$ | 15 to 20 |

Ca | $0.75\times 0.75$ | $150\times 150$ | 20 | 50 | $48\pm 51$ | 30 |

## 2.3.

### Postprocessing

AO-OCT volumes were registered in all three dimensions with subcellular accuracy—a process accelerated by our custom 3-D B-scan registration algorithm that registers individual fast B-scan images to a reference volume using 3-D cross correlation.^{22} Assessment of cellular structural information was enhanced by averaging registered AO-OCT volumes to reduce noise while preserving retinal content. Temporal dynamics of inner retinal layers were assessed on three different timescales: 0.38 to 2.25 s (fast), 0 to 16 min (intermediate), and 1-year interval (slow). Each timescale required different postprocessing as described below.

## 2.3.1.

#### Fast dynamics

To characterize the fast temporal dynamics (over seconds), we developed a temporal autocorrelation method (Appendix A) that quantifies temporal change in the spatial intensity pattern within an estimation window with center pixel coordinates $\overrightarrow{{r}_{c}}=(X,Y,Z)$ and dimensions $({N}_{x},{N}_{y},{N}_{z})$. Given a sequence of $T$ volume images we have a $T$-element time series in which each value is an $({\mathrm{N}}_{x}\times {\mathrm{N}}_{y}\times {\mathrm{N}}_{z})$-element vector. We compute a correlation coefficient function $\rho (\overrightarrow{{r}_{c}},\mathrm{\Delta}t)$ [Eq. (4) of Appendix A] between pairs of vectors at times $t$ and $t+\mathrm{\Delta}t$. This function is observed to decrease monotonically with $\mathrm{\Delta}t$ [see Figs. 1(f) and 1(l)]. The time constant $\tau $ [Eq. (5) of Appendix A] quantifies the rate of decrease. The method is considerably more complex than the conventional definition of the autocorrelation coefficient [Eq. (1) of Appendix A] in order to mitigate the effects of several key sources of error that are known to degrade the accuracy and robustness of correlation measurements: (1) biases generated by static retinal structure in the images [Eq. (2) of Appendix A], (2) bias and uncertainty generated by typical sources of measurement noise [$\beta $ term in Eq. (4) of Appendix A], (3) biases generated by information loss caused by eye motion, (4) error in fitting an exponential decay to the correlation [Eq. (5) of Appendix A], and (5) bias generated by residual eye motion after image registration [Eq. (6) of Appendix A]. Based on prior studies^{23}^{–}^{26} we assumed optical roughness of the imaged retinal layers to be greater than the imaging wavelength, i.e., that the speckle patterns in our images were fully developed. This assumption implies that time constant $\tau $ must be independent of tissue optical roughness, a point we revisit in Sec. 4.

Correlation coefficients were computed from volume videos acquired with imaging protocol A and characterized over the range of 0.38 to 2.25 s (determined by the protocol’s 5.3 Hz volume acquisition rate and 2.25 s video length). To evaluate the trade-off between signal-to-noise (SNR) ratio and spatial resolution, we compared time constants for different retinal layers and tissue types using estimation windows of three different sizes.

The first estimation window size (window #1) covered the entire lateral extent of the volume image (typically $300\times 300\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ or equivalently $300\times 300\text{\hspace{0.17em}\hspace{0.17em}}{\mu \mathrm{m}}^{2}$) and was 7 pixels ($6.6\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$) deep, thus including most of each retinal layer (NFL, GCL, or IPL) in depth without extending beyond it. This inclusion of hundreds of thousands of pixels yielded the most accurate time constant estimates (see Fig. 10 in Appendix A). On the other hand, each layer is composed of cellular structures of different tissue types that might have their own unique time constants, which this large window size could not differentiate.

To assess the dynamics of specific tissue types, we used a smaller estimation window size (window #2). It consisted of a $4\times 4\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ $(X,Y,Z)$ stack ($4\times 4\times 6.6\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{3}$). We used this smaller window to compute a separate value of the time constant $\tau $ for each $XY$ location at a fixed $Z$ (depth) corresponding to the center of each retinal layer. This permitted visualization of spatial variation in temporal activity across the lateral extent of each layer. This small window size results in noisier measurements and a reduced $\tau $ because noise decorrelates. To improve the SNR ratio of our method without sacrificing tissue specificity, we averaged $\tau $ across pixels of the same nominal tissue type, which were semiautomatically determined based on differences in tissue reflectance: NFL bundles, GCL somas, GCL vasculature, and IPL synapses and processes (without vasculature).

The final estimation window size (window #3) was based on the smallest size of GCL somas in our images ($7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$) and used to more precisely evaluate their temporal dynamics. We selected pixels within a $7\times 7\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ volume (corresponding to $7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in each lateral dimension and $6.6\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in depth) centered on the soma. The volumetric center of the soma was manually identified using a customized graphic user interface display window that presents real-time *en face* ($XY$) and cross-sectional ($XZ$ and $YZ$) slices of the AO-OCT volume image with cursor position superimposed.^{15} The $7\times 7\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ window size captured intracellular dynamics while avoiding contributions from adjacent structures in the GCL such as glial processes, vasculature, and extracellular space.

## 2.3.2.

#### Intermediate dynamics

To characterize the intermediate temporal dynamics (across minutes), we constructed time-lapse image sequences of the same retinal patch from AO-OCT images acquired at different time points using imaging protocol B. A time-lapse sequence was generated for each time interval (0 to 5, 5 to 11, and 11 to 16 min for subject S1; 0 to 4, 4 to 8, 8 to 13 min for subject S2) and the motion of each pixel was quantified using temporal speckle contrast, defined as the ratio of the standard deviation (SD) of the reflectance amplitude to its mean (see Appendix B). We tested this method on GCL somas and macrophage-like cells that were observed $5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ above the ILM.

## 2.3.3.

#### Slow dynamics

We characterized the slow temporal dynamics (across a year) by reimaging the same retinal patch 1 year later using imaging protocol B. Image volumes from the two time points were registered to each other using a two-step process: first, rigid displacements of the volumes were corrected using the MATLAB (The MathWorks, Natick, 2017) function *imregtform*, which iteratively optimizes image similarity using Mattes’ metric.^{27} Second, nonuniform pixel-level displacements (image warp) were corrected using the MATLAB function *imregdemons*, which iteratively optimizes local image similarities using diffeomorphic demons algorithm.^{28}^{,}^{29} We then analyzed the cellular-level changes between registered volumes over the intervening year. This involved visually inspecting the registered images in rapid alternation and detecting difference in cell locations.

## 3.

## Results

## 3.1.

### Fast Temporal Dynamics of Inner Retinal Layers, Isolated RNFBs, and Individual GCL Somas

Figure 1 shows the cellular structures and corresponding fast temporal dynamics (correlation coefficients) of the inner retinal layers (NFL, GCL, and IPL) of the two subjects as obtained from AO-OCT. Individual GCL somas and RNFBs are clearly delineated in the intensity images [Figs. 1(d) and 1(j)] and [Figs. 1(c) and 1(i)], respectively. Vasculature (capillaries, arterioles, and venules) are also evident in all three retinal layers. Temporal correlation coefficients $\rho (\overrightarrow{{r}_{c}},\mathrm{\Delta}t)$ of the entire layers computed using window #1 ($300\times 300\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ stack) are shown in Figs. 1(f) and 1(l). NFL and IPL exhibit similar temporal dynamics, whereas GCL is clearly faster. The full-layer correlation decay time constants $\tau $ [Fig. 3(a)] confirm that GCL dynamics are $\sim 33\%$ faster than those of NFL and IPL.

A two-way analysis of variance (ANOVA) tested for variations in $\tau $ with retinal layer and subject. We found a main effect of retinal layer to be significant, $F(\mathrm{2,470})=2040$, $p<0.001$, where 2 is the degrees of freedom of the three retinal layers and 470 is the degrees of freedom of the 476 measurements with eight total number of levels (two subjects, three layers, and three interactions). However, this main effect was qualified by a significant interaction between the retinal layer and subject, $F(\mathrm{2,470})=12$, $p<0.001$. There was no main effect of subjects, $F(\mathrm{1,470})=0.1$, $p=0.75$. Bonferroni-adjusted comparisons indicated that the time constant of GCL was significantly faster than that of NFL [$p<0.001$, 95% confidence interval (CI) of the $\text{difference}=-0.38$ to $-0.32$ for subject S1; $p<0.001$, 95% CI of the $\text{difference}=-0.43$ to $-0.37$ for subject S2] and IPL ($p<0.001$, 95% CI of the $\text{difference}=-0.39$ to $-0.34$ for subject S1; $p<0.001$, 95% CI of the $\text{difference}=-0.39$ to $-0.33$ for subject S2) of both subjects. Repeated measures on the same retinal patches gave the same results with retina layer (see Appendix G).

We assessed specific tissue-type dynamics by computing the autocorrelation using a small estimation window centered on each $XY$ pixel location of each retinal layer (window #2—a $4\times 4\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ stack as described in Sec. 2). The resulting spatially resolved time constants for the three layers and two subjects are shown in Figs. 2(a)–2(f) as grayscale maps with corresponding histograms in Fig. 2(g). Time constants averaged over pixels of each retinal layer are plotted in Fig. 3(b) and over specific tissue types within a layer in Fig. 3(c).

As expected, interlayer and intersubject differences show the same trends as obtained with the larger window size (window #1). However, because the smaller amount of signal pooling inherent in the smaller window #2 produced measurements with more (temporally decorrelated) noise, time constants were reduced overall by about 20%. This prevents comparison across different window sizes, but the smaller window produces spatially resolved time constant measurements and allows us to compare the dynamics of different tissue types in the same layer.

A two-way ANOVA tested for variations in $\tau $ with tissue type and subject. Both main effects were significant: tissue type {$F(\mathrm{3,1}\mathrm{e}5)=3.4\mathrm{e}3$, $p<0.001$} and subject {$F(\mathrm{1,1}\mathrm{e}5)=143$, $p<0.001$}. However, these main effects were qualified by a significant interaction between the two, $F(\mathrm{3,1}\mathrm{e}5)=19$, $p<0.001$. Bonferroni-adjusted comparisons indicated that time constant of GCL vasculature was significantly faster than that of RNFBs {$p<0.001$, 95% CI of the $\text{difference}=-0.53$ to $-0.46$ for subject S1; $p<0.001$, 95% CI of the $\text{difference}=-0.55$ to $-0.47$ for subject S2}, that of IPL synapses and processes (without vasculature) {$p<0.001$, 95% CI of the $\text{difference}=-0.41$ to $-0.35$ for subject S1; $p<0.001$, 95% CI of the $\text{difference}=-0.37$ to $-0.30$ for subject S2}, and that of GCL somas {$p<0.001$, 95% CI of the $\text{difference}=-0.10$ to $-0.03$ for subject S1; $p<0.001$, 95% CI of the $\text{difference}=-0.09$ to $-0.02$ for subject S2} of both subjects.

We assessed individual GCL soma dynamics by computing the autocorrelation using a small estimation window centered on each soma (window #3—a $7\times 7\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ stack as described in Sec. 2). Figure 4 shows the results, color-coded and superimposed on the corresponding reflectance amplitude image. As seen in the Figs. 4(c) and 4(f) histograms, the time constant distribution is nearly unimodal with a mean and SD of $0.63\pm 0.12\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ for subject S1 and $0.74\pm 0.12\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ for subject S2 but positively skewed ($\text{skewness}=1.2$ and 1.1) and leptokurtic ($\text{kurtosis}=7.0$ and 5.3).

## 3.2.

### Intermediate Temporal Dynamics of Macrophage-Like Cells at the ILM, Somas in the GCL, and Vasculature Perfusion

Figures 5 and 6 illustrate the sensitivity of our method to temporal dynamics of macrophage-like cells, GCL somas, and vasculature perfusion over time durations of minutes. Figures 5(a), 5(b), 6(a), and 6(b) show time-lapse image triplets inserted into separate color channels of a single RGB image. Subtle movements of cells appear as color changes at the level of individual pixels ($1\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}=1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$), whereas white/black pixels indicate absence of motion (see figure caption, for details). The sequence of time-lapse images of the macrophage-like cells, as shown in Figs. 5(c) and 6(c) and the associated videos (Videos 1 and 2), further highlights the dynamics of these cells. A more quantitative and sensitive assessment of motion in terms of speckle contrast is given in Figs. 5(d) and 6(d), revealing activity over the entire footprint of the macrophage-like soma and processes. Figures 5(e) and 6(e) depict the expected result that speckle contrast highlights those vessels that are perfused and demonstrates good mapping to the vascular structure in the intensity images in Figs. 5(b) and 6(b).

Delineation of the macrophage-like cells for subject S1 was more difficult due to uneven ILM topography, a more reflective ILM and NFL located closer to the macrophage-like cells, and reduced image quality. The stronger ILM reflection in the younger subject (S1) is consistent with that expected in younger eyes.^{30} Despite these difficulties, we observed similar motion of the macrophage-like cells in both subjects.

## 3.3.

### Slow Temporal Dynamics of Macrophage-Like Cells, RNFBs, and GCL Somas

Figures 7 and 8 illustrate the capability of our method to assess slow temporal dynamics of macrophage-like cells and GCL somas by depicting pairs of images of the same patch of retina acquired 1 year apart. The associated videos (Videos 3 and 4) visually highlight the differences by presenting the images in rapid alternation. In contrast to the subtle motility-related changes that occur in macrophage-like cells over minutes (Figs. 5 and 6), the changes over a year are vast [Figs. 7(a), 7(e), 8(a), and 8(e)]. Macrophage-like cells appear to have been replaced between images. These cells are too active to be assessed at this timescale. In contrast, the intricate web of RNFBs in Figs. 7(b), 7(h), 8(b), and 8(h) appears stable over the whole year. Figures 7(c), 7(d), 7(g), 7(h), 8(c), 8(d), 8(g), and 8(h) demonstrate that we can also reimage the same patch of GCL somas a year later with striking one-to-one correspondence of somas and vasculature. Remarkably, there is little change in the cellular details of the GCL over this interval, indicating that the cell network is highly stable.

## 4.

## Discussion

## 4.1.

### Fast Temporal Dynamics (Seconds)

## 4.1.1.

#### Characterizing fast temporal dynamics

In the first part of this study, we characterized the fast temporal dynamics of the inner retina on three different spatial scales: (1) entire retinal layers, (2) structures of specific tissue types within the retinal layers, and (3) individual GCL somas. We quantified the dynamics in terms of correlation coefficients and time constants using a new correlation method (Appendix A) in conjunction with three averaging windows (denoted windows #1, #2, and #3) that defined the three spatial scales.

The largest averaging window (window #1: $300\times 300\times 6.6\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{3}$) permitted us to separate the contributions of the individual layers (NFL, GCL, and IPL) and to achieve exceedingly small 95% CIs, as shown in Figs. 1(f) and 1(l). Across the two subjects, average time constants, $\tau $, were $1.12\pm 0.02\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$, $0.74\pm 0.01\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$, and $1.10\pm 0.02\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ ($\text{mean}\pm 95\%\mathrm{CI}$) for the NFL, GCL, and IPL, respectively [Fig. 3(a)]. As evident in these figures, the GCL was significantly more dynamic ($\sim 33\%$ faster) than the NFL and IPL ($\mathrm{p}<0.001$).

Dynamics more specific to tissue type were obtained with smaller window #2 ($4\times 4\times 6.6\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{3}$) with the results given in Figs. 2, 3(b), and 3(c). As expected, the fast temporal dynamics of blood flow in the vasculature and their corresponding shadows produced the smallest time constants (darkest portions of the $\tau $ images). The $\tau $ over only the segmented vasculature in GCL was $0.46\pm 0.30\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ ($\mathrm{mean}\pm \mathrm{SD}$), close to the shortest time duration our method can resolve (0.38 s due to the 5.3-Hz volume acquisition rate in imaging protocol A). This suggests that the actual $\tau $ for the vasculature may fall outside our measurement range. The RNFBs were the most stable, appearing white in the figure and yielding the largest time constants ($1.05\pm 0.22\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ over segmented RNFBs), followed by the IPL synapses and processes ($0.91\pm 0.23\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ with vasculature contributions excluded) and the GCL somas ($0.62\pm 0.20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ with vasculature contributions and intercellular space excluded).

Histograms of the three time-constant images obtained using window #2 are shown in Fig. 2(g) with modes for NFL, GCL, and IPL of $0.89\pm 0.19$, $0.55\pm 0.13$, and $0.89\pm 0.17\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ for subject S1 and $0.87\pm 0.20$, $0.59\pm 0.15$, and $0.91\pm 0.17\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ for subject S2, respectively. While some of this variance is attributable to noise, the vast majority of the pixels in the spread had an intensity $>5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ above the noise floor. Thus we interpret these as signal and attribute them to differences in temporal dynamics of the tissue. Our method is therefore sensitive enough to measure local variations in dynamics within a single retinal layer.

Finally, we assessed dynamics of individual GCL somas using window #3 ($7\times 7\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$) centered on each GCL soma. As seen in Fig. 4, the mean time constants (0.63 s for subject S1 and 0.74 s for subject S2) are consistent with those from the GCL soma measurements shown in Fig. 3(c) using window #2. The $\sim 3\times $ difference between the least and most active GCL somas ($\tau \sim 0.4$ to 1.0 s for subject S1 and 0.5 to 1.3 s for subject S2) is notably larger than the 95% CI of our measurement ($\sim 0.05$ for subject S1 and $\sim 0.07$ for subject S2), indicating that we can detect differences in activity between somas.

We tested for correlations between soma activity and other soma parameters measurable in our AO-OCT volumes, namely size and reflectance. Soma size is of particular interest as it is a distinguishing feature of GC subtype (e.g., midget GC somas are smaller than parasol GC somas)^{10}^{,}^{20}^{,}^{21} at retinal eccentricities outside the fovea—such as the 12-deg eccentricity in this study. Figure 9 shows the resulting correlations. The time constant shows a weak positive correlation with soma radius, but it is significant for only one of the two subjects (${R}^{2}=0.01$, $p=0.19$ for subject S1 and ${R}^{2}=0.03$, $p=0.006$ for subject S2). Thus, smaller somas do not particularly exhibit slower or faster dynamics compared to larger somas, at least over the temporal range that we tested (0.38 to 2.25 s).

The time constant exhibits a weak but significant correlation with soma reflectance (${R}^{2}=0.19$, $p<0.001$ for subject S1 and ${R}^{2}=0.08$, $p<0.001$ for subject S2). Thus, somas with greater activity (smaller $\tau $) are generally less reflective (amplitude/pixel measured), perhaps suggestive of differences in the concentrations and distributions of organelles that move about within these cells. Finally, the strongest correlation we observed was between soma radius and reflectance. This positive correlation was moderate to strong and statistically significant in both subjects (Pearson, ${R}^{2}=0.17$, $p<0.001$ for subject S1 and ${R}^{2}=0.48$, $p<0.001$ for subject S2), meaning the larger somas are generally more reflective than the smaller somas. This result is consistent with our previous finding.^{10}

## 4.1.2.

#### Comparison to other retinal studies using speckle decorrelation

We know of no other study that has reported temporal correlation measurements of the inner retinal layers in the living human eye. It is, therefore, difficult to compare our measurements to the literature due to differences in the state of the tissue (*in vivo* versus *ex vivo*), species, retinal tissue type, and imaging and processing methods. These differences also confound attribution of our correlation measurements to subcellular activity. Nevertheless, a few comparisons to the literature are made.

The study most similar to ours, by Thouvenin et al.,^{31} measured the intracellular dynamics using correlation with full-field OCT of *in vitro* macaque and mouse retina. For macaque, temporal dynamics of 1 to 2 s or more were prevalent in RNFBs and IPL, whereas quicker dynamics of less than 1 to 2 s were dominant in GCL. They also observed that the dynamics of somas were faster than those of their surroundings. Our *in vivo* measurements in human showed the same trends with faster dynamics in the GCL than in the NFL and IPL and higher activity in GCL somas than in their surroundings. However, our measured dynamics were consistently twice as fast as theirs, perhaps due to differences in the measuring systems, experimental protocols, or state of the tissue (*in vivo* compared to *ex vivo*).

Lee et al.^{32} measured the intracellular dynamics of retinal GCs in extracted mouse retina using dynamic light scattering (DLS) OCT. They quantified the dynamics in these cells in terms of a diffusion coefficient that they reported as 1 to $4\text{\hspace{0.17em}\hspace{0.17em}}{\mu \mathrm{m}}^{2}/\mathrm{s}$. To compare, we followed Berne and Pecora^{33} by estimating an equivalent diffusion coefficient using the measured time constants of our GCL somas in Fig. 4 (see Appendix F). Based on this, the diffusion coefficient of our GCL somas was $8.0\pm 1.3\text{\hspace{0.17em}\hspace{0.17em}}{\mu \mathrm{m}}^{2}/\mathrm{s}$ (${\sigma}_{i}^{2}/2$), comparable but higher than that reported by Lee et al. We also computed the diffusion coefficients of RNFBs and IPL synapses and processes using the time constants in Fig. 2 (resulting in 5.3 and $6.2\text{\hspace{0.17em}\hspace{0.17em}}{\mu \mathrm{m}}^{2}/\mathrm{s}$, respectively) but have no measurements in the literature to compare to.

Speculation surrounds the attribution of these tissue dynamics. Thouvenin et al.^{31} suggested that the dynamics they observed arise from the active transport of organelles in the cells and possible cell membrane fluctuations and surface remodeling. In particular, the cytoplasm may provide a strong dynamic, owing to the trafficking of many organelles there. Similar to Thouvenin et al., Lee et al.^{32} suggested that their dynamics attribute to the intracellular motion of relatively large organelles (0.1 to $10\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in size, which covers the primary cell organelles). Activity-induced osmotic swelling may be another, which can occur locally within subcellular components (soma, dendrites, and axons). In photoreceptor cell outer segments, for example, growing evidence points to swelling as the dominating physical response of these cells.^{34}^{–}^{36} Given the general similarity of our *in vivo* measurements to the *ex vivo* ones of Thouvenin et al. and Lee et al. and the similarity of imaging methods, we expect our measurements to be sensitive to the same intracellular dynamics.

In other studies, Huang et al.^{37}^{,}^{38} reported much slower ($\tau \sim 34\pm 16\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$) and ($\tau \sim 59\pm 16\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$) speckle pattern dynamics in *in vitro* rat RNFBs and attributed them to axonal activity in the microtubules. This notable difference—their estimates were $>15\times $ slower than Thouvenin et al.’s and $>30\times $ slower than ours—suggests that the underlying nerve fiber bundle mechanism probed by these studies is either: (1) dramatically different in the two species; (2) slowed when removed from the eye, the extent of which depends on the method of extraction; or (3) actually two different mechanisms, perhaps because Huang et al.’s temporal sampling was 25 to $50\times $ coarser than Thouvenin et al.’s and ours (5 and 10 s compared to 0.01 and 0.19 s). More recently, a similar difference in $\tau $ was reported *in vivo* for IPL in mouse ($\tau =39\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$) using OCT by Zhang et al.^{39} Their temporal sampling was also considerably coarser than ours (12 s compared to 0.19 s).

Finally, two similar AO-OCT systems (including the one in this study) were used with similar scan pattern and sampling to measure organelle motility dynamics *in vivo* in human retinal pigment epithelial (RPE) cells of the outer retina.^{40} Their reported average time constant of $<5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ for RPE cells is notably slower ($\sim 5$ times) than for any of the inner retinal structures we measured ($\sim \text{\hspace{0.17em}}1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$). Similarity of systems and protocol indicate that the speed difference cannot be attributed to our use of AO-OCT. On the other hand, we have developed a more rigorous postprocessing method to handle retinal motion and the much weaker reflections from inner retina. This difference might explain the difference in results.

## 4.1.3.

#### Influence of eye motion and optical roughness assumption

We mitigated key sources of error that typically affect correlation estimates. These errors included (1) biases generated by static retinal structures, (2) bias and uncertainty caused by typical sources of measurement noise, (3) biases from information loss (image gaps) due to eye motion, (4) errors in fitting an exponential decay to the correlation, and (5) bias generated by residual eye motion (after image registration). While we were careful to address the most serious errors in our imaging study, errors could have still accrued. Two of concern include the effects of residual eye motion and our assumption of fully developed image speckle (i.e., that optical roughness of the retina is greater than the imaging wavelength). We discuss these sources of error now.

We corrected for effects of eye motion at subcellular resolution using 3-D B-scan registration and then further reduced them by time averaging over different combinations of AO-OCT volumes separated by the same $\mathrm{\Delta}t$. As stated above, we also accounted for eye-movement-related information loss within the estimation window. While we cannot be certain that these were sufficient, the correlation plots in Figs. 1(f) and 1(l) provide two lines of evidence to suggest that they were. First, the 95% CIs of the correlation coefficients (shaded colored bands) are exceedingly small (95% CI of $\rho =0.02$, on average). The significant differences we measured between the three layers’ correlation coefficients cannot be attributed to eye motion because they were imaged simultaneously (hence eye motion must be identical for all three layers). Second, the magnitude and pattern of eye motion is known to vary between subjects. Assuming this is true for our two subjects, a dominant effect of eye motion would have led to a difference in trends of their correlation coefficients. This was indeed observed and corrected by removing the effects of residual eye motion that manifest primarily as subpixel errors [Eq. (6) and Appendix E].

To simplify the interpretation of our results, we assumed the scattering properties of all of the layers and types of retinal tissues we examined to be “optically rough.”^{23}^{–}^{26} Optical roughness refers to retinal scatter within the coherence volume of the AO-OCT beam (nominally $2.4\times 2.4\times 4.7\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{3}$) that is dominated by optical path length differences greater than the AO-OCT wavelength (790 nm), thus leading to fully developed speckle. This assumption is commonplace in characterizing tissue because of the abundance of submicron-sized organelles that densely populate cells and are known to scatter light.^{25}^{,}^{41}^{–}^{43} Optical roughness or lack thereof could significantly affect the correlation coefficient estimate. In the presence of small tissue displacements (e.g., residual eye motion), speckle for optically rough tissue (e.g., organelle-filled cells) decorrelates more rapidly than for optically smooth tissue (e.g., surface membranes), regardless of whether the displacements are corrected in postprocessing (e.g., Appendix E). Subpixel displacement of the retina could therefore affect the correlation coefficient and time constant estimates differently depending on the scattering properties of the tissue. While it is our understanding that the three retinal layers (NFL, GCL, and IPL) and tissue types [RNFBs, GCL somas, GCL vasculature, and IPL synapses and processes (without vasculature)] that we examined are approximately optically rough, we did not attempt to measure optical roughness.

## 4.2.

### Intermediate Temporal Dynamics (Minutes)

We examined the temporal dynamics of macrophage-like cells, GCL somas, and vasculature perfusion over a time duration of minutes using time-lapse imaging and temporal speckle contrast analysis. Direct visualization and color coding of the time-lapse videos enabled us to detect micron-scale motion of macrophage-like cell processes [Figs. 5(a), 5(c), 6(a), and 6(c)]. To the best of our knowledge, these are the first observations of macrophage-like cell dynamics in the living human retina. The cells’ few stout processes and their apparent random distribution in a narrow region just anterior to the ILM suggest that they might be hyalocytes, a subtype of macrophage-like cell that typically resides in the cortical vitreous, either adjacent to or abutting to the retinal surface.^{44}^{–}^{47} Hyalocytes share a common origin^{44}^{,}^{48} and similar dynamics with microglial cells.^{46}^{,}^{47}^{,}^{49} Both are scavenging cells that continuously probe their local microenvironment. The cellular motion we observed was consistent with that reported for fluorescence-labeled microglia in *ex-vivo*^{50} and *in-vivo*^{51} experiments using mice, and so in our earlier report^{10} we interpreted the cells to be microglial or possibly astrocytes. However, microglial cells primarily populate the GCL and the inner and outer plexiform layers, so interpretation as hyalocytes is more plausible. Whatever their type, our methods have the sensitivity and temporal resolution to detect and track these cells *in vivo* as they scavenge about the retinal surface.

Our time-lapse results of the GCL [Figs. 5(b) and 6(b)] reveal—as expected—a highly stable GCL soma mosaic and vasculature over the time duration of minutes. These images show clear demarcation of the vasculature network, including small capillaries, but fail to differentiate perfused from nonperfused vessels. However, the temporal speckle contrast metric reveals subtler dynamics in the time-lapse videos [Figs. 5(e) and 6(e)]. The largest values are detected in the vasculature due to blood flow (color-coded as white), but we also observe an elevation or bias that permeates the entire GCL (GCL somas and extracellular space), color-coded as orangish red in the figures. This bias is greater than the system sensitivity as measured in the vitreous at the corresponding macrophage-like cell layer [Figs. 5(d) and 6(d)], color-coded as dark violet. The same layer demonstrates motion of individual macrophage-like cells, revealing activity over the entire footprint of the cell’s soma and processes and is consistent with our direct visual inspection of the time-lapse videos.

Interestingly, this speckle contrast metric permitted us to identify one nonperfused capillary. This capillary—located at the bottom of the speckle contrast image in Fig. 6(e) (yellow arrows)—is identifiable by its strikingly dark appearance compared to the other vessels. Dysfunction of this vessel is not evident in the corresponding color-coded time-lapse image [Fig. 6(b)].

## 4.3.

### Slow Temporal Dynamics (1 Year Interval)

In pairs of images acquired 1 year apart and inspected visually in rapid alternation, macrophage-like cells caused the most obvious image changes, with the same cells likely not present in both images (Figs. 7 and 8). This degree of activity is consistent with the scavenger role of macrophages, which are known to migrate across and through the retina. As they have been implicated in the pathogenesis of numerous retinal diseases, their numbers are believed to fluctuate as a function of retinal health. We can now measure their numbers and track them longitudinally.

The cells that compose the GCL appear highly stable over the 1-year interval except for abrupt changes associated with GCL soma loss and remodeling. For the GCL patches shown in Figs. 7(c) and 7(g) of subject S1 and in Figs. 8(c) and 8(g) of subject S2, we identified 831 and 589 somas, respectively, that were present at both times. We also identified one soma in subject S2 that was present in the first image but not the second. This missing soma is more salient in magnified view in Figs. 8(d) and 8(h), which also reveal migration of neighboring somas into the void created by the vanished soma that thus represents a form of retinal remodeling at the cellular level. Loss of 1 out of 590 GCL neurons is consistent with the histological reports of aging-related loss (0.19 to 0.72%/year^{52}^{–}^{56}) and within the range of loss rates that we have reported in an ongoing AO-OCT study of five different retinal locations in each of four normal subjects.^{57} The ability of our method to detect the loss of a single GCL soma demonstrates potential for extremely early detection of onset of GC-affecting diseases such as glaucoma.

Interestingly, we observed an unidentified dark globoid in subject S1 [Figs. 7(d) and 7(h)] that formed over the 1-year interval and appeared to have displaced adjacent GCL somas. The $24\text{-}\mu \mathrm{m}$ diameter feature generates bright reflections at its top and bottom boundaries. Its size and reduced internal reflectance are consistent with a displaced soma from the inner nuclear layer, e.g., a Müller cell (which are known to sometimes displace to the GCL^{58}^{,}^{59}). However, unlike neighboring GCL somas, no internal structure is evident, so it may instead be a fluid-filled vacuole (microcyst). The nearest vessel that could supply fluid is $\sim 20\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ away. A second smaller globoid is evident at the far left of the image (second red arrow). In subsequent imaging sessions of this same retinal patch, the larger globoid disappeared over 12 months, whereas the smaller one increased in size. We have since identified similar globoids in the GCL of other subjects.

## 5.

## Conclusion

We have developed a noninvasive method based on AO-OCT and a suite of novel postprocessing methods that measures both structural and physiological activities in retinal tissue down to the level of individual cells in the living human eye. The method was successfully applied to quantify the temporal dynamics of entire inner retinal layers, of specific tissue types, and of individual cells across three different timescales. Detecting physiological dynamics in this way offers the exciting possibility of longitudinally tracking very early cellular changes associated with disease onsets that cannot currently be detected clinically. This new capability also advances the prospects for noninvasively mapping functional aspects of neural circuitry in the living human retina.

## Appendix A: Temporal Autocorrelation Method

Temporal autocorrelation is already used in Doppler OCT and OCT angiography. These methods are designed to detect rapid changes that occur on the scale of microseconds to milliseconds,^{60}^{,}^{61} a range fast enough for measuring blood flow (e.g., velocity ranging from $\sim 1$ to $35\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}/\mathrm{s}$^{62}^{–}^{64}). For our application of measuring subcellular organelle motility, the changes we sought to detect are entirely diffusive (random motion; no flow) and occur over a much longer time period (couple of seconds). These translate into a variance rate from 5 to $29\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{2}/\mathrm{s}$ for our AO-OCT method (see Appendix F). Such slow dynamics (well below that detected by Doppler OCT and OCT angiography) exposes our method to eye motion artifacts and therefore requires careful attention to minimize these artifacts.

Starting with these initial requirements, we developed a correlation method using the conventional definition of the autocorrelation coefficient:

## Eq. (1)

$$\rho (\overrightarrow{r},\mathrm{\Delta}t)=\frac{{\u27e8A(\overrightarrow{r},t)A(\overrightarrow{r},t+\mathrm{\Delta}t)\u27e9}_{T}}{{\u27e8{A}^{2}(\overrightarrow{r},t)\u27e9}_{T}},$$^{33}(static structure; nonexponential decay), averaging and Rician- or Rayleigh-noise corrected correlation estimation (measurement noise

^{65}

^{–}

^{69}), and masked correlation

^{70}(information loss). We customized each of these to our application (imaging in the living human retina) and then combined them to create a unique sequence that mitigates our key sources of error. New are our methods to handle eye motion, including information loss and residual subpixel-level motion. These were implemented as follows.

First, to avoid biases generated by static retinal structure, time-invariant contributions in the AO-OCT volumes were removed by subtracting the time average of each pixel:

## Eq. (2)

$${A}^{\prime}(\overrightarrow{r},t)=A(\overrightarrow{r},t)-{\u27e8A(\overrightarrow{r},t)\u27e9}_{T},$$Second, we reduced the uncertainty and bias caused by noise. For the former, we performed both temporal and spatial averaging. Temporal averaging was realized by averaging $\rho $ over all possible volume combinations of $\mathrm{\Delta}t$ within a given volume video and across all volume videos acquired of the same retinal patch. For our study, the temporal sample size ranged from 23 to 253, depending on $\mathrm{\Delta}t$. Further improvement was gained by computing the mean correlation between ${A}^{\prime}(\overrightarrow{r},t)$ and ${A}^{\prime}(\overrightarrow{r},t+\mathrm{\Delta}t)$ across a 3-D spatial estimation window $w$ centered on pixel location $\overrightarrow{r}$. This entailed first removing spatially invariant contributions by subtracting the local spatial average: ${A}^{\u2033}(\overrightarrow{r},t)={A}^{\prime}(\overrightarrow{r},t)-\u27e8{A}^{\prime}(\overrightarrow{r},t){\u27e9}_{\overrightarrow{r}\in w}$, where $\u27e8{\u27e9}_{\overrightarrow{r}\in w}$ denotes the spatial averaging over $w$. With these changes, Eq. (1) becomes

## Eq. (3)

$$\rho (\overrightarrow{r},\mathrm{\Delta}t)={\u27e8\frac{\u27e8{A}^{\prime}(\overrightarrow{r},t)\xb7{A}^{\prime}(\overrightarrow{r},t+\mathrm{\Delta}t){\u27e9}_{\overrightarrow{r}\in w}}{\sqrt{\u27e8{[{A}^{\prime}(\overrightarrow{r},t)]}^{2}{\u27e9}_{\overrightarrow{r}\in w}{\u27e8[{A}^{\prime}(\overrightarrow{r},t+\mathrm{\Delta}t)]}^{2}{\u27e9}_{\overrightarrow{r}\in w}}}\u27e9}_{T}.$$In our study, we used $w$ sizes of $4\times 4\times 7$, $7\times 7\times 7$, and $300\times 300\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$, denoted as window #2, #3, and #1 in the main text. Because the average speckle size in our AO-OCT images was approximately equal to the nominal optical resolution of the AO-OCT ($2.4\times 2.4\times 4.7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in retinal tissue) and image sampling was $\sim 1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\text{pixel}$ in all directions, the effective number of independent samples occupying each window size was roughly 4, 13, and 23,000, respectively.

Noise can also create a bias in the correlation estimate that artificially attenuates $\rho $ ^{71}^{,}^{72} and reduces the corresponding estimated time constant, $\tau $ (defined below). This can be particularly problematic for additive pixel-correlated noise such as the photon noise, read noise, and relative intensity noise that are all present in our AO-OCT volume images. To correct this noise bias, we modified Eq. (3) as

## Eq. (4)

$$\rho (\overrightarrow{r},\mathrm{\Delta}t)={\u27e8\frac{{\u27e8{A}^{\u2033}(\overrightarrow{r},t)\xb7{A}^{\u2033}(\overrightarrow{r},t+\mathrm{\Delta}t)\u27e9}_{\overrightarrow{r}\in w}}{\sqrt{{\u27e8{[{A}^{\u2033}(\overrightarrow{r},t)]}^{2}\u27e9}_{\overrightarrow{r}\in w}\u27e8{[{A}^{\u2033}(\overrightarrow{r},t+\mathrm{\Delta}t)]}^{2}{\u27e9}_{\overrightarrow{r}\in w}}}\xb7\beta [\mathrm{SNR}(\overrightarrow{r},t)]\u27e9}_{T},$$^{26}

^{,}

^{41}

^{–}

^{43}

^{,}

^{73}

^{,}

^{74}More specifically, we first generated two correlated complex vectors ${\overrightarrow{\mathit{s}}}_{1}$ and ${\overrightarrow{\mathit{s}}}_{2}$, each of length equal to the number of speckles within an averaging window of set size. We controlled the correlation between the two vectors using Cholesky decomposition and also controlled the signal strength by multiplying both vectors by an adjustable constant. We also generated two uncorrelated noise vectors ${\overrightarrow{\mathit{n}}}_{1}$ and ${\overrightarrow{\mathit{n}}}_{2}$. These were added in the complex domain to the two signal vectors, and their absolute values (i.e., $|{\overrightarrow{\mathit{s}}}_{1}+{\overrightarrow{\mathit{n}}}_{1}|$ and $|{\overrightarrow{\mathit{s}}}_{2}+{\overrightarrow{\mathit{n}}}_{2}|$) were used to compute correlation coefficients. To model temporal averaging, we generated 1000 pairs of vectors and then computed the mean and SD of their correlation coefficients. While Fig. 10 shows only the case for the desired correlation equal to 1 for simplicity, we confirmed that this bias correction is valid for any desired correlation. Similar approaches have been proposed, for example, for estimating and correcting systematic error in diffusion tensor magnetic resonance imaging,

^{65}polarization-sensitive OCT,

^{66}

^{,}

^{67}and OCT angiography.

^{68}

^{,}

^{69}Here our method was designed for correcting errors in computing correlation coefficients based on AO-OCT amplitude signals.

Figure 10(a) shows the simulated correlation coefficient $\rho (\overrightarrow{r},\mathrm{\Delta}t)$ both without (red trace) and with (green trace) $\beta $ correction for different levels of spatial averaging ($n=4$, 13, and 23,000). The true value of $\rho (\overrightarrow{r},\mathrm{\Delta}t)$ was set to 1. This figure reveals two key findings. First, the notable differences between the red and green traces substantiate the importance of correcting for Rayleigh distributed noise bias, especially at low SNR. Second, the reduction in the widths of the green and red shaded regions as $n$ increases indicate that more precise estimates result from increased spatial averaging. Further improvement is possible by also performing temporal averaging, which reduces standard error of the mean (SEM) of the correlation coefficients. We employed temporal averaging in processing our data but not in the Monte Carlo simulation. Figure 10(b) shows the estimate of $\beta $ we used in our study, obtained with the largest window size ($n=\mathrm{23,000}$).

Note that the use of fixed windows for spatial averaging exposes our correlation estimates to eye motion biases. An eye movement can cause the imaging beam to skip over a region of retina, yielding a gap in the volume image within a given estimation window. This corresponds to error source #3: biases generated by information loss caused by eye motion. To avoid this window bias, we automatically masked those unimaged pixels in the averaging window.

Next, we estimated the time constant, $\tau $, by summing the time-averaged correlation coefficients across $\mathrm{\Delta}t$:

## Eq. (5)

$$\tau (\overrightarrow{r})=\sum _{\mathrm{\Delta}t=0s}^{2.25s}\rho (\overrightarrow{r},\mathrm{\Delta}t).$$This general expression for $\tau $ avoids the assumption of an exponential decay (error source #4), which our data did not follow. Note that the discrete integral in Eq. (5) theoretically underestimates $\tau $ because it is bounded at 2.25 s, the maximum time duration of our AO-OCT videos. While $\tau $ is sensitive to the maximum time duration (quantified in Appendix D), we used 2.25 s as it captures the period of most rapid change in the NFL, GCL, and IPL and also avoids unwanted disturbances, such as eye blinks and tear film breakup.

Finally, our image registration algorithm corrects motion artifacts as small as a single image pixel, but this leaves subpixel-level artifacts that can bias the correlation coefficient and time constant, which is also known as “spatial decorrelation noise” in Doppler OCT and OCT angiography.^{74}^{–}^{76} Because speckle noise dominates our AO-OCT images, we described this decorrelation bias as^{74}^{–}^{76}

## Eq. (6)

$${\sigma}_{\rho}^{2}=\mathrm{exp}[-\sum _{i=x,y,z}{\left(\frac{{\epsilon}_{i}}{{w}_{i}}\right)}^{2}],$$We computed the SEM to assess confidence limits of our $\rho $ and $\tau $ estimates. SEM of $\tau $ was determined by summing the squares of the $\rho $ SEMs and taking the square root.

## Appendix B: Temporal Speckle Contrast

We used temporal speckle contrast, an established motion metric based on time-varying speckle,^{73}^{,}^{77} to quantify intermediate AO-OCT image temporal dynamics (across minutes). Temporal speckle contrast is defined as the ratio of the SD of the reflectance amplitude to its mean:

## Eq. (7)

$$C(\overrightarrow{\mathrm{r}})=\frac{\sqrt{{\u27e8{[A(\overrightarrow{\mathrm{r}},t)-\u27e8A(\overrightarrow{\mathrm{r}},t){\u27e9}_{T}]}^{2}\u27e9}_{T}}}{{\u27e8A(\overrightarrow{\mathrm{r}},t)\u27e9}_{T}},$$^{26}

^{,}

^{41}

^{–}

^{43}This maximum value for the reflectance amplitude is equivalent to 1 for the corresponding reflectance intensity, a difference attributable to their different probability density functions (Rayleigh versus exponential).

^{25}

^{,}

^{41}

^{,}

^{42}

^{,}

^{73}Of practical significance, use of Eq. (7) on our AO-OCT data underestimates the true speckle contrast, a consequence of dewarping the B-scans in postprocessing to correct for nonlinearity of the scan pattern.

## Appendix C: Effectiveness of Our Method to Remove Structural Correlation Bias from Time Constant Measurements

We evaluated the effectiveness of our method [see Appendix A, Eq. (2)] to remove structural correlation bias (time-invariant contributions) from our time constant measurements. To determine, we computed the temporal correlation coefficient of NFL using window #1 ($300\times 300\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ stack) on the same dataset shown in Fig. 1(f) (subject S1), whereas varying the averaging time period $T$ used in Appendix A, Eq. (2). Figure 11 shows the resulting temporal correlation coefficient and corresponding time constant for $T=0$, 1, 2, 5, 10, and 15 min. As evident in the plots, short averaging periods ($T<5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{min}$) overestimate the time constant, but this error decreases asymptotically and becomes negligible for averaging periods $T>10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{min}$. We obtained similar results for GCL and IPL. Thus, we conclude that the $T=15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{min}$ averaging period used in our study is sufficient to remove structural biases in the layers, tissues, and cells we examined.

## Appendix D: Confirm Correlation Coefficient Reaches Zero for Long Time Periods

The correlation coefficients determined in Figs. 1(f) and 1(l) did not reach zero for the 2.25-s interval shown. Since in theory the coefficients must for a long enough time period, we tested this limit by recomputing the correlation coefficients after incorporating additional AO-OCT volume images that extended the time interval to 1000 s, almost $500\times $ longer than the 2.25 s used in our study. To extend, we concatenated videos that were acquired consecutively of the same retinal patch and recorded 30 to 240 s apart (see Table 2). Figure 12 (top row) shows our extended correlation results for the inner retinal layers (NFL, GCL, and IPL). As evident for both subjects, the layers decorrelate to less than 0.1 after $\sim 30\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ and reach zero after about 60 to 330 s depending on the layer. Thus, we confirm our method can detect zero correlation.

Next we computed time constants from the Fig. 12 (top row) correlation coefficient traces using Eq. (5) in Appendix A and plotted them in Fig. 12 (bottom row) on a log-log scale. As shown, the time constant strongly depends on the integration time period, varying from 0.3 s with the shortest integration period ($T=0.3\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$) to 8.4 s with the longest ($T=330\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$). This monotonic increase attributes from the shallow, nonexponential decay trace of our correlation coefficients that exhibit appreciable energy at low temporal frequencies ($\text{periods}\text{\hspace{0.17em}}>\text{\hspace{0.17em}}5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$). As such increasingly longer integration periods capture increasingly lower temporal frequencies that increase the time constant. As this dependence is fundamental to the correlation theory used, a common approach to circumvent it is to select a fixed integration period for all measurement comparisons. For this study, we selected $T=2.25\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$ as this captures the period of most rapid change in the correlation coefficients of NFL, GCL, and IPL.

## Appendix E: Determine Decorrelation Bias of Residual (subpixel) Eye Motion for Correcting Time Constant Measurements

Our image registration algorithm corrects motion artifacts as small as a single image pixel. We therefore expect ${\epsilon}_{i}^{2}$ (residual displacement error) in Eq. (6) of Appendix A to be limited by the sample spacing, i.e., the $1\text{-}\mu \mathrm{m}/\text{pixel}$ spacing used in our study. This limitation has been demonstrated in previous studies under the assumption of fully developed speckle.^{74}^{–}^{76} With this assumption, we estimated the size of ${\epsilon}_{i}^{2}$ by comparing correlation differences measured for two different sample spacings: 1 and $1.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\text{pixel}$. This approach follows that used in speckle metrology to measure surface roughness by purposefully changing speckle size, wavelength, or illumination angle by a known amount.^{26} To evaluate, we measured the same retinal patches of the same subjects using imaging protocols A and C (Table 1). Protocol C gave the same SNR ratio and A-scan exposure duration as protocol A but with coarser spatial sampling ($1.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{A}$-scan instead of $1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{A}$-scan). Because image registration accuracy is limited to pixel size, the coarser sampling should increase the residual displacement error (measured in microns, not pixels).

Figure 13 (top row) shows the resulting temporal correlation coefficients of the inner retinal layers (NFL and GCL) computed using the two protocols with window #1 ($300\times 300\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ stack and $150\times 150\times 7\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ stack, respectively). Time constants are listed in each plot. The coarser sampled measurements have time constants that are 27% and 14% lower in subjects S1 and S2, respectively. Assuming this reduction is due entirely to residual displacement errors as described by Eq. (6) of Appendix A, we solved for ${\epsilon}_{i}^{2}$, resulting in a ${\epsilon}_{i}$ of $0.86\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ (subject S1) and $0.60\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ (subject S2) for $1\text{-}\mu \mathrm{m}/\text{pixel}$ sampling. For both subjects, ${\epsilon}_{i}$ is near the $1\text{-}\mu \mathrm{m}/\text{pixel}$ sampling, indicating our image registration algorithm corrected eye motion at the pixel resolution of our method, as we had expected. Next, we estimated ${\sigma}_{\rho}^{2}$ for both subjects using the ${\epsilon}_{i}$ values with Eq. (6) and then corrected this bias in our original correlation coefficients traces [see Fig. 13 (bottom row), with and without the ${\sigma}_{\rho}^{2}$ correction]. As evident in the plots, correction of residual displacement errors results in shallower decaying correlation coefficients and increased time constants that are insensitive to sample spacing (1 and $1.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\text{pixel}$).

## Appendix F: Estimate Random Motion of Intracellular Scatterers from Time Constant Measurements

Temporal correlation measurements are sensitive to fluctuations in light scattered from moving particles (e.g., intracellular organelles) that occupy the coherence volume of the AO-OCT beam ($2.4\times 2.4\times 4.7\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{3}$). Following Berne and Pecora,^{33} we modeled this diffusive motion as a simple random walk and from which the correlation decay can be described by an exponential: $\mathrm{exp}[-\mathrm{\Delta}t\xb7\sum _{i=x,y,z}{({\sigma}_{i}/{w}_{i})}^{2}]$, where $\mathrm{\Delta}t$ is the temporal sampling interval of our AO-OCT ($0.19\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}/\text{volume}$), ${w}_{i}$ is the speckle size of our AO-OCT system [2.4, 2.4, and $4.7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ for $i=x,y,z$ in retinal tissue, and ${\sigma}_{i}^{2}$ is the random motion variance of scatterers in the retinal tissue that occurs per second in each direction ($i=x,y,z$)]. We relate this correlation decay expression to the time constant by $\mathrm{exp}(-\mathrm{\Delta}t/\tau )\sim \mathrm{exp}[-\mathrm{\Delta}t\xb7\sum _{i=x,y,z}{({\sigma}_{i}/{w}_{i})}^{2}]$ and solve for ${\sigma}_{i}^{2}$ to obtain an estimate for the random motion in each direction as described by ${\sigma}_{i}^{2}\sim \sum _{i=x,y,z}{w}_{i}^{2}/3\tau $ ($\mu {\mathrm{m}}^{2}/\mathrm{s}$). From this expression and our AO-OCT’s measured range of $\tau $ (0.38 to 2.25 s), we determined the range of random scatterer motion that we could measure: 5 to $29\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{2}/\mathrm{s}$ [or 2.5 to $15\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{2}/\mathrm{s}$ for the equivalent diffusion coefficient (${\sigma}_{i}^{2}/2$)].

## Appendix G: Repeated Measures of the Correlation Coefficients

In Secs. 3.1 and 4.1.1, significant differences were found between retinal layers, sublayer tissues, and individual cells. For example, the two-way ANOVA test for variations in $\tau $ with retinal layer and subject showed that residual sum of squares (equivalent to repeatability errors) ${\sigma}_{e}^{2}$ was just 10% of total sum of squares (${\sigma}_{\tau}^{2}$). Thus, a small fraction (10%) of the total variance is attributed to repeatability error.

To further substantiate this finding, we repeated the study by measuring the same retinal patches of the same subjects at four different time points using imaging protocol A (see Table 2) Each time point consisted of five to six videos acquired within 6 min. Results were analyzed for variations in $\tau $ with retinal layer and subject at each time. We tracked fast temporal dynamics up to 2.6 Hz (Nyquist frequency of the 5.3-Hz volume acquisition rate) over a 1-deg retinal field of view and with $1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{A}$-scan lateral spacing. For simplicity, we computed the temporal correlation coefficients of the entire layers using window #1. Figure 14 shows the time constants measured at the four time points. Despite fewer numbers of independent samples than in the main study (Sec. 3.1), we found the main effect of retinal layers to be significant for all the time points ($p<0.001$). The main effect of subjects and the interaction between subjects and layers were significant in two of the four time points but not for the same points. Bonferroni-adjusted comparisons indicated that the time constant of GCL was faster than that of NFL and IPL of both subjects at all time points ($p<0.001$). Thus, we conclude that differences between retinal layers are significant and repeatable. Differences between subjects are too small for our repeatability test, as designed, to measure.

## Disclosures

Donald Miller, Kazuhiro Kurokawa, and Furu Zhang have a patent on adaptive optics-optical coherence tomography technology and stand to benefit financially from any commercialization of the technology. Otherwise, none of the authors are aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this article.

## Acknowledgments

We thank Zhuolin Liu for early assistance on the project and Timothy Turner for software development. This study was supported by HHS | NIH | National Eye Institute, Grant Nos. R01EY029808 and R01EY018339.

## References

*in-vivo*cellular resolution neuronal and vascular retinal imaging,” Neurophotonics, 6 (4), 041105 (2019). https://doi.org/10.1117/1.NPh.6.4.041105 Google Scholar

*In vivo*measurement of organelle motility in human retinal pigment epithelial cells,” Biomed. Opt. Express, 10 (8), 4142 –4158 (2019). https://doi.org/10.1364/BOE.10.004142 BOEICL 2156-7085 Google Scholar

*In situ*characterization of the human hyalocyte,” Arch. Ophthalmol., 112 (10), 1356 –1362 (1994). https://doi.org/10.1001/archopht.1994.01090220106031 AROPAW 0003-9950 Google Scholar

*Ex vivo*dynamic imaging of retinal microglia using time-lapse confocal microscopy,” Invest. Ophthalmol. Vis. Sci., 49 (9), 4169 –4176 (2008). https://doi.org/10.1167/iovs.08-2076 IOVSDA 0146-0404 Google Scholar

*in vivo*show dynamic process motility at rest,” Invest. Ophthalmol. Vis. Sci., 58 (8), 316 (2017). IOVSDA 0146-0404 Google Scholar

## Biography

**Kazuhiro Kurokawa** received his PhD in engineering from the University of Tsukuba, Japan in 2013. Currently, he is a research associate at the Indiana University School of Optometry. His research primarily focuses on improving the capability of adaptive optics optical coherence tomography (AO-OCT) for visualizing and quantifying structure and function of individual retinal cells in the living human eye.

**James A. Crowell** received his PhD in cognitive psychology from U.C. Berkeley, working with Dr. Martin Banks of the School of Optometry on heading perception, signal detection theory, and ideal observers. He worked for 14 years as a visualization programmer supporting psychological research at the University of Illinois’ Illinois Simulator Laboratory. Currently, he implements ever-faster adaptive optics systems in the Miller Lab at the IU School of Optometry and contributes to papers and presentations.

**Furu Zhang** received his PhD in vision science from Indiana University under the mentorship of Professor Donald Miller. During his PhD, he focused on expanding the capabilities to image physiological activities of photoreceptors in living human retina using AO-OCT. Currently, he is a postdoctoral fellow at the U.S. Food and Drug Administration, working on developing advanced optical systems for retinal imaging and more sensitive optical biomarkers for retinal diseases.

**Donald T. Miller** earned his BS in applied physics from Xavier University and his PhD in optics from University of Rochester. Since 1998, he has been a professor at Indiana University School of Optometry. He has pioneered the development of high-resolution optical systems for imaging the back of the eye based on the combined technologies of adaptive optics and optical coherence tomography. He is a fellow of SPIE and the Optical Society of America.