Phase-Sensitive Measurements of Depth-Dependent Signal Transduction in the Inner Plexiform Layer

Non-invasive spatially resolved functional imaging in the human retina has recently attracted considerable attention. Particularly functional imaging of bipolar and ganglion cells could aid in studying neuronal activity in humans, including an investigation of processes of the central nervous system. Recently, we imaged the activity of the inner neuronal layers by measuring nanometer-size changes of the cells within the inner plexiform layer (IPL) using phase-sensitive optical coherence tomography (OCT). In the IPL, there are connections between the neuronal cells that are dedicated to the processing of different aspects of the visual information, such as edges in the image or temporal changes. Still, so far, it was not possible to assign functional changes to single cells or cell classes in living humans, which is essential for studying the vision process. One characteristic of signal processing in the IPL is that different aspects of the visual impression are only processed in specific sub-layers (strata). Here, we present an investigation of these functional signals for three different sub-layers in the IPL with the aim to separate different properties of the visual signal processing. Whereas the inner depth-layer, closest to the ganglion cells, exhibits an increase in the optical path length, the outer depth-layer, closest to the bipolar cell layer, exhibits a decrease in the optical path length. Additionally, we found that the central depth is sensitive to temporal changes, showing a maximum response at a stimulation frequency of around 12.5 Hz. The results demonstrate that the signals from different cell types can be distinguished by phase-sensitive OCT.


INTRODUCTION
Non-invasive, cellular imaging of the neuronal activity in humans is compelling for clinical diagnosis and basic science (1). Currently, only optical methods may achieve the necessary resolution for this purpose. In addition, the retina is the only site, where the central nervous system is directly accessible for optical imaging. Besides being part of the central nervous system (CNS) and thereby providing diagnostic information on neuronal diseases, the initial information processing of visual impressions, which is mainly performed by the bipolar, amacrine, and ganglion cells, can be studied here (2)(3)(4). Each of these cell classes can be further divided into different functional cell types, which are sensitive to specific features of the visual signal. Currently, around 13 different types of bipolar and around 30 types of ganglion cells are known. They are responsible for processing different aspects for the visual information, e.g., brightness, color, as well as spatial and temporal variations, thus splitting the visual information into different feature channels (2,5).
In recent years, there has been extensive research to identify these different cell types and assign them to different functional tasks (2)(3)(4)(5). Although many types have already been classified, there are still gaps in the knowledge about these different cells.
To access these cell classes one can use retinal imaging with functional contrast. Suitable functional signals in the inner layers of the retina could be detected with different imaging methods [reviewed in (6)]. Techniques capable of identifying different cell types make, for example, use of calcium imaging (7,8). However, these methods require insertion of genetically encoded calcium indicators and are thus only available for animal studies.
Several other, non-invasive techniques could be and partially were applied to the living human eye. For example, electroretinogram (ERG) measurements provide information about the overall activity of the retina (9), but cannot spatially resolve contributions from different cell types in the retina. Furthermore, neurovascular coupling using laser Doppler flowmetry (10,11) or optical coherence tomography angiography (OCTA) (12) can measure increased blood flow when the retina is stimulated. This has already been used to measure the sensitivity of the retina to temporal changes of the visual signal in cats (10,11) and mice (12). The techniques may also achieve a relatively high resolution by measuring blood flow in the capillaries. Still, these techniques do not allow distinguishing of the functional contributions of individual functional cell types.
Furthermore, it has been shown that light-stimulated activation of the retinal function may lead to a change in reflectivity of the IPL that can be measured with optical coherence tomography (OCT). These changes in reflectivity due to functional activation could be demonstrated in different animals, for example, ex vivo in frog retina preparations (13), in paralyzed tree shrews (14), and in a few studies even in vivo in human (15). However, such measurements of changing reflectivity exhibit a high noise level, making it particularly challenging to extract the functional information from in-vivo measurements.
One technique, that emerged in the past years for non-invasive detection of the functional changes is phase-sensitive OCT. This technique can detect variations in the optical path length between the different retinal layers with a resolution of a few nanometers (16). Additionally, due to the tomographic character of OCT, such variations can be assigned axially and laterally to specific locations with a resolution better than 10 µm, which makes it in principle possible to assign functional changes to single cells. This cellular phase-based functional contrast has so far been demonstrated for the photoreceptor outer segments (OS) with several OCT techniques including Fourier-domain scanning (17)(18)(19)(20), line-field (21), and full-field (16,22) systems. Recently, by using a full-field swept-source OCT (FF-SS-OCT) system, we were able to show functional phase-changes in the inner plexiform layer (IPL) (23), where the synaptic connections of the bipolar, amacrine, and ganglion cells are located. However, low contrast, an order of magnitude smaller optical path length changes, and heart-beat induced pulsatile retinal motion make this endeavor extremely challenging in the inner neuronal layers. In the first demonstration of functional signals in the IPL, the resulting noise was minimized by extensive averaging of the phase over several voxels in the axial dimension. Due to this averaging, the contributions from all different functional bipolar/ganglion cells were integrated and no differentiation of the functional change in the IPL was possible. Unfortunately, the structural sensitivity of FF-SS-OCT has been insufficient to image individual ganglion or bipolar cells, although single ganglion cell imaging was demonstrated for other OCT systems by averaging hundreds of images (24).
Here, we use a recently demonstrated phase evaluation that takes advantage of the statistical properties of random and systematic phase changes to minimize phase noise (25) and use the stratification of the synaptic connections to differentiate signals from different cell types in the living human retina. With this phase evaluation, phase averaging in the axial dimension can be avoided to a large extent, so that the functional changes could be assigned to their origin in the retina with much higher accuracy allowing the evaluation of sub-layers. At the same time, the reduced noise levels enable a better investigation of the lateral dependency of the signal in the IPL. Although, single cell functional imaging of the inner neurons is still not possible due to the limited structural sensitivity, we could use the increased axial resolution and the improved signal-to-noise ratio (SNR) to distinguish the contributions of different functional cell types to the change in optical path length. To this end, we exploit the interconnection behavior of the different cell types in the IPL. Since certain functional bipolar and ganglion cell types form their synaptic connections only in a specific depth of the IPL, the corresponding properties of the visual impression are also processed at a specific depth of the IPL (2-4). Basically, 5 different sub-layers of approximately equal thickness can be distinguished histologically (26) and with visible light OCT (27). Thus, it can be expected that this stratification of the IPL enables a differentiation of functional cell types without resolving individual cells by investigating the phase changes in different depths or sub-layers of the IPL.

RESULTS
Here, the functional response of sub-layers to different visual stimuli was investigated in order to see if we can distinguish the response of different bipolar and ganglion cell types. During OCT imaging, the retina of a volunteer with no known diseases or ametropia was stimulated with a white LED combined with a transmission mask that projects a square-shaped pattern at   (Figure 1A). At this position, very little distortion of the actual stimulation pattern and the functional OCT signal was observed (23) due to the wiring of the photoreceptors by the Henle fibers to the inner layers of the retina (28). In addition, at this position in the retina, the size of the dendritic fields of the ganglion cells are in the range between 100 and 200 µm for non-midget cells, or even lower for the midget cells (29). The size of the square shaped stimulation pattern is 1 mm × 1 mm, so that the resolution of the stimulated area in the IPL is expected to be sufficient to analyze the stimulation pattern. The irradiance on the retina was about 57 nW/mm 2 , which corresponds to a luminance of 140 cd/m 2 and is therefore in the photopic vision regime (30). In order to investigate the functional phase change in different sub-layers of the IPL, the lower edge of the IPL was segmented (see Section 5 for detailed information). Phase differences were calculated in depths with constant distance to the segmented lower edge of the IPL. Due to the limited axial resolution of the OCT system, it was only possible to evaluate the optical path length change of three different sub-layers in the IPL, in which the 5 histologically differentiable layers partially overlap (see Figures 1B,C).  Figure 2J. The response of the photoreceptors OS matches the essentially uniform stimulation pattern. The slightly washed out edges probably result from aberrations of the projection of the stimulation pattern on the retina and ocular motion. In contrast, the functional imaging of the IPL shows a pronounced signal at the edge of the stimulation pattern, while the optical path length changes in the center of the stimulation pattern are weaker for all three layers (Figures 2G-I). By this inhomogeneous behavior, the contrast at the edge regions of the stimulation pattern is increased and a prominent contour becomes visible. Furthermore, the functional changes in the sub-layers of the IPL also differ. The greatest change in the optical path length is in the inner sublayer (closest to the ganglion cell layer). Here, a change in the optical path length of approx. 25 nm is observed in the edges and an increase of approx. 15 nm in the center. In the central layer, only a change in the edges of the stimulation pattern is observed (approx. 10 nm), while no change is observed in the center. Additionally, in the outer layer, the optical path length decreases during stimulation. Here, it is not possible to distinguish between edges and center of the stimulation pattern. The time courses of the three layers are shown in Figure 2K (solid lines). The curves were obtained by averaging the phase differences over the central area of the stimulation pattern. Residual motion artifacts, caused for example by vascular pulsation, were corrected by subtracting a background that was obtained by averaging the phase change in an area in a safe distance (at least 200 µm) from the stimulation site.
In  is therefore about the same strength as the signal from the inner sub-layer.
Since the central layer is the only one that showed an increased response to flicker stimulation, it was investigated more closely to evaluate the sensitivity to different stimulation frequencies.
In order to do so, the stimulation was performed with a continuous illumination and 11 different flicker frequencies of square waveform ranging from 1 to 55 Hz with a duty cycle of 50 %. Consequently, the number of photons reaching the retina was kept constant for all frequencies and halved compared to the constant stimulus. In order to quantify the sensitivity of the phase changes in the central IPL to the temporal frequencies, only the phase changes after 8 s in the central part of the stimulation pattern were evaluated. Signal changes in the edge regions of the stimulation pattern were neglected for this investigation, since they were not effected by different stimulation frequencies in our observations. To reduce the influence of noise or remaining artifacts caused by the arterial pulsation, the optical path length changes over the last 625 ms were averaged. Figure 3 shows the response of the central layer of the IPL to the different stimulation frequencies. As expected, continuous or low-frequency stimulation with 1 Hz only resulted in a minimal change in the optical path length. The signal increased for higher frequencies (2 − 6.7 Hz) and reached a maximum of around 15 nm for stimulation frequencies of 12.5 − 16.6 Hz. At even higher frequencies, the response decreased again and reached an optical path length change comparable to a continuous stimulation at 50 Hz.

DISCUSSION
For a functional evaluation of the phase in the human IPL in vivo, it is particularly important to have sufficient stability of the phase. Besides the general requirements for a phase-stable OCT system, as it is needed for the functional evaluation of the photoreceptors, a robust phase evaluation is crucial for a stable measurement of functional changes in the IPL, since a fast decorrelation of the speckle pattern is observed there, which is likely caused by metabolic processes in the inner retinal layers.
In contrast to other techniques that attempt to more directly assess retinal function and have been mostly used ex vivo or in animals, FF-SS-OCT relies on the changes in the phase of the backscattered light. The technique relies heavily on suitable postprocessing of the data. As such, even on a state-of-the-art GPU (Radeon VII) accompanied by an Intel i7 8700 and 64 GB of RAM, the processing of the camera's data to suitably segmented OCT volumes takes almost 7 min, when excluding disc I/O. Afterwards almost 2 min are required for the phase evaluation (25), which is done on the CPU (Intel i7) and requires sufficient RAM to achieve this processing time.
The proposed technique is robust enough to be applied to the living human retina as it was demonstrated here. Previous techniques that classified and measured the different ganglion cell types are invasive and can only be applied in animal models or are even limited to cell cultures. An important disadvantages is still that it does not allow functional imaging of individual cells but rather has to rely on the stratification of the dendrite connections in the IPL.
The most important parts to obtain the results presented here is the phase-stable and fast data acquisition, combined with the phase evaluation taking statistical properties into account. In contrast, the exact algorithms used for registration and segmentation could be replaced by other algorithms, however, given the low SNR of the used OCT volumes, many existing solutions fail.
The observed functional changes can be well-explained by the stratification behavior of bipolar and ganglion cell types in the IPL. The OFF-bipolar cells, for instance, only connect in the outer layers (closest to the bipolar cells) with the OFF-ganglion cells, while the synaptic connection of the ON-bipolar and ganglion cells form in the inner layers (closest to the ganglion cells) (2, 3). Since the ON-cells get activated by light onset and the OFF-cells get inactivated, those two cell types react in opposite ways to the same visual stimulus. The increase of the optical path length observed here correlates with the stratification of the ON-bipolar and ganglion cells, while the observed decrease in optical path length coincides with the connection of OFFbipolar and ganglion cells. Furthermore, we observed a high sensitivity to flickering light in the central layer of the IPL, whereas the response of the innermost and outermost layers was not affected by changes in time of the stimulation. This behavior fits well to the known stratification of sustained and transient cells. Sustained bipolar cells, which process slow visual changes, stratify on the inner and outer edges of the IPL, while transient bipolar cells, which process fast visual changes, tend to stratify in the middle of the IPL (2). Since the functional response in the inner and outermost layer does not change under flicker stimulation, those signals are probably caused by cell types that are sensitive to a different property in the visual signal such as intensity or chromatic information.
The temporal dependence of the IPL response of the middle layer with a maximum between 10 and 20 Hz and a return to the background above 50 Hz fits well to whole retina measurements using flicker ERG, electrical ex-vivo measurements on individual cells (31,32), and measurements of increased blood flow (10). In all these experiments, a decrease in the functional contrast was observed for temporal frequencies at about 30 Hz and above. The frequency dependence measured here also follows the temporal contrast sensitivity curve of the human eye (33). It should be noted, however, that only slow changes in optical path length are measured here and that potential high-frequency changes in optical path length cannot be detected with the used volume rate (8 Hz). Furthermore, based on results obtained for photoreceptors by Tomczweski et al. (22), we expect the amplitude of the fast response to be more than one order of magnitude lower and, at these frequencies, it will be difficult to differentiate from heart-beat induced pulsation (34). Finally, while the observed response is linked to the neuronal activity, the time course of the observed path length changes may differ and show a slower temporal response.
A further interesting effect is the edge sensitivity of the functional changes across the stimulated area, which can be explained by the center-surround processing within the receptive fields of some bipolar and ganglion cell types (35,36). In this processing scheme, stimulation of the marginal (surround) region of the receptive field acts antagonistically to stimulation in the central region. Areas in which the center of the receptive field and only part of the surrounding are stimulated thus produce a greater functional change than areas in which both the entire center and surround are stimulated. Such a processing scheme makes the cells particularly sensitive to edges in the visual signal (37). However, either the center activation has a higher influence than the surround inhibition, or additional cells contribute that lack the center-surround processing within their receptive field, causing a functional signal in the center of the stimulation pattern.

CONCLUSION
Here, we demonstrated that FF-SS-OCT combined with suitable post-processing provides sufficient phase stability and spatial resolution to analyze the neuronal signal processing in the human retina. Although, the functional changes of the optical path length in the IPL detected in this study can so far not be directly correlated with a specific physiological origin or molecular process, they are in good agreement with the existing knowledge of the visual signal processing of bipolar and ganglion cells.
It seems promising that phase-sensitive FF-SS-OCT can be used in the future to study the function of different, so far unclassified cell types in vivo. This could help monitoring innovative clinical treatments for ophthalmic diseases such as stem cell or optogenetic therapy in the retina (38,39). Given the currently time-consuming nature of the measurements presented here, future studies will be required to simplify the procedure, confirm robustness and reproducibility of the signals in multiple individuals, and to evaluate the physiological parameters that can be extracted.
Differentiating functional groups of ganglion and bipolar cells does not need to resolve individual cells, but is possible by studying signals from different depths in the IPL, where their synaptic connections are stratified. This work thus demonstrates the potential of phase-sensitive OCT to study complex neurological function and signal processing in vivo, which opens avenues for functional studies in humans.

Setup
All measurement were done with a Michelson-interferometer based FF-SS-OCT [see e.g., (40)(41)(42)]. It uses a tunable light source (Superlum Broadsweeper, BS-840-1) with 51 nm tuning range centered at λ 0 = 841 nm, resulting in an axial resolution of 8.4 µm (FWHM) in air. The light source illuminates an area of 2.6 mm × 1.5 mm on the retina. The total radiant flux on the retina from the sample illumination was about 5 mW. The light backscattered by the retina was superimposed with a collimated reference beam and detected by an ultra-high speed camera (Fastcam SA-Z, Photron) at a frame rate of 60 kHz at 640 × 368 pixels. During the wavelength scan the camera acquired 512 images per volume, which results in a volume rate of 117 Hz at full duty cycle. The data had to be stored with 8bit precision on an internal storage of about 8 GB. With the image size used and the number of frames required for each volume, this ultimately restricts the number of volumes in one imaging sequence to about 70. To achieve 8 s of measurement time, a volume was triggered only every 125 ms, corresponding to a volume rate of 8 Hz. Longer time intervals were hardly useful, because macroscopic and microscopic motion frequently destroyed the phase correlation between successive volumes.
During measurements, the retina was stimulated with a white LED and a transmission mask, which projected a square pattern on the retina. The stimulation light was coupled via a cold mirror into the sample arm and an infrared long pass filter in front of the camera ensured that no stimulation light disturbed the OCT imaging. The stimulation irradiance was about 57 nW/mm. The conversion from the irradiance E e [ W m 2 ], which was obtained from the spectrum of the used LED, to the illuminance E v [ lm m 2 ] is done by the following equation: Here, V(λ) is the luminous efficiency function (43) and 683 lm W is a conversion factor. The conversion to luminance L v [ cd m 2 ] is done according to where A p is the pupil area for a dilated pupil (50 mm 2 ) and D = 17 mm is the distance from the pupil to the retina.

Measurements
All measurements shown here were carried out on a single female volunteer with no known diseases or ametropia. Written informed consent was obtained. Compliance with the maximum permissible exposure (MPE) of the retina and all relevant safety rules was confirmed by the ethics board of the University of Lübeck. All methods and measurements were performed in accordance with the relevant guidelines and regulations. The study was approved by the ethics board of the University of Lübeck (ethics approval Ethik-Kommission Lübeck 16-080). To optimize the light detection, the pupil of the volunteer was dilated by mydriatic eye drops to a diameter of approximately 8 mm.
Before the beginning of each measurement session, the subject was dark-adapted for at least 20 min and all measurements were performed in a darkened room. Between measurements, a break of at least 5 min enabled a complete regeneration of the stimulated area in the retina. All stimulation (continuous and flicker) started at the 5th recorded volume, and continued for the full duration of the measurement, so that the first 4 volumes served as baseline.

Data Processing
After the reconstruction of the OCT data, numerical corrections for axial motion, dispersion (44), and aberrations (41,45) were done for all 70 volumes in a dataset.

Registration
Afterwards, the movement between successive volumes was compensated by a non-rigid registration of the data. The registration utilizes correlation of sub-areas of the different volumes to determine their relative displacement and finally interpolates a three-dimensional displacement map between all volumes.

Segmentation
In order to investigate the functional changes in different layers of the IPL, the lower edge of the IPL was segmented in a volume after aligning the volumes according to their registration. The utilized segmentation was loosely based on the work of Kafieh et al. (46). First, data points with a particularly high axial gradient were selected from a coarse-grained volume. Afterwards, a 3Dgraph based algorithm called diffusion maps (47) was used. This algorithm maps the selected data points into a new coordinate system, where they are rearranged according their connectivity, i.e., based on their transition probability after multiple timesteps given only their transition probability in a single timestep. The single time-step transition probability of the graph was thereby determined from the local proximity of the data points to each other and the local curvature of the surface, assuming that transitions along a surface curvature are more likely than those that deviate from the local curvature. The mapped data points are then clustered by a density-based algorithm, called DBSCAN (48), separating them into different layers and finally interpolated to give the respective layers in the high resolution dataset. Based on the segmentation, the curvature of all 70 complexvalued volumes was removed using the Fourier shift theorem and the individual A-scans were aligned, resulting in the lower edge of the IPL being flat.

Phase Evaluation
After segmenting, the phase changes of the volumes were evaluated for different sub-layers. This phase evaluation was done by using an adaption of the extended Knox-Thompson algorithm (49)(50)(51)(52), which is described in more detail elsewhere (25). For the evaluation, the cross-spectrum U x, y, z, t, t is calculated, which contains all possible temporal phase differences between the times t and t + t for each position (x, y, z), U(x, y, z; t, t) = U(x, y, z, t)U * (x, y, z, t + t), where U x, y, z, t represents the complex-valued OCT volumes. Following, the spatial phase-differences were calculated between two layers Z 1 and Z 2 , each consisting of two voxels in depth, with constant distance to the segmented lower IPL edge. In order to receive a better SNR the phase-information for each layer were averaged over the 2 voxels in depth, according to In the upper layers of the retina, the phase evaluation is especially challenging, since metabolic processes and blood flow cause strong molecular movements here. For this reason, decorrelation of the speckles occurs very fast and makes the phase evaluation difficult.
To minimize this kind of noise, an ensemble average of the phase differences over several adjacent lateral points is formed by convolving the phase with a lateral Gaussian function according to U Z 1 −Z 2 (x, y; t, t) ≈ G σ 2 (x, y) * U Z 1 −Z 2 (x, y; t, t). (5) Here, * denotes a circular convolution and G σ 2 (x, y) is a Gaussian function with standard deviation of σ = 4 pixels. The resulting, averaged cross-spectrum U Z 1 −Z 2 (x, y; t, t) was then used to find a phase function φ x, y, t that is consistent with all entries of the cross-spectrum (25).
From the resulting en-face phase images, the area in which a functional change occurs were determined. The phases in the area were averaged using an adapted mask to obtain a more stable ensemble average. Additional movement in the IPL, for example by vascular pulsation and other metabolic processes resulted in a continuously changing and non-reproducible background, which is superimposed to the functional signal. For correction and removal of the background noise, the phase changes in a second area that was not stimulated was subtracted from the functional changes (see Figure 4). The mask for the functional change was obtained manually for each measurement since the stimulated area could change slightly. The mask for the background noise was the same for all measurements.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics board of the University of Lübeck. The participant provided written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
CP worked on the setup, collected, analyzed, and interpreted the data, worked on the post-processing of the data, and wrote the manuscript. HS worked on the setup, collected the data, and worked on the post-processing. KG worked on the post-processing. SH, LP, and DM worked on the setup and collected the data. YM interpreted the data. GH obtained the funding and interpreted the data. DH obtained funding, worked on the post-processing, and interpreted the data. All authors contributed to the article and approved the submitted version.