Forward multiple scattering dominates speckle decorrelation in whole-blood flowmetry using optical coherence tomography.

Quantitative blood flow measurements using optical coherence tomography (OCT) have a wide potential range of medical research and clinical applications. Flowmetry based on the temporal dynamics of the OCT signal may have the ability to measure three-dimensional flow profiles regardless of the flow direction. State-of-the-art models describing the OCT signal temporal statistics are based on dynamic light scattering (DLS), a model which is inherently limited to single scattering regimes. DLS methods continue to be applied to OCT despite the knowledge that red blood cells produce strong forward multiple scattering. Here, we postulate that forward multiple scattering is the primary mechanism causing the rate of speckle-decorrelation derived from data acquired in vivo to deviate from the rate of decorrelation determined in phantom experiments. We also postulate that multiple scattering contributions to decorrelation are only present when the sample exhibits velocity field inhomogeneities larger than the scale of a resolution volume and are thus absent in rigid bulk motion. To test these hypotheses, we performed a systematic study of the effects of forward multiple scattering on OCT signal decorrelation with phantom experiments under physiologically relevant flow conditions and relative bulk motion. Our experimental results confirm that the amount of forward multiple scattering affects the proportionality between lateral flow and decorrelation. We propose that multiply scattered light carries information from different locations in the sample and each location imprints scattering dynamics on the scattered light causing increased decorrelation rates. Our analysis confirms that the detection of forward scattered light inside the vessel lumen causes an increase in the rate of decorrelation which results in an overestimation of blood flow velocities at depths as shallow as 40 µm into whole blood for OCT systems with typical numerical apertures used in retinal imaging.


Introduction
Optical coherence tomography (OCT) measures the path length of back-scattered light to obtain images of tissue microstructure in vivo [1,2]. OCT has become a standard of care in clinical ophthalmology as it can non-invasively measure and visualize retinal pathology. OCT angiography generates volumetric images and is used to map arteries and veins in the posterior eye. Angiographic maps are, however, essentially binary, indicating regions with or without flow [3,4]. There is evidence that ocular diseases such as diabetic retinopathy, macular degeneration, and glaucoma manifest through changes in the retinal and choroidal vasculature before physical symptoms develop [5][6][7]. This suggests that quantitative flow measurements in the retina could benefit diagnosis, risk stratification, disease monitoring, and assessment of treatment effectiveness.
Challenges in producing reliable measurements have hindered the development and clinical use of OCT flowmetry. Doppler signal analysis is the current research standard for measuring flow rates. However, there are drawbacks with this method when imaging retinal vasculature in vivo. This phase-difference contrast is highly sensitive to the angle between the flow and the beam, is only sensitive to the axial component of the flow, and requires phase stable acquisition [8][9][10][11][12]. The majority of retinal vasculature lies roughly perpendicular to the beam axis which makes the Doppler frequency shifts in these vessels very small with respect to the the large decorrelation of the phase induced by the lateral flow component. Although systems have been built to overcome these challenges [13][14][15][16], they require complex setups with additional illumination or detection optics which are not compatible with the constraints of clinical ophthalmic OCT systems. It has been suggested that methods based on signal decorrelation instead of the Doppler effect may provide more robust methods for calculating blood flow velocities in vasculature [17,18].
Here we denote all OCT flowmetry techniques based on the dynamics of the OCT signal as speckle-decorrelation methods: complex-amplitude techniques, such as dynamic light scattering OCT (DLS-OCT) and Doppler variance consider the decorrelation of the full complex OCT signal, while amplitude-based techniques make use of the OCT signal intensity or modulus [intensity-based DLS-OCT and split spectrum amplitude decorrelation angiography (SSADA)]. These analysis methods have been proposed as alternatives to Doppler OCT because of the ability to use the same established OCT hardware and the potential to provide -with an understanding of the impact of axial velocity gradients on the decorrelation [19]-angle-insensitive quantitative flowmetry. Several amplitude-based decorrelation techniques have been employed for imaging retinal blood flow [18,[20][21][22], but quantification has not been achieved. Complex-amplitude based speckle-decorrelation flowmetry gives access to both the lateral and axial components of the flow through the general DLS framework [23,24]. In this regard, DLS-OCT is the state of the art for measuring flow rates with OCT, but its accuracy is limited to single scattering regimes. Despite attempts to validate these methods in vitro, previously reported DLS-OCT studies have been limited to single-scattering liquids with scatterers that are small with respect to the size of red blood cells or with rigid bulk motion phantoms, which do not allow confirmation of decorrelation behavior with blood flow in vivo [23,25]. Measurements performed in vitro with blood have only been collected as an overall calibration tool with loosely-related decorrelation techniques such as SSADA [17,[20][21][22]. However, strong forward multiple-scattering (FMS) is known to affect the structural OCT signal [26][27][28]. Well-recognized decorrelation tail artifacts produced by FMS have been identified when imaging whole blood [17,23,29,30], but the impact of such artifacts inside vessels has not been examined in detail. The linear relationship between rate of decorrelation and flow velocity with DLS-OCT has been well established in a single-scattering regime [18,19,23,24,[31][32][33], but when imaging whole blood this relationship has been observed to exhibit a higher proportionality constant than predicted by the DLS-OCT model [34].
In this work we show how FMS plays a critical role in the decorrelation of the OCT signal inside the lumen of blood vessels at depths as shallow as 40 µm after light enters whole blood, implying that many current state-of-the-art DLS-OCT measurements result in an overestimation of blood flow. We demonstrate significant differences in decorrelation between flow with a single-scattering liquid and a liquid exhibiting FMS in vitro for physiologically relevant retinal blood flow rates. We further isolate the mechanism causing this increase in decorrelation by performing bulk-motion phantom-based experiments. These measurements confirm that FMS affects the signal decorrelation only in the presence of velocity field inhomogeneities, i.e. flow velocity gradients, and plays no role in the case of rigid-body translation. Together, these results highlight the importance of developing a theoretical OCT flowmetry model to account for FMS when measuring blood flow rates; we believe that the results from this work will be critical for informing these future models.

Single-scattering model
As is true with any form of coherent imaging, speckle -granularity in an image due to the microscopically heterogenous nature of the material from which the coherent light is scatteredis present in OCT tomograms [35]. A specific speckle pattern is given by the microscopic distribution of scatterers and analysis of the temporal fluctuations of speckle can be used to determine the dynamics of the scatterers in the sample. In DLS theory, the motion of flowing particles is derived from the time autocorrelation function of the signal [36] because the rate of change of the speckle pattern is related to the speed of the scatterers and the size of the resolution volume. Combining DLS theory with OCT, it is possible to calculate the autocorrelation function of the signal to spatially resolve flow profiles in three dimensions [23,31].
Theoretical models that use the autocorrelation of the signal to determine flow velocity are in good agreement with experimental measurements for flow with small scatterers in a singlescattering regime [19,23]. Here, we define the coordinate system (x, y, z), where z represents the depth coordinate, and x, y represent the plane perpendicular to the beam axis. The normalized first order autocorrelation function of the complex signal for a depth z with time step τ is defined where F is the complex backscattering signal of the sample, · · · represents an ensemble average, and * the complex conjugate. In DLS-OCT, F is given by the Fourier transform of the fringe signal. DLS has the potential to accurately determine both diffusive and deterministic motion of the suspended particles, with the condition that all decorrelation contributions are accounted for in the theoretical model. Gradients in the local axial velocity of the scatterers have been shown to impact the rate of decorrelation for both intensity and complex-signal calculations. In single scattering regimes, accounting for axial velocity gradients in g (1) (the autocorrelation of the complex signal) and g (2) (the autocorrelation of the intensity) have resulted in the recovery of accurate flow rates [19]. Under these conditions, in a single-scattering regime assuming Gaussian illumination and detection beam profiles as well as Gaussian wavenumber spectrum, the first order autocorrelation as a function of the decorrelation time τ can be written where i is the imaginary unit, n is the refractive index, k c the central wavenumber, D the diffusion constant of the scatterers, (w x ,w y ) the confocal diffraction-limited lateral point-spread function (PSF) e −2 beam radius, and w z the e −2 radius of the transform-limited axial PSF in the tomogram intensity. The first exponential term in Eq. (2) represents the Doppler component [18], the second term is a quadratic phase induced by axial velocity gradients [19], the third term describes decorrelation due to diffusion [23,37,38], the fourth term describes the increased decorrelation due to axial velocity gradients [19], and the last three terms describe decorrelation due to translational motion [18,31] -note the confocal effect contribution discussed below was missing in Srinivasan et al. [18] and first introduced by Weiss et al. [31]. The motivation for including the gradient terms and their implications on accurate flowmmetry are described by Uribe-Patarroyo et al. [19]. From the confocal effect on the lateral resolution of the system -and when illumination and detection paths are the same-we definew x,y = w x,y / √ 2, where (w x , w y ) are the diffraction-limited beam waist sizes of the illumination path only. If the lateral resolution is determined from e.g. the lateral speckle size (radius) in the tomogram intensity, the beam radiiw x,y are measured directly. We denote the e −2 radius of the confocal, practically-realized resolution volume in the intensity signal -affected by defocus, aberrations and axial aberrations (such as imperfect dispersion compensation or polarization mode dispersion)-with the effective lateralŵ x,y and axialŵ z parameters. We highlight that from this definition, althoughw x,y,z are system specific and spatially invariant,ŵ x,y,z are expected to vary spatially.
Similar to Ref. [19], ì v is defined as the mean voxel velocity and has components (v x , v y , v z ), ì ∇v z is the standard gradient operator applied to v z , and |ì ∇v z | 2 is the magnitude of a modified gradient operator |ì ∇v z | 2 = ŵ x ∂v z ∂x 2 + ŵ y ∂v z ∂y 2 + ŵ z ∂v z ∂z 2 which multiplies each derivative by the effective resolution volume size. The dependence on local axial velocity field inhomogeneities within the PSF volume means that speckle-decorrelation based flowmetry is sensitive to the angle between the flow and the OCT beam, however the impact on the rate of decorrelation is understood and therefore has the potential to be accounted for. We also emphasize this is a local effect that is not affected by large-scale inhomogeneities of the velocity field, only by the derivative of the velocity field inside the OCT resolution volume.

Multiple-scattering regime considerations
We define FMS as multiple scattering consisting of one of three scenarios: a) a forward scattering event followed by a back-scattering event after which light is detected by the OCT system; b) a back-scattering event followed by a forward scattering event after which light is detected by the OCT system; c) a forward scattering event followed by a back-scattering event followed by a forward scattering event after which light is detected by the OCT system. Events a) and b) are an anti-symmetric pair and should occur with the same probability, while event c) is symmetric. Forward scattering is defined loosely as scattering in the mostly forward direction. It is expected that the probability of events a) and b) would decay with the distance between the location of the two scattering events because the mostly forward scattered light has some probability of propagating laterally (through side-scattering) and therefore exiting the -generally weak-confocal gate. Our definition of FMS does not include side-scattered light or forward scattered light that exhibit significantly longer optical path lengths that appear to originate at deeper depths in the OCT tomogram. We also denote FMS as a non-local effect, in the sense that the different scattering events do not have to occur inside the resolution volume. We hypothesize that due to the combined coherent-confocal gating of the OCT system, most of the light scattered at a large angle from the propagation axis of the beam is rejected and the back-scattering of forward scattered light, rather than side-scattered light, is more likely to be detected. Thus, the contribution of FMS on the OCT signal becomes more significant as the phase function of the scattering particles becomes less isotropic. Due to the large size of red blood cells (RBCs) as compared to the OCT wavelength, light is scattered preferentially in the forward direction [39].
We postulate that RBCs produce significant amounts of FMS even within the normal range of hematocrit levels and at the numerical apertures (NAs) generally attainable in clinical ophthalmic OCT systems. However, there is no current framework that includes FMS in DLS-OCT analysis. In the experimental section of this work we show how important it is that the effects of FMS in the signal decorrelation inside the vessel lumen be understood and corrected in order to implement OCT flowmetry accurately in whole blood.
An important outcome from this work will be an understanding of the basic mechanism by which FMS alters the temporal dynamics of the OCT signal. In order to elucidate this, we designed the bulk-motion experiments section of this work. Rigid-body lateral bulk motion of the sample with respect to the OCT beam (or vice-versa, motion of the OCT beam with respect to the sample) is equivalent to the case of imaging a static sample with a full-field coherent imaging configuration. It is well established in full-field coherent imaging that the speckle size only depends on the numerical aperture of the system [35,40]. This is a result derived from the frequency content in the mixing of random scattering contributions at the lens aperture and the frequency filtering imposed by the finite size of the pupil of the optical system regardless of the nature of the process generating the individual random scattering contributions -single or multiple scattering [35]. In the case of lateral motion of the OCT beam with respect to a solid sample by means of, e.g., galvanometer mirrors, a change in the lateral scanning speed is equivalent to a change in the lateral sampling of the coherent imaging system mentioned above. Therefore, the lateral speckle size in an OCT B-scan necessarily depends only on the numerical aperture and the lateral scanning speed, regardless of the scattering process, focal plane position [31,33] or other optical aberrations [19].
However, this dependence only holds as long as the equivalence between an imaging system and the raster-scan system is true. A sample exhibiting diffusive motion or non-zero derivatives of the velocity field violates this: there is no equivalent imaging system for such a sample. As a result, multiple scattering can only affect the decorrelation of the OCT signal in the presence of non-local velocity field inhomogeneities or diffusive motion. Here we focus on the regime in which translational motion dominates diffusive motion. The above leads us to suggest that velocity field inhomogeneities are necessary to observe an increase in decorrelation through the following mechanism: because multiply scattered light carries information from different locations in the sample, and each location imprints its local scatterer dynamics on the scattered light, the presence of velocity field inhomogeneities in the flow produces a faster decorrelation when disparate contributions are added coherently. This has strong similarities to the mechanism by which local axial velocity gradients increase the signal decorrelation in single scattering regimes: signals from scatterers with dissimilar axial velocities add coherently, which results in faster decorrelation. However, in the case of multiple scattering we expect that all velocity field inhomogeneities, not only axial velocity gradients, contribute to the faster decorrelation. We expect that the scattering angle will deviate from exactly the forward direction, therefore inducing a Doppler shift with components from, in general, both axial and lateral velocities. Although we make no attempt to develop a theoretical model of the OCT signal in the presence of multiple scattering in the present work, our reasoning above defines the expected functional relationship of an additional factor to Eq. (2) to account for multiple scattering: such a factor would depend on the Jacobian of the velocity field ì v, the phase function of the scatterers, and the NA of the system. Finally, we want to emphasize the difference between decorrelation effects from local axial velocity gradients and velocity field inhomogeneities. The former, originating from differences in the relative axial motion of the back-scattering particles contributing to a single resolution volume, are known to produce increased decorrelation in single scattering regimes due to Doppler effects. The latter is defined in the context of our hypothesis on the increased decorrelation in FMS regimes. Axial velocity gradients contribute only in non-perpendicular flow conditions, while FMS is hypothesized to contribute in presence of velocity gradients in general. Given that the origin of axial velocity gradient is in dissimilar Doppler signals inside a resolution volume, they are a local effect. This makes them relevant to velocity changes at pathlength scales related to the optical path difference of detected multiply scattered light. In this work, we make every effort to avoid local axial gradient effects that arise in single scattering and focus on the interaction between FMS and non-local velocity field inhomogeneities.
Based on the above, we conclude that static bulk motion experiments cannot be used to gain insight into the effects of FMS on the decorrelation of the OCT signal. For this reason, we divide the methods and results sections of this work in two: Sec. 3 deals with the experiments assessing the effect of FMS in flow conditions, and Sec. 4 demonstrates the relationship between multiple scattering and velocity field inhomogeneities.

Materials and methods
In the following we consider an OCT system configured to acquire N t M-mode A-lines (A-lines at the same sample location) at N x different lateral locations, also known as a hybrid M-mode-B-scan. The complex OCT signal (after Fourier transformation of the fringe signal) F is a function of z (depth), x is the in-plane lateral coordinate (perpendicular to the flow direction), y is the out-of-plane lateral coordinate, t (time), f (frame), and polarization channel p = 1, 2 which corresponds to the channels of the polarization-diverse receiver. The frame index refers to the acquisition of multiple hybrid M-mode-B-scans without out-of-plane scanning.
If the OCT system is not phase stable, we perform a weighted least-squares linear fit in depth x, t, f , p)} to correct for phase shifts [41]. The weights for the linear fit are given 10 log 10 (|h|) − n dB , where n dB = 10 log 10 (n) is the intensity of the system's noise floor in dB determined by acquiring a tomogram with the sample arm light blocked and averaging |F| 2 across all A-lines. Each autocorrelation function for a given frame f , and polarization channel from the polarization-diverse receiver p was calculated as where τ is an integer time-difference index and the time index t starts at 0. We have N t samples in the final realization of g (1) , however we scanned multiple smaller correlation windows of size w t similar to the Bartlett's method in time-series analysis [42]. We corrected each realization of the autocorrelation function for the contribution from noise using the expression for the noise-corrected complex correlation coefficient introduced in Ref. [43], obtaining where the signal-to-noise ratio (SNR) is defined as SNR = |F| 2 /n − 1. We then performed a weighted average of the different realizations and the two polarization channels to improve the estimation of g (1) . The weights were defined to be the mean intensity of the signal inside each correlation windowĪ(z, Thus, the autocorrelation function becomes Instead of fitting the autocorrelation function to the model of Eq. (2), since contributions from Brownian motion are negligible in our regime, we defined a threshold g c = exp (−1), and calculated the time difference τ c where the autocorrelation dropped below the threshold after linear interpolation of g (1) . Given the behavior of the absolute value of the decorrelation is Gaussian, decorrelation time from a given threshold is related to the decorrelation time for any threshold. Instead on relying on a custom fitting procedure (see e.g. Ref. [23]), we have observed that since the behavior is Gaussian, the calculation of a decorrelation time from a threshold captures all the important information in a simple, quantitative manner. Given we avoid flow rates in which diffusion is expected to play a major role, we believe the use of a threshold is more efficient and conveys all the necessary information. The decorrelation time in time-difference units was converted into ms by considering the time scale t s defined as the inverse of the A-line frequency of the OCT system. We inverted τ c to calculate the rate of decorrelation τ −1 as This gives a pixel-by-pixel value for the decorrelation rate for each frame. With this equation, under lateral bulk motion conditions, the relationship between velocity and decorrelation rate reduces to v x = 1/2w x τ −1 [32]. Experiments in this Section were acquired with a wavelength-swept frequency domain OCT system with a 103 nm 3-dB unidirectional sweep centered at 1290 nm and a 54 kHz repetition rate. The system uses a frequency shifter for the removal of depth degeneracy [44]. The digitizer was configured to acquire 1536 samples per A-line, resulting in A-lines with 1024 pixels (px) in depth after zero padding and complex conjugate mirror removal, for a total depth range of 6 mm and a pixel scaling of 5.9 µm/px in air. The time resolution of the system t c is 18.5 µs. The sample arm was equipped with a 4 mm e −2 beam diameter fiber collimator [Schäfter + Kirchhoff GmbH, Hamburg, Germany, 60FC-4-M20-08] and a 36 mm focal length 5x scan lens [Thorlabs, Inc., Newton, NJ, USA, LSM03] with a theoretical e −2 beam radius of w x = w y ≡ w xy = 12 µm in air, resulting in a confocal lateral resolution ofw xy = w xy / √ 2 = 8 µm. This scanning lens was chosen in order to match the typical NA of clinical ophthalmic systems, which are known to have beam waist diameters in the range 15-30 µm [45]. We configured the galvanometer scanner to acquire hybrid M-mode-B-scans at 128 different lateral locations, each M-mode consisting of N t = 128 A-lines. We set w t = 64, thus obtaining N t /w t = 2 autocorrelation functions per M-mode which we averaged. The total data to calculate each flow frame was therefore 16,384 A-lines.
For the phantom, a polyethylene tube [Scientific Commodities Inc., Lake Havasu City, AZ, USA, BB31695-PE] with inner diameter 280 µm was embedded in a 1% weight/volume agar gel with 0.1% intralipid in volume. The length of the tube was placed on a polymer base to mimic the scattering properties of the choroid located beneath the retinal vasculature. In order to measure a Doppler signal within the same acquisition, we employed one of two techniques. With polystyrene bead flow, we used the phantom depicted in Fig. 1, where hybrid M-mode-B-scans were recorded at both locations indicated with arrows within one data collection. The rate of decorrelation was calculated when the tubing was horizontal and the Doppler signal was calculated in the region where the tubing was at an angle of ∼10 degs to the horizontal. To avoid coagulation induced by the introduction of physical angles in the tubing, when working with whole blood we removed the telecentricity of the scan lens by increasing the distance between the lens and the galvanometer scanning mirrors such that they were no longer at the back focal plane. This resulted in a varying angle of incidence as a function of beam lateral position, reaching an angle of ∼10 deg on the periphery of the scanning region as shown in Fig. 1. The focal plane was located within the vessel. Similarly to the experiments that used beads, this enabled us to acquire B-scans at different out-of-plane positions allowing us to extract the Doppler signal as a baseline measurement from which to compare the rate of decorrelation measured for the same flow. The Doppler signal was calculated using the Kasai autocorrelation using the full N t samples as the correlation window [32] and only the polarization channel with highest intensity. The angle of the tube with respect to beam (vertical axis) was used to calculate the absolute flow velocity from the Doppler axial velocity v z . Because it was difficult to guarantee a completely perpendicular angle, in order to rule out any potential effects from local axial velocity gradients, we limited the quantitative analysis to the center of the lumen.
First we assessed differences in the rate of decorrelation for similar flows between the erythrocytes and spherocytes (spherical RBCs) in vitro to determine if the rotation and deformation dynamics of non-spherical erythrocytes could be a significant source of decorrelation in our flow velocity regime. For the measurement of erythrocytes, we used whole swine blood at 45 % hematocrit, which is within the normal range for humans. We collected the precipitated RBCs and re-diluted them with 0.7% phosphate-buffered saline to make a spherocyte solution at approximately the same hematocrit. The spherical shape was confirmed with phase contrast microscopy.
Challenges in maintaining laminar flow with blood in small capillaries in vitro due to coagulation, led us to identify a more controlled surrogate to study the effect of RBC-induced FMS on the decorrelation in flow conditions. The chosen surrogate was a solution of 3 µm polystyrene beads [Polysciences Inc., Warrington, PA, USA] at a concentration of 1.68 × 10 −3 particles/µm 3 , which is similar to the concentration of RBCs in whole blood (45% hematocrit corresponds to 5.2 × 10 −3 particles/µm 3 , assuming an effective diameter of 5.5 µm to account for ellipsoid shape of the RBCs). This surrogate solution was estimated to result in ∼14 dB higher back-scattered signal intensity due to their higher refractive index (1.57 in contrast to 1.38 for RBCs at 1300 nm [46,47]), but similar scattering cross sections (σ s = 7.3 µm 2 for RBCs, σ s = 17.3 µm 2 for the 3 µm beads, as calculated by Mie theory), and similar scattering anisotropy (g = 0.98 for RBCs, g = 0.94 for the 3 µm beads). This resulted in 75% as much forward-scattering to total-scattering ratio as with whole blood. It also provided strikingly similar multiple scattering behavior in the tomograms to that seen in retinal flowmetry, such as similarly sized decorrelation tails. Solutions of 6 µm polystyrene beads were also considered, since they are closer to the true RBC long-axis diameter of 7.5 µm. However, due to the much higher refractive index of polystyrene than RBCs at the system wavelength, the 3 µm beads exhibit a scattering cross section much closer in area to that of RBCs. (σ s = 93 µm 2 for the 6 µm beads). Beads with a 6 µm diameter would provide ∼24 dB higher back-scattered signal intensity than RBCs, much higher attenuation, and only 73% as much forward-scattering to total-scattering ratio as whole blood. Additionally, the stock concentration of 6 µm beads was too dilute to create fully developed speckle [1.05 × 10 8 particles/mL] making an accurate determination of the signal decorrelation impossible. We also imaged 0.5 µm beads [1.82 × 10 10 particles/mL] to use as a reference for the rate of decorrelation in a single-scattering (SS) regime (g = 0.47, σ s = 0.009 µm 2 ) which we denote as SS beads. The differential scattering cross sections scaled by the appropriate concentration for all samples used or considered in these experiments are provided in Fig. 2. When imaging polystyrene spheres, we controlled the flow through the tube via an inline flow sensor and a pressure regulator [Elvesys Microfluidics Innovation Center, Paris, FR]. We also studied the impact of the thickness of FMS media traversed (blood layer thickness) on the rate of decorrelation by imaging 3 µm beads at the same concentration in tubes with 125 µm, 200 µm and 280 µm diameter lumen. In all cases the flow velocities were referenced with the Doppler measurement as described above.

Experimental results
Before starting the comparison between single and multiple scattering decorrelation regimes, it is necessary to verify that there are no additional sources of decorrelation that can be explained in a single-scattering model present in RBCs in our flow regime. For this reason, we first assessed the decorrelation rate as a function of flow velocity in erythrocytes and spherocytes at a hematocrit relevant for whole blood concentrations. Spherocytes, being fairly spherical cells, should show reduced decorrelation contributions from rotational, tumbling and vibrational motions. If the decorrelation is similar between the two samples, it is safe to assume these decorrelation sources are minimal in our flow regime. Figure 3(a) shows a subset of structural OCT cross-sections (out of a total of 4 flow rates for each sample) of flow through a 280 µm lumen. The intensity of the tomogram provides the luminance, and the rate of decorrelation is overlaid as the hue. The hue is an isoluminant map described in Ref. [48]. The limits for both variables are indicated in the two-dimensional colorbar. The cancellation of the telecentricity of the scan lens had the undesired side-effect of reducing the tomogram signal, which explains the lower tomogram intensity compared to subsequent figures. The white square indicates the central region used to mitigate effects from local axial velocity gradients for the quantitative analysis. This region was 100 µm × 100 µm (32 px × 22 px in the x × z plane) in size. When comparing the mean decorrelation within this region with the mean Doppler signal for the same region, we see that the decorrelation rate for a given flow velocity was virtually the same for both samples. We estimated an uncertainty of ±11% in the Doppler velocity from error in the estimation of the Doppler angle, however this is a systematic error that would affect all values in the same way. As expected with the hybrid M-mode-B-scan acquisition, the signal from the static portions of the phantom does not decorrelate in time; this is seen as the absence of hue outside the lumen. Figure 3(b) shows the average autocorrelation function |g (1) | as a function of τ for the central region of interest and confirms that the autocorrelation function for erythrocytes and spherocytes exhibit a similar shape. While we aimed for the v z component to be zero in the setup, it was changing slightly during and between experiments which made direct comparison of the complex g (1) less intuitive. Plotting |g (1) | as a function of τ allows for the visualization of the different decorrelation rates without contributions from the varying Doppler signal. Note that the shape of the autocorrelation function is Gaussian with the rate of decorrelation increasing with increasing flow rate. Figure 3(c) shows the expected linear relationship between rate of decorrelation and velocity. From this we can infer that, at flow velocities relevant for mid-to large retinal vessels, the potential decorrelation due to tumbling, tank treading, and other shape transformations is not important. However, this should be reevaluated for flow velocities in capillaries. With confirmation that these effects are minimal, we can assume that the use of spherical particles is a good surrogate to study the decorrelation effects of FMS. Additionally, given the fact that whole blood decorrelation has been observed to be higher than that of intralipid in a single-scattering regime in qualitative experiments [17], it is possible to start recognizing that FMS and its relation to velocity field inhomogeneities are the primary culprits for the increased rate of decorrelation previously identified. Confirmation of this hypothesis was the goal of the next experiment. We now compare the decorrelation behavior of 0.5 µm beads at a concentration consistent with a single scattering regime with 3 µm beads at the concentration described in Sec. 3.1. Figure 4(a) shows exemplary cross sections of the same lumen with the single scattering 0.5 µm beads and the FMS 3 µm beads. Figure 4(b) shows the absolute value of the autocorrelation function for both the 3 µm and 0.5 µm beads respectively. In Fig. 4(c) we compare the mean rate of decorrelation to that measured in the whole blood. From this we confirm that while the relationship between the Doppler velocity and the decorrelation time is linear, the presence of FMS in the 3 µm beads and blood samples results in a faster decorrelation of the signal in time. This manifests as a larger slope in the relationship between the rate of decorrelation and the Doppler velocity with increased FMS. The slope for the 3 µm beads is the same as that of whole blood. Using Eq. (2), we calculated a prediction for the slope of the 0.5 µm SS beads to be 0.09 1/µm. This defines the expected relationship in a SS regime between speed in mm/s and the corresponding decorrelation rate in ms −1 for this system. We expect that the main sources of discrepancy between model and experiment in a SS regime to be imperfect knowledge of the exact NA of the system and statistical bias in the calculation of g (1) due to the relatively small correlation window [49]. In a final flow experiment, we wanted to evaluate the dependence of the increased decorrelation due to FMS with the blood layer thickness. This is important since there is variety in the size of retinal micro-vasculature. We imaged the same 3 µm bead phantoms in tubes with inner lumen diameters of 200 µm and 125 µm. The results are shown in Fig. 5 where, again, we show a subset of B-scans with the decorrelation overlaid as the hue. The region of interest at the center of the tube was scaled to be a constant fraction of the lumen diameter, in order to have the same immunity to axial velocity gradient contributions. Their sizes were 72 µm-by-72 µm (24 px × 17 px in x × z) for the 200 µm diameter tube, and 45 µm-by-45 µm (16 px × 11 px in x × z) for the 125 µm diameter tube. The mean decorrelation for the tubes with smaller lumen is plotted with the mean decorrelation for the 3 µm beads in the 280 um diameter tube already shown in Fig. 4 and shows the decorrelation as a function of Doppler velocity. As expected from the conclusions of Jansen et al. [34], the dependence on lumen size suggests that the more multiply scattering media crossed by the beam, the greater the impact of FMS on the rate of decorrelation. The inclusion of the single scattering reference from Fig. 4 gives us information about how deeply within a vessel the assumption of a SS regime holds. As seen from the 125 µm data, even at a depth of 62.5 µm (the center of the tube) ±22.5 µm (the depth range of the region of interest) the effects of FMS are important to consider; decorrelation is already 25% faster than in single-scattering media. We were not able to perform an experiment in a lumen small enough to obtain the same behavior with 3 µm beads as in single scattering. This is an important experiment to carry out in the future, using whole blood, in order to determine the maximum depth over which decorrelation in whole blood follows the DLS model. This finding also indicates that, in order to accurately calculate blood flow velocity from the rate of decorrelation, the size of the vessel must be measured using the structural OCT image.

Materials and methods
To understand the mechanism behind the increased rate of decorrelation observed within the lumen, we turned to bulk motion experiments to exclude effects from Brownian motion and to isolate the dependence of the increased decorrelation on non-local velocity field inhomogeneities. We first designed an experiment to confirm our assertion in Sec. 2 that Fourier optics predicts that the decorrelation rate in a single scattering regime and in a forward multiple-scattering regime are the same for bulk translational motion. Next, we moved to an experiment where two solid phantoms were moved rigidly at two distinct speeds with respect to the OCT beam. This experiment had the goal of elucidating the effect of FMS on the decorrelation rate in the presence of dissimilar velocities at a non-local level, as well as determining the role of the thickness of the FMS-producing upper layer in the increased decorrelation of the bottom layer.
The bulk-motion data were acquired with a wavelength-swept frequency domain OCT system with a MEMS-VCSEL source [Thorlabs Inc. Newton, NJ, USA], with a 100 nm 10-dB unidirectional sweep centered at 1310 nm and a 100 kHz repetition rate (this a time scale of t s = 10 µs). The digitizer was configured to acquire 3072 samples per A-line, which resulted in an A-line with 2048 pixels in depth after zero padding and complex conjugate mirror removal, for a total depth range of 11 mm and a pixel scaling of 5.4 µm/px in air. The sample arm was equipped with the same fiber collimator and 5x focusing lens used in the flow experiments. This source had the added advantages of a faster repetition rate and phase stability, although the phase correction described in Sec. 2 was applied as well to correct for galvanometer-induced phase noise.
The first experiment had the goal of confirming that the decorrelation rate for bulk translational motion in SS and FMS regimes are the same. To this end, we fabricated two solid agar phantoms: a 1 % agar gel with 1 % intralipid in volume phantom (IL phantom exhibiting single scattering: g = 0.4, σ s = 0.001 µm 2 ) and a 1 % agar gel with 5 % 3 µm-beads in volume (FMS phantom: g = 0.94, σ s = 17.3 µm 2 ). Instead of moving the samples laterally, we left them fixed and moved the OCT beam with the galvanometer scanners at a constant speed. Compared with the experiments described in Sec. 3, given the homogeneous speed across the sample, we did not need to acquire hybrid M-mode-Bscans but instead configured the galvanometers to scan with uniform speed to calculate the autocorrelation function using a single B-scan as input; the autocorrelation was calculated across adjacent A-lines in a B-scan. We acquired 2048 A-lines per B-scan, but in the analysis omitted the first 1024 A-lines to avoid the galvanometer flyback and settling time as well as the last 24 A-lines, thus N t = 1000. The rest of the analysis followed as explained in Sec. 3.1. This included calculating a single g (1) per z, p and f , and calculating a weighted average among polarization channels. In this case w t = 100, therefore we averaged the resulting 10 different realizations. A single best frame in terms of sample homogeneity was chosen for the analysis. The rate of decorrelation τ −1 was calculated as in Eq. (6) giving the rate of decorrelation as a function of depth. The rate of decorrelation for the top 124 px (500 µm) -avoiding specular reflections-was averaged to give a single value for the rate of decorrelation for the analyzed region of the sample.
The second experiment was designed to highlight the impact of FMS layer thickness on the rate of decorrelation. We fabricated three phantoms: two "top" phantoms (SS and FMS phantoms), and a fixed "bottom" FMS phantom. A central section was carved out of both of the top phantoms in order to have a variable thickness of scattering material depending on the lateral position of the phantom. The top phantom was fixed in place and imaged with the OCT beam scanning laterally at a constant speed, reproducing the first experiment. In contrast, the bottom phantom was placed underneath the top phantom on a rotation stage moving at a constant rate which led to constant lateral motion in the imaging region. The net effect of this setup was to have two solid phantoms undergoing rigid translation at two distinct velocities, i.e. a non-local velocity inhomogeneity. Both phantoms were visible in the same OCT tomogram so that it was possible to determine the behavior of the decorrelation of the bottom phantom as a function of the scattering regime and thickness of the top phantom. We note that the choice to rotate, instead of translate, the lower phantom had the goal of easing data collection and ensuring consistency of the motion. We configured the galvanometer scanning and the rotation of the bottom phantom to closely match the direction of motion. That is, the scan direction of the galvanometers matched the direction of the tangent of the rotating phantom at the site of imaging. We scanned a region 1 mm in x (in-plane) and 2 mm in y (out-of-plane). The in-plane scanning served to provide a robust estimation of the decorrelation time, and the out-of-plane scan enabled us to see the effect of different thicknesses of the top layer. Each B-scan consisted of 4096 A-lines, providing sufficient sampling to accurately determine the signal decorrelation in this configuration. We selected N t = 2000 in each analyzed B-scan with w t = 100, and averaged a total of 38 z positions. We acquired B-scans at 8 y locations in the 2 mm range indicated above, and repeatedly acquired 8 B-scans at each location for a total of 64 B-scans. The acquisition time for the 8 y locations was 2.6 seconds, which given the mass of the rotating stage guaranteed the stability of the rotational velocity during the acquisition.
As a final step, because it was difficult to guarantee precisely the same rotational speed and off-axis position of the imaging site between experiments with different types of phantoms, we normalized the decorrelation rate of the bottom phantom τ b , measured as a function of top phantom thickness d in each experimental run, by its decorrelation rate in absence of the top layer . This allowed us to compare the behavior of the decorrelation time as a function of top layer thickness d even if the absolute velocity of the bottom phantom was different. This was necessary due the difficulty in reproducing exact absolute velocities between the different experiments.

Experimental results
The results from the first solid-phantom experiment are shown in Fig. 6 where we show a small region of interest of the tomogram of each (SS and FMS) phantom, as well as the calculated decorrelation rate encoded in the hue. Visually, the lateral speckle size is similar between the two samples. We filter the decorrelation rate in depth with a gaussian of radius of 20 px and a standard deviation of 12 px. The horizontal bands are due to the expected depth variation in the estimation of the decorrelation function given its stochastic nature. The average measured decorrelation rate over the region of interest is τ −1 SS = 2.6 ± 0.13 ms −1 and τ −1 FMS = 2.48 ± 0.18 ms −1 for the SS and FMS phantoms, respectively. Using Eq. (2), we calculated a predicted value of 2.12 ms −1 . However, in the present work, the relationship between decorrelation rates in a single and multiple scattering regime is the important metric. We highlight the fact that the difference between the two measured decorrelation rates is within one A-line acquisition time (t s ) for the OCT system.
The results from Fig. 6 confirm our assertion that FMS cannot increase the decorrelation rate in rigid bulk motion. The results from the second solid-phantom experiment are shown in Fig. 7. It is immediately clear that the decorrelation rate of the bottom phantom depends on the type of scattering produced by the first phantom and, in the case of FMS, on the thickness of the FMS-producing material. Figures 7(a)-(b) confirm that the decorrelation rate of the bottom phantom is the same independent of the presence of a singly-scattering top layer. Note that the decorrelation of the bottom layer is necessarily higher than that of the top phantom because it is moving at a higher speed due to the added velocity from the rotation. Also note that, as expected from the results of Fig. 6, the decorrelation in the upper layer is the same independent of the composition of that layer. More importantly, because we did not observe an increase of the decorrelation rate with depth in Fig. 6, we can conclude that the mere presence of a top Fig. 7. Bulk motion phantom experiment with two distinct velocities to determine the impact of FMS layer thickness on the rate of decorrelation. (a) Intensity B-scans with decorrelation rate overlay for intralipid (SS) upper layer, no upper layer, and 3 µm polystyrene beads. Note that the rate of decorrelation in the bottom layer for the single-scattering and no-scattering layers are similar, while the decorrelation in the bottom layer when the FMS layer is on top is higher. This demonstrates that forward multiple-scattering in the presence of the non-local velocity field inhomogeneity causes an increase in the signal decorrelation and thus the measured flow rate. We show the normalized decorrelation rate in order to directly compare the behavior of the decorrelation time as a function of top layer thickness d even though the absolute velocity of the bottom phantom was not exactly the same between experiments. Scale bars are 100 µm in x and z. (b) The plot shows the dependence on the rate of decorrelation in the lower layer due as a function of the thickness of the upper layer for the intralipid layer in a SS regime (blue) and the 3 µm polystyrene beads in a FMS regime (red). The 95% confidence interval for the fit is indicated by the dashed lines. (c) Schematic of layered phantom setup, not to scale. FMS layer does not increase the decorrelation rate of layers below it. We conclude that when the top FMS layer is present above the bottom phantom [ Fig. 7(c)], the speckle in the bottom phantom decorrelates faster due to FMS from the upper layer and the two distinct velocities of the phantoms. These two experiments demonstrate that having multiple velocities within the sample is a necessary requirement to see decorrelation contributions due to FMS. Figure 7(d) shows the dependence of the signal decorrelation on the thickness of the upper layer. As the amount of FMS material crossed by the OCT beam increases, the decorrelation rate in the bottom sample increases. Whereas increasing the thickness of the SS layer does not produce such an effect. Due to the large amount of samples (N t = 2000), the SNR correction of g (1) is expected to work fairly well, and therefore we do not expect that the drop in the signal of the bottom layer will affect the determined decorrelation rate. If this were the case, the SS phantom would exhibit the same trend. Note that due to experimental constraints, the distance between the top and bottom layers was about ∼2.1 mm in air. Even at this distance, the FMS managed to significantly affect the decorrelation rate of the bottom layer. These results confirm that FMS in the presence of non-local velocity field inhomogeneities causes an increase in the decorrelation of the OCT speckle signal. This also confirms the behavior identified in Fig. 5, i.e. the more FMS media crossed, the faster the signal in that region decorrelates. The effect of FMS on the signal decorrelation will thus impact the gradient term in DLS-OCT theory and not the bulk motion components.

Discussion
Strong FMS is known to affect the structural OCT signal. However, prior studies have not confirmed that the decorrelation behavior seen with blood flow agrees with the current DLS-OCT theory. We discussed the effects of velocity field inhomogeneities and FMS at the center of the vessel as defined in Sec. 2.2. Our observations of the decorrelation pattern for flow with a component along the OCT optical axis (the geometry used for acquiring the Doppler reference velocity) indicate the local axial velocity gradient contribution scales in a similar way with the amount of FMS present. However, it is clear that a model accounting for the increase of decorrelation at the center of the vessel needs to be developed before the local axial gradient effects at the edges of the lumen are considered. Our results indicate that in order to obtain true flow velocities from OCT flowmetry, structural information about the vessel diameter has to be included in the model because the amount of FMS material traversed by the beam increases the effects of FMS on the signal. It is also clear that if the determination of the flow profile is of interest (given that blunt, non-parabolic flow profiles have been observed in capillaries [50,51]), the decorrelation rate increase with depth within the lumen observed in our data must be quantitatively analyzed to determine the functional relationship. We note that due to this increased decorrelation, the structure at the bottom of the vessels is distorted, and the exact location of the vessel wall is difficult to discern. We also point out that the increased decorrelation from FMS in the presence of velocity field inhomogeneities implies that the speckle evolves in blood flow at a faster rate than it translates, hindering en-face speckle-tracking techniques for determining blood flow.
In the flow experiments, the decorrelation observed below the lumen (decorrelation tail) (seen in Fig. 5 for the 125 µm tube, not shown for other diameters) has much faster decorrelation than the decorrelation measured within the lumen. From our results, we expect that the coherent addition of the static signal from back-scattering with the dynamic forward-scattering signal inside the vessel should produce a faster decorrelation due to the larger difference in velocities compared to two scattering events happening inside the flow. It is unclear how feasible it is to determine blood flow from the choroidal vasculature, given it is known that the retinal pigment epithelium layer produces multiple scattering [52]. However, we expect that for sufficiently high NA, FMS effects in blood can be minimized; in Ref. [30] an objective lens with ∼6× higher NA than a typical NA in clinical ophthalmic systems was used, and they found that the decorrelation tail had generally lower decorrelation than the vessel lumen. We speculate that as the size of the depth of field decreases, the FMS light that is accepted back and detected decreases, thus reducing the effect of FMS on the rate of decorrelation. However, higher-resolution retinal imaging is challenging in the clinic as the effect of optical aberrations or defocus on the FMS contribution would need to be quantified, especially given the weakening of the confocal gate with defocus. Therefore, although increasing the NA can mitigate the decorrelation effects of FMS, it is not a solution applicable to all scenarios.
Given the relationship between the degree of polarization (DOP) or similar depolarization metrics with the presence of multiple scattering in the OCT signal [52], we expect that evaluating the behavior of the DOP with FMS could provide clues on how to determine the amount of FMS in a given location of the tomogram without significant a priori information. However, the DOP is a cumulative metric, and therefore local depolarization metrics should be used to assess the amount of FMS at a local level. It is not known if the ratio between overestimation in whole blood and in beads would be maintained at all diameters. However, assuming it is, given the overestimation factors in flow velocity with the beads for the [125, 200, 280] µm diameter tubes were [1.25, 2.1, 2.3], the decorrelation overestimation factor for whole blood would be the same. A natural extension of the experiment shown in Fig. 7 is to vary the distance between the two phantoms to understand how the change in velocity per unit distance (magnitude of the velocity gradient) affects the FMS contribution to the decorrelation. Making the top phantom with RBCs could also be possible in order to have a more accurate characterization of the FMS effects in whole blood.

Conclusion
With this study of the effects of multiple scattering on the decorrelation of the OCT speckle signal, we saw that FMS by RBCs prevents the direct transformation from decorrelation rate to flow velocity in whole blood. We can conclude that FMS is a primary cause for increased decorrelation rate in blood flow, and depending on the vessel diameter, can dominate the decorrelation. We showed the effects of FMS on the rate of decorrelation are dependent on the scattering properties derived from concentration and size of particles in the flow. We confirmed that FMS affects the decorrelation of the OCT speckle signal in the presence of velocity field inhomogeneities due to velocity changes of the scatterers at length scales related to the optical path difference of detected multiply scattered light.
Accurate quantitative measurements of blood flow in vivo with speckle-decorrelation based OCT methods at depths as small as ∼40-50 µm into FMS medium are not possible at the typical NA of retinal systems without accounting for FMS. To develop truly quantitative speckledecorrelation based methods, the DLS-OCT model must be extended to include more than one scattering event in strongly forward scattering media such as blood. From this work, we expect the functional relationship in Eq. (2) to have additional factors that depend on the Jacobian of the velocity field ì v, the phase function of the scatterers, and the NA of the system. With this correction and the considerations in Sec. 5, quantitative blood flow measurements using OCT flowmetry could be possible. While our motivation here centered on retinal imaging for diagnosis and tracking of disease progression with speckle-decorrelation based flowmetry, all micro-angiographic applications for which OCT is currently employed could benefit from quantification of blood flow velocities. Current applications include dermatology [53], neurology [54], and oncology [55]. We believe that the conclusions derived from this work will be useful in the development of a new DLS-OCT model that considers FMS decorrelation contributions. This model will enable accurate flowmetry with OCT from the signal decorrelation and other structural parameters, such as vessel diameters, determined from the OCT signal itself.