Three-dimensional, three-vector-component velocimetry of cilia-driven ﬂuid ﬂow using correlation-based approaches in optical coherence tomography

: Microscale quantiﬁcation of cilia-driven ﬂuid ﬂow is an emerg-ing area in medical physiology, including pulmonary and central nervous system physiology. Cilia-driven ﬂuid ﬂow is most completely described by a three-dimensional, three-component (3D3C) vector ﬁeld. Here, we generate 3D3C velocimetry measurements by synthesizing higher dimensional data from lower dimensional measurements obtained using two separate optical coherence tomography (OCT)-based approaches: digital particle image velocimetry (DPIV) and dynamic light scattering (DLS)-OCT. Building on previous work, we ﬁrst demonstrate directional DLS-OCT for 1D2C velocimetry measurements in the sub-1 mm/s regime (sub-2.5 inch/minute regime) of cilia-driven ﬂuid ﬂow in Xenopus epithelium, an important animal model of the ciliated respiratory tract. We then extend our analysis toward 3D3C measurements in Xenopus using both DLS-OCT and DPIV. We demonstrate


Introduction
Cilia driven-fluid flow is an important physiological process in numerous organ systems. Ciliary flow is responsible for clearance of mucus from the respiratory tract, movement of cerebrospinal fluid (CSF) in the ventricles of the brain, determination of left-right patterning in the embryonic node, and movement of ova in the Fallopian tubes [1]. Because ciliary flow results from the shearing action of many cilia along a complex geometrical surface, ciliary flow lacks certain symmetries, such as unidirectionality and axisymmetry [2][3][4]. While these symmetries can simplify the quantification and analysis of other types of flow, such as Poiseuille flow in arteries, they are typically not applicable in the context of cilia-driven fluid flow. As such, the most complete description of the flow generated by cilia is a three-dimensional, three-component (3D3C) vector field v(r). v(r) describes the velocity field v(r) = (v x (r), v y (r), v z (r)) at every point in space r = (x, y, z). That is, for a steady-state flow field described in three spatial dimensions, the fluid motion at each location in three-dimensional space is described by a three-component vector.
A number of optical techniques have shown utility in imaging microfluidic ciliary flow [1].
These techniques include brightfield, epifluorescence, confocal fluorescent microscopy, and optical coherence tomography (OCT). In particular, OCT has shown great promise for visualizing and quantifying ciliary physiology due to its combination of intrinsic contrast from scattering, relatively high resolution (∼1-10 µm), and inherent parallelization of axial scanning without mechanical movement of samples [4][5][6]. Additionally, within OCT, multiple velocimetry techniques have been implemented using OCT-based data, including Doppler [7], particle tracking velocimetry [4], digital particle image velocimetry (DPIV) [5,8], and dynamic light scattering (DLS)-based approaches [9]. Traditional Doppler velocimetry is perhaps the most common type of flow measurement in OCT but inherently only recovers a single component of velocity (axial or v z ). We previously demonstrated the implementation of particle tracking to quantify ciliary flow [4] as well as investigate ciliary physiology and pathophysiology [10]. Particle tracking, however, requires sparse seeding of particles, a condition that typically leads to relatively sparse measurements of the velocity flow field. More recent work has demonstrated measurements of total speed from the OCT signal in addition to the Doppler signal [9]. In particular, these newer techniques employ either the autocorrelation or cross-correlation of the OCT signal, and in some cases can be extended to generate directional velocity measurements. In this manuscript, we will focus on an autocorrelation-based technique we previously demonstrated, directional DLS-OCT [11], as well as a related cross-correlation based technique, DPIV.
DPIV refers to an image processing method in which regions of an image are tracked from frame to frame in order to estimate velocity ( Fig. 1(a,b)) [12]. In DPIV, images are divided into smaller regions of interest (ROIs), and ROIs of adjacent image frames are cross-correlated against one another. By calculating a correlation function g(x 0 , z 0 ; δ x, δ z, τ) between regions from frame to frame and estimating the spatial location that maximizes the cross correlation ( Fig. 1(b)), features can in effect be tracked in order to estimate a displacement, and thus a velocity.
If performed on two-dimensional imaging data, DPIV yields 2D2C velocity data. In principle, DPIV can also be performed on serial volumetric data to yield 3D3C velocity estimates [13,14]. However, if serial volume acquisition is not sufficiently fast, ensembles of particles may decorrelate due to non-translational motion including diffusion, shearing, and rotation. These non-translational motions then lead to a degradation of signal-to-noise ratio (SNR) [15]. Additionally, many imaging acquisition protocols and image processing algorithms are already available and optimized for two-dimensional data [16], and we wished to leverage these preexisting algorithms.
In order to circumvent the needs for repeated volume measurements, multiple measurements of lower dimensional quantities can be combined to synthesize higher dimensional data. In the simplest example, multiple one-component velocity measurements made at different angles can be combined into a single two-or three-component measurement [17]. Stereoscopic PIV, for example, uses two cameras angled with respect to one another to compute 2D3C data from 2D2C data [18]. Stereoscopic PIV as implemented with a brightfield microscope is not truly cross-sectional, however, as velocity measurements may arise from the movement of defocused particles. In contrast, OCT offers truly cross-sectional data, and thus the ability to localize velocimetry measurements in three-dimensional space. As such, OCT offers the possibility of generating 3D3C data from cross-sectional 2D2C data. We describe this process in further detail in the methods section.
An alternative approach to feature tracking methods such as DPIV is to quantify the decorrelation rate of a time varying signal to estimate speed [9,11,19,20]. In dynamic light scatteringbased OCT (DLS-OCT), the temporal decorrelation of a light-scattering signal is fit to a model autocorrelation function in order to recover the axial velocity and total speed of the underlying Fig. 1. In this study, we use two distinct correlation-based approaches to directional OCT velocimetry: DPIV (a,b) and DLS (c,d). (a) In DPIV, two frames of intensity information I(x 0 , z 0 ), acquired with an interframe time of δt, are compared against each. Each image is segmented into smaller regions of interest (ROI) I roi (x, y), centered around the point (x 0 , y 0 ), and cross-correlation is calculated in the spatial domain (δ x, δ z). The process is then repeated for multiple ROIs. (b) The corresponding correlation function g(x 0 , z 0 ; δ x, δ z, τ) for each region exhibits a peak where δ x = v x δt and δ z = v z δt. Estimation of this peak location yields the velocity estimate for (v x , v z ) (c) In DLS-OCT, the temporal signal at a single point (x 0 , z 0 ) is correlated with itself, except now the complex field signal E(x 0 , z 0 ,t), instead of the intensity, is used, and the correlation is calculated in the time domain. Additionally, scatterers are seeded significantly more densely such that speckle formation occurs. (d) The correlation function g(x 0 , y 0 ; τ) exhibits a peak at τ = 0, but the rate of decorrelation (γ), i.e. the inverse of the width of the Gaussian, is proportional to the difference in velocity of the scanner and fluid flow. By modulating the scan speed v scan x , we can minimize the decorrelation rate γ and estimate v f low x as the minimum of that curve. The process is then repeated at multiple spatial locations (x 0 , z 0 ). fluid ( Fig. 1 (c,d)). DLS is typically implemented in highly scattering media. In this highly scattering regime, each focal volume is occupied on average by multiple particles such that pixel intensities are random quantities governed by speckle statistics, and the decorrelation of these speckle patterns gives an estimate of total flow speed ( Fig. 1(d)).
We previously showed that by imparting a variable scan bias to our measurements, we could extend DLS-OCT to yield directional velocimetry measurements [11], a technique we term directional DLS-OCT. In our previous demonstration, the velocity component data was only two-component, and the measurements were only made along a single axial line (1D). However, by (a) scanning along an additional orthogonal plane to get 3C measurements, and (b) repeating these measurements at multiple spatial locations to get 3D measurements, 3D3C data can in principle be synthesized from lower dimensional measurements.
Here, we demonstrate two OCT-based, correlation-based methods of generating 3D3C vector flow fields without serial volume acquisition. First, we utilize existing DPIV implementations in order to measure 2D2C data, and we extend the measurements to synthesize 3D3C data using multiple 2D2C planes at two different angles. In the process, we investigate the signal-to-noise ratio (SNR) properties of this process to optimize those two angles. Secondly, building on our previous work in DLS-OCT, we show 1D2C DLS-OCT measurements in vivo in the Xenopus laevis model of epithelial ciliary flow. We then extend these measurements by implementing an orthogonal scan plane (yielding 3C measurements) at multiple spatial locations (yielding 3D measurements) in order to generate 3D3C data. We validate DPIV-OCT and DLS-OCT in a calibrated flow phantom and then demonstrate both techniques to quantify Xenopus epithelial flow. Finally, we show the possibility of measuring 3D3C flow in two additional ciliated systems using DPIV. We measure CSF flow in Xenopus tropicalis using endogenous particulate matter in the ventricles of the developing brain, and we measure ex vivo ciliary flow in mouse trachea using exogenous particle tracers.

OCT Data acquisition and processing
All OCT measurements were made using a spectral domain OCT system (Thorlabs Telesto) with 1325 nm center wavelength. Images were acquired at a 28 kHz linerate with an estimated sensitivity of 96 dB. The optical power incident on the sample was measured at 1 mW.

Animal measurements
All Xenopus and mouse procedures were reviewed and approved by Yale's Institutional Animal Care and Use Committee, which is Association for Assessment and Accreditation of Laboratory Animal Care-accredited. For epithelial flow measurements, Nieuwkoop-Faber (NF) [21] Stage 32 Xenopus laevis tadpoles were immobilized with benzocaine and imaged in a 6 mm diameter, 1 mm deep cylindrical well [22] as previously described [10]. Polystyrene microspheres (Bangs Labs, density 1.05 g/cm 3 ) of various sizes were used as a contrast agent. For DPIV measurements, 5 µm microspheres were diluted 1:100 in 1/9x Modified Ringer's (MR) solution. For DLS-OCT measurements, 500 nm microspheres were diluted 1:100 in 1/9x MR with 0.1% Tween added to prevent bead aggregation. Additionally, for DLS tadpole epithelial measurements, a liquid immersion adapter was used to minimize common path signal due to reflection off the fluid surface. For CSF measurements, stage 46 Xenopus tropicalis embryos were immobilized with benzocaine, but no exogenous scattering agents were added.
For tracheal flow measurements, the trachea from a one-month old mouse pup was prepared as in [23]. The pup was euthanized using intraperitoneal urethane. The trachea was removed by making two cross-sectional cuts: a cranial cut near the larynx, and a caudal cut above the level of the carina. The lumen of the excised trachea was then flushed with phosphate buffered saline (PBS) and filled with 5 µm microspheres diluted 1:100 in PBS. Temperature was maintained at 36-38 • C using a heat-controlled stage.

2D2C DPIV
DPIV processing of log-scale intensity image data was performed using PIVlab [16]. Contrast limited adaptive histogram equalization was employed as an image preprocessing method for CSF flow images. Window sizes ranged between 32-128 pixels, and were chosen such that each window had ∼5 particles on average. A two-pass algorithm was employed to use half the original window size on a second-pass attempt at peak identification.
DPIV can lead to false matches especially in areas of low back-scattered intensity. Specifically, in our CSF and mouse trachea flow measurements, particles were significantly sparser than in epithelial flow measurements. Because we typically correlated over 20 frames and thus had n = 19 measurements at each spatial location, we were also able to calculate a standard deviation σ and mean µ of our velocity measurements at every spatial location. Thus, for the CSF and tracheal measurements, velocity measurements were discarded at locations where the standard deviation was greater than 100 µm/s and 20 µm/s, respectively.

Synthesis of 3D3C data from 2D2C DPIV
After generating 2D2C data using DPIV, we synthesized 3D3C data using the following method. An OCT B-scan is in a plane that is located in a global xyz coordinate system ( Fig.  2(a)). The angular orientation of any given image plane in that global volume can be described by the angle θ that the image plane forms with the global xz-plane ( Fig. 2(b), top left). In principle, one could measure DPIV data in two orthogonal planes that correspond to the global xzand yz-planes (θ 1 = 0 • and θ 2 = 90 • , respectively). In this framework, the global v x component corresponds to the transverse velocity component (the component orthogonal to the optical axis) measured at θ 1 = 0 • . Likewise, the global v y component corresponds to the transverse velocity component measured at θ 2 = 90 • . Lastly, the global v z component is the average of the axial components measured at each θ . In practice, however, DPIV becomes less reliable when there is a significant out of plane flow in a given image plane (see Appendix for details on SNR tradeoff) [15]. Thus, we chose to measure flow in two non-orthogonal planes that were still primarily directed along the expected direction of flow θ f .
In order to generate 3C data from two non-orthogonal measurements ( Fig. 2(b), top left), we first define each plane with respect to the global x-axis as θ 1 and θ 2 . We define the unit vector associated with these axes as e 1 and e 2 . When we generate a DPIV measurement in plane 1, we estimate an axial (v z ) and lateral (v 1 ) velocity ( Fig. 2(b), bottom). Similarly, in plane 2 we generate an axial (v z ) and lateral (v 2 ) measurement. The key step is to note that the lateral flow measurements v 1 and v 2 in each plane are just the projections of the global velocity vector (v x , v y ) onto the image plane axes e 1 and e 2 . More explicitly: Solving for v x and v y yields: Equations (3) and (4) are valid for the points that two planes intersect. For two imaging planes, there is only one (x 0 , y 0 ) line of data intersection ( Fig. 2(b), bottom). In order to increase the number of intersecting sampling points, we acquired parallel imaging planes with an offset in the xy-plane, typically approximately 10 for each angle, leading to a grid of 100× A-scan depths ( Fig. 2(b), top right). As described above, the global v z estimate was calculated by averaging the measurements of the local v z measurement from each plane.

Directional DLS velocimetry
The methodology for recovering the lateral flow velocity by directional DLS-OCT is described in depth in [11]. Given the time-varying complex-valued interferometric signal E(x 0 , y 0 ,t) at a single spatial location (x 0 , y 0 ) (Figs. 1(c) and 2(c)), we can compute the autocorrelation g(x 0 , y 0 ; τ), where τ is the correlation time lag. The decorrelation rate γ of the signal can be estimated by fitting the autocorrelation to the function |g(τ)| = exp[−(γτ) 2 ]. γ is a measure of the total decorrelation time of the signal, which is proportional to the speed of particles flowing through a focal volume. By additionally imparting a scan bias v scan x either in the direction of or opposite flow, we increase or decrease the effective decorrelation rate γ as measured by our system ( Fig. 1(d)). By measuring γ over a range of scan velocities, we recover a curve that follows the form: We can then fit this decay-rate scan bias curve γ(v scan x ) to solve for the parameter B. The B parameter corresponds to the scan bias speed that best matches the flow rate along the scan axis, and thus is equal to v x . A typical set of scan biases was, for example, ± 0.24, 0.42, 0.56, 0.84, 1.05, and 1.4 mm/s, with 5 repeats at each scan bias.
The effect of the directional DLS-OCT procedure is to give us the projection of the lateral velocity along the plane of our scan bias. By scanning along the perpendicular axis with the same scheme, we can additionally estimate the perpendicular velocity component v y . Similar to DPIV, we can in principle set up our two scan axes to coincide with the global x and y axes. Also similar to DPIV, however, SNR considerations [24], combined with prior knowledge that flow is primarily directed along the global x-axis, led us to rotate our scan axes by 45 • (Fig.  2(c), red X). We measure v 1 and v 2 and similarly reconstruct the global v x and v y measurements using Eqs. (3) and (4).
In each of our measurements, we can estimate the Doppler signal simultaneously, giving an estimate of v z . Overall, then, each set of scans yields v = (v x , v y , v z ) along an A-scan. Repeating this procedure at multiple locations ( Fig. 2(c), gray X's) then yields a 3D3C velocity field estimate.

Three-dimensional volume rendering
In order to optimally display and visualize three-dimensional data, we used two subroutines downloaded from the MATLAB Central File Exchange. The first subroutine was a 3D Canny edge-detection algorithm to perform volumetric image segmentation [25]. Edge detection was performed after simple intensity-based thresholding and allowed us to identify the structural features of the tadpole, ventricle, and tracheal morphology. The second subroutine enabled vector field visualization using a 3D quiver plot [26].

Validation of 3D3C OCT velocimetry using calibrated flow phantoms
We validated 3D3C DPIV-OCT and DLS-OCT by measuring velocity fields in a flow phantom. Using a syringe pump running at a bulk rate of 100 µL/min, we pumped a bead solution through a 0.5 mm (z-direction) x 5 mm (y-direction) rectangular cross-section capillary tube ( Fig. 3(a)). A rectangular tube was chosen to avoid any refraction artifacts due to scanning laterally (y-direction) over the curved surface of a cylindrical tube. Because the phantom crosssection has a large aspect ratio, flow was expected to approximate Poiseuille flow between two planes, characterized by flow directed along the x-direction (v x ) with a parabolic profile varying along the z-axis. We predicted a maximum speed of approximately 0.9 mm/s. The tube was additionally tilted to impart a small z-velocity component.
For DPIV measurements, we acquired a series of 20 B-scans (representative B-scan shown in Fig. 3(b)) at 11 different planes. Each B-scan was 3 mm in length, with 1024 A-scans. The entire acquisition was 450,000 A-scans taken over 16 s. Particles were seeded sparsely enough such that the average number of particles per focal volume was much less than 1. Log-intensity images were cross-correlated in order to yield 2D2C images as described in Section 2.3 and reconstructed into 3D3C measurements as described in Section 2.4. For our DPIV protocol, we had to choose angles θ 1 and θ 2 for our scan plane acquisition. Based on the considerations of our SNR analysis in the Appendix, we chose our imaging planes to be ±15 • oblique to the direction of flow along the x-axis.
For DLS-OCT, particle density seeding was increased such that there were multiple particles per focal volume. A representative M-mode image in Fig. 3(c) shows the intensity trace over time, although we used the complex signal for correlation calculations. For our scan protocol, we scanned at 6 different scan biases from 0.23-2.0 mm/s, with 5 repeats each in a 3x3 grid of xy points, leading to a total of 2,556,000 A-scans over 90 s of acquisition. Note that although both measurements employ correlation analysis, DLS-OCT data was correlated along the temporal axis, while DPIV data was correlated in the spatial domain at a single, fixed time delay corresponding to one inter-frame period (Fig. 1).
As shown in Fig. 3(d) and Visualization 1, the 3D3C reconstruction using DPIV yielded the expected flow profile oriented in the direction of the length of the capillary tube. From the side, oblique, and superior views of the flow field, one can see that the flow resembles planar Poiseuille flow, with the geometric tilt of the tube also apparent. Flow speeds reach approximately 1 mm/s, in good agreement with the expected 0.9 mm/s peak flow. In Fig. 3(e) and Visualization 2 is a 3D3C reconstruction using DLS-OCT with similar views. The recovered velocity field again shows a parabolic flow profile resembling planar Poiseuille flow with a peak of approximately 1 mm/s.
The possible sources of error in this validation measurement include an element of pulsatility (i.e. non-stationarity) in flow due to motor stepping or capacitance in the tube system, small imperfections in the geometry of the rectangular tube leading to non-idealized Poiseuille flow, and a small settling velocity of the beads along the z-axis (estimated <10 µm/s).

Demonstration of directional DLS-OCT in vivo
Because we previously used our directional DLS-OCT protocol to measure flow only in a calibrated phantom profile, we first wanted to confirm that we could extend the measurements towards in vivo flow using our basic 1D2C protocol before moving towards 3D3C acquisition. As such, we tested our methodology by making measurements of v x and v z on epithelial ciliary flow in Xenopus laevis tadpoles.
Measurements were made by allowing the epithelial cilia in Xenopus to drive fluid flow of a physiological solution seeded with polystyrene microspheres (Fig. 4, right). After initial testing, we had to make two changes to our experimental setup in comparison to our original phantom measurements. First, it was qualitatively noted that 100 nm beads, which we had used in our previous studies, seemed to abrogate ciliary flow. In contrast, 500 nm beads did not appear to diminish flow. As such, we selected a 500 nm bead solution, with 0.1% Tween added to prevent bead aggregation. Secondly, it was noted that common path reflections caused a non-overlapping ghost image to appear. We reasoned that the ghost image would also be er- roneously incorporated in our autocorrelation computations. Thus, we used a liquid immersion focal spacer with a tilted air-glass interface in order to avoid a direct reflection (Thorlabs Focal Spacer SD-OCT probe) and to eliminate the ghost image. After making these two changes to our experimental protocol, we then measured epithelial ciliary flow. Figure 4 shows a B-scan of the embryo corresponding to a coronal cross-section of the embryo (see Fig. 5(a) for global orientation). Here, our scan protocol consisted of 5 scan biases between 0.35 -1.4 mm/s, with 5 repeats each. We acquired a total of 230,000 Ascans over 8.2 s. Our reconstructed flow profile along a single line highlighted in red in Fig.  4, is shown on the left. The flow exhibited multiple characteristics consistent with previous measurements of epithelial flow. Flow was directed from head to tail, reaching a peak speed of ∼500 µm/s. Flow was primarily transverse (v x ) in nature, with axial flow (v z ) approximately an order of magnitude lower. Moreover, the measured flow speed was minimal near the glassliquid interface (Fig. 4, blue line), consistent with a no slip boundary from the focal spacer. As we measured deeper into the fluid, the flow speed increased to a maximum near the surface of the embryo, dropping back to zero in the body of the animal itself (Fig. 4, green line).
We found that when using our procedure to obtain velocimetry measurements, in a certain subset of animals, the resulting flow was unphysical in that it did not conform to the qualitative criteria above, such as head-to-tail flow, and maximal speed near the surface. Indeed, DLS-OCT techniques have previously been described to fail in vivo when the model correlation function does not properly fit the data [9]. Additionally, the technique suffers from high error measurement when decorrelation due to diffusivity approaches decorrelation due to advective flow [24]. In our case, we found that measurements failed subjectively in 4/9 animals. We next measured the 3D3C flow field produced by the epithelium of Xenopus laeivs tadpoles using both DPIV and directional DLS-OCT. For our DPIV protocol, we again had to choose angles θ 1 and θ 2 for our scan plane acquisition. As previously noted, flow is primarily head-totail in nature. Thus, we chose our angles θ 1 and θ 2 to be ±15 • with respect to the cranial-caudal axis of the animal (Fig. 5(a)). For each angle we acquired 11 planes, with 20 repeats per plane. Each B-scan was 6 mm in length, with 1024 A-scans per B-scan, for a total of 450,560 A-scans taken over 16 s. For DLS-OCT, we performed 1D3C measurements over a 3x5 grid along the length of the animal. At each point, we acquired 6 scan biases from 0.24-1.4 mm/s, with 5 repeats each, for a total of 1,692,000 A-scans taken over 60 s. Figure 5(b) and Visualization 3 show the results of our DPIV velocimetry measurements. The 3D3C vector field again shows many of the same qualitative features we previously measured in Xenopus epithelial flow. Flow near the surface is primarily directed head to tail. Flow is nearly zero far from the surface of the embryo and reaches a maximum near the surface, consistent with shear-driven fluid flow. Peak flow speeds reach ∼500 µm/s. Figure 5(c) and Visualization 4 show the results of DLS-OCT 3D3C velocimetry. Once again, the epithelial flow exhibits qualitatively similar behavior, with flow primarily directed head to tail and with highest magnitude near the surface of the embryo. Maximal flow speeds in this particular embryo reach nearly 1 mm/s. Of note, as we observed when acquiring 1D2C DLS-OCT measurements, certain data sets resulting from our DLS-OCT protocol did not exhibit the expected features of flow.

Measurement of CSF flow
We next measured cerebrospinal fluid (CSF) flow in NF Stage 46 Xenopus tropicalis tadpoles using DPIV. During development, cerebrospinal fluid circulates through the ventricular system of the brain. OCT-based volumetric measurements ( Fig. 6(a)) provide an intensity-based contrast that can be used to delineate the boundaries of these ventricles within the embryo ( Fig.  6(b)) based on morphological appearance [27]. As previously described in NF Stage 47-49 X. laevis embryos, CSF circulates amongst the various chambers of the ventricles (Fig. 6(c)), from the fourth ventricle and midbrain ventricle, through the cerebral aqueduct, to the third and lateral ventricles [28]. While previous studies of CSF flow were all performed by injecting exogenous contrast agents [29, 30], we observed that endogenous particulate matter provided a native scattering contrast agent in OCT that can be used to perform DPIV (Fig. 6(e) and Visualization 5).
Using this native contrast, we quantified 3D3C CSF flow in the ventricular system of Stage 46 tadpoles. Noting that flow was primarily oriented along the rostral-caudal axis, we used two non-orthogonal planes oriented at ±15 • from the rostral-caudal axis. We acquired 13 planes with 20 repeats each. Each B-scan was 1.2 mm in length and contained 1024 A-scans, for a total of 532,480 A-scans acquired over 20 s.
From the video and subsequent DPIV processing of a single oblique cross-section ( Fig. 6(e) and Visualization 5), one can see that there is recirculatory flow in the ventricles. Our 3D3C reconstruction is shown in Fig. 6(d) and Visualization 6. As previously described [28], flow was directed primarily along the length of the ventricular system, exhibiting a complex recirculatory pattern. In order to better visualize the recirculation, we computationally reconstructed the native axial and sagittal planes of the embryo from our 3D3C data set. The axial cross section shows non-zero vorticity in its plane of motion. Overall, speeds were significantly slower here than epithelial ciliary flow, approximately ∼20-80 µm/s. Notably, due to a sparsity in sampling flow in the third and lateral ventricles, we were unable to reconstruct flow in those chambers.

Measurement of tracheal cilia-driven fluid flow in ex vivo mouse trachea
Using our DPIV protocol, we estimated the 3D3C flow field generated by the ciliated epithelium in an ex vivo trachea from an adolescent mouse. In adolescent mice, the tracheal wall  is sufficiently thin to enable cross-sectional flow imaging of the lumen through the wall itself ( Fig. 7(a,b)). As described in [23], the flow in the excised trachea can be quantified by the imaging the movement of beads in a physiological solution filling the lumen of the trachea (Fig. 7(c) and Visualization 7). Because flow is primarily oriented along the length of the trachea, in a caudal-cranial (tail-head) orientation, we imaged in non-orthogonal planes at ±15 • to the cranial-caudal axis. We acquired 13 planes with 20 repeats each. Each B-scan was 6 mm in length and contained 2048 A-scans for a total of 1,064,960 A-scans acquired over 50 s.
As can be seen in Visualization 7, flow very near to the edges is tail-to-head in nature, but recirculation occurs in the middle of the trachea. We have previously observed that recirculation occurs due to geometric effects. As the ciliated surface moves fluid from head to tail, if the outlet of the trachea is unobstructed, flow remains uniderctional. The presence of an air-fluid meniscus at the cranial aspect of the trachea, however, may effectively block outflow, thus causing fluid to recirculate back down the center of the trachea. As a result, the portions of the flow may appear to be in the opposite (head-to-tail) direction. We typically minimize recirculation in the context of diagnostic tests by ensuring the trachea remains patent. In the context of general flow imaging, however, it may be advantageous to be able to quantify such recirculatory patterns.
The 3D3C reconstruction is shown in Fig. 7(d) and Visualization 8. From the top view, the 3D3C reconstruction again shows flow in a tail-to-head direction near the ciliated surface, and recirculatory in the center of the trachea. Consistent with previous reports, [23], the average flow speed is ∼20 µm/s, with maximum speeds reaching ∼50 µm/s. The primary tail-to-head flow pattern is not as evident on the side view. This lack of easily visualizable head-to-tail flow is possibly due to the fact that head-to-tail flow is localized to a relatively small region near the tracheal surface, as reported previously [23]. DPIV effectively performs spatially averaging when the images are subdivided into smaller regions of interest. Thus the technique may not be as ammenable to highly spatially-localized flow patterns.

OCT-based approaches are useful for quantifying 3D3C flow in ciliated systems
Measurement of cilia-driven fluid flow is important for studying multiple organ systems. Cilia are important for clearance of mucus in the bronchi, circulation of fluid in the ventricles of the brain, establishment of the right-left axis in the embryonic node, and for movement of ova in the Fallopian tubes. Each of these ciliated organ systems has a unique geometry, ranging from tubular (bronchi, Fallopian tubes) to ellipsoid (node) to irregularly shaped (ventricles). Because the direction of flow cannot be inferred from the shape of the enclosing structure, full characterization of cilia-driven fluid flow requires 3D3C velocity estimation.
Here we showed that two OCT-based techniques, DPIV-OCT and DLS-OCT, are capable of quantifying 3D3C flow. Both techniques performed well on a pressure-driven flow phantom, correctly obtaining the direction of flow, and yielding a parabolic profile consistent with planar Poiseuille flow. Additionally, we found that both techniques are capable of measuring biological cilia-driven fluid flow in the Xenopus animal model system. Consistent with previous results from particle tracking velocimetry, we found that flow was primarily directed in a head-tail direction, with highest speeds of 500-1000 µm/s near the epithelial surface and decreasing farther from it.
We also showed that DPIV can be used to quantitatively characterize CSF flow patterns in the brain. Cilia-dependent CSF flow has been characterized previously, and several genes are known to disrupt proper development of flow in the ventricles due to their effects on cilia [30][31][32]. In particular, Xenopus provides an important animal model system for studying CSF flow because the ventricles can be imaged without dissection of the brain [27], and many tools exist that allow rapid and facile genetic manipulation . Previous imaging studies of Xenopus CSF, however required injection of exogenous contrast agents, typically visualized with epifluorescence or brightfield microscopy. In contrast, endogenous particulate matter in the CSF is sufficient to perform OCT-based DPIV and quantify fluid flow. We are currently uncertain of the nature of the scattering material in the developing embryo brain, although 600 nm extracellular membrane particles have been reported in rat embryonic ventricular fluid [34]. The material was present in approximately 50% of the embryos we imaged at NF Stage 46. Further investigation may determine the exact nature of this light-scattering material.
We lastly showed that the use DPIV to quantify the 3D3C flow field in mouse trachea. In the trachea, the airway surface liquid (ASL), which consists of the mucus layer and the periciliary layer, is known to be cleared in a tail-to-head fashion [35]. Because of the typically thin layer size (7 µm periciliary layer and 8-80 µm mucus layer) [1], imaging the ASL directly can be challenging. µOCT, a variant of OCT with higher spatial resolution, has recently been used to quantify ciliary dynamics and tracheal mucus velocity in swine trachea [36]. In other studies, however, the mucus layer is either chemically lysed or washed away prior to quantifying ciliary flow dynamics [37, 38], and we previously used a similar protocol [23]. Here, we showed that DPIV-OCT can be used to quantify the vectorial flow field in mouse trachea. Although in typical experiments we ensure that the trachea is patent and flow is unidirectional, here we were able to observe the primary tail-head flow direction of flow in addition to more complex recirculatory flow patterns arising secondary to outflow blockage. Overall, we believe that this type of assay will be useful for quantifying the effects of physical, chemical, and genetic changes to tracheal ciliary flow. Such studies may include the deleterious effects of hyperoxia in a neonatal setting as well as the ability of pharmacological agents to improve ciliary function [23].

Comparison of DPIV versus DLS-OCT
Based on the underlying physics of signal formation, DPIV and DLS-OCT are both suited for different imaging regimes. DPIV is better suited for slower, less densely scattering fluids, and DLS-OCT for faster, more densely scattering fluids. DPIV relies on image cross-correlation, which in turn is dependent on the same pattern being identified in subsequent frames. In DPIV, pattern decorrelation diminishes SNR (see Appendix for additional details). As a result, fast flows can cause patterns to move through the plane of motion, leading to a decrease in SNR. Additionally, for implementation of DPIV in a coherent imaging modality such as OCT, it is helpful, though not strictly necessary, to use a scattering solution with less than one particle per focal volume on average. In more densely scattering regimes, decorrelation can occur more rapidly due to speckle formation [40], imposing a requirement for higher imaging speed.
In contrast, the decorrelation time of speckle patterns, a limiting factor in DPIV, is actually utilized in DLS-OCT to provide an estimate of speed. DLS-OCT is typically restricted to regimes where there are many particles per focal volume. Additionally, DLS-OCT requires many independent speckle decorrelation events, and thus is optimal at faster flow speeds that lead to more rapid decorrelation. In this sense, the two techniques are optimized for different flow regimes, but DPIV is capable of operating in a wider range of scattering conditions. For these reasons, although both techniques were suitable for relatively fast (∼500 µm/s) epithelial ciliary flow, DPIV was better suited for the sparsely scattering CSF, and the slower tracheal flow speeds.
Both techniques as implemented here have an implicit assumption of temporal stationarity. In both of our techniques, we reconstructed a 3D3C fluid flow field from lower dimensional data. Because the underlying data is not acquired simultaneously, the ability to reconstruct the fluid flow field requires that the fluid is undergoing steady state flow. That is, although the individual ensembles of particles may change, the vectorial flow field must be invariant in time, with dv/dt = 0. In the case of cilia-driven flow, this stationarity condition is typically met [4].
Both techniques also make use of the spatiotemporal correlation function g(x 0 , z 0 ; δ x, δ z, τ). In the case of DPIV, the spatial correlation function at a fixed time lag δt is calculated, while in DLS-OCT, the temporal correlation function is calculated at a fixed δ x = δ z = 0 (i.e. autocorrelation). In both cases, in order to empirically calculate the correlation function, statistical averaging over some domain needs to be performed. In the case of DPIV, averages are calculated over spatial data (i.e. local regional binning), while in DLS-OCT, averages are calculated over temporal data. As a result, DPIV implicitly assumes local spatial stationarity, while DLS-OCT assumes local temporal stationarity.
One consideration regarding the assumption of stationarity is whether the addition of tracer particles affects either the stationarity of flow or the flow field itself. Specifically, we observed that 100 nm beads appeared to affect ciliary flow. Beads have the potential to affect the flow pattern either by interacting with the ciliated surface or by changing the rheological properties of the surrounding fluid. As to the second point, our beads occupied less than 0.1% of the volume, leading to an estimated <0.25% change in viscosity [39]. Regarding the overall effect of particles, we qualitatively observed that the flow pattern was unchanged upon addition of beads of any size ≥500 nm over a period of ∼1-2 min. This observation guided us to use 500 nm beads for our DLS-OCT acquisitions, and 5 µm beads for our DPIV acquisitions. As such, we believe that the addition of tracer particles in our experiments was unlikely to affect the actual flow field or the assumption of stationarity.
Additionally, although both techniques are acquired by the same OCT system with equivalent spatial resolution and sampling, because DPIV uses spatial averaging, with typical window sizes of 32-128, the final spatial resolution of velocity measurements is lower. Analogously, because DLS-OCT uses temporal averaging, the temporal resolution of the technique is somewhat lower than DPIV. In practice, DLS-OCT requires longer data acquisitions to reliably make velocity estimates, approximately 400,000 data points at a single location or 2 seconds of signal integration per scan bias. Thus, we found it impractical to acquire data at more than 15 or so (x 0 , y 0 ) points given our current data processing framework. Because of this overall constraint, our final reconstructed 3D3C data was significantly coarser in spatial sampling than DPIV.
Lastly, as mentioned previously, in a certain subset of data sets acquired using DLS-OCT, the expected dynamics of the autocorrelation function do not appear to be well-modeled using the traditional DLS-OCT autocorrelation function. This phenomenon has been previously reported [9]. As a result, in large subset of data sets, we observe non-physical vectorial field reconstructions, including predictions of no flow where flow is clearly visible from OCT video acquisition. Possible reasons may include violations of temporal stationarity (i.e. the flow field changes over time), a non-Gaussian focus, scanner noise, and propagation of noise (i.e. shot noise) of the original signal.
A summary of some of the comparisons between DPIV-OCT and DLS-OCT is included in Table 1:

Conclusion
We demonstrated that complex three-dimensional, three-component (3D3C) microfluidic ciliadriven fluid flow can be quantified using OCT-based, correlation-based velocimetry methods. We showed that both DLS-OCT and DPIV-based approaches yield reasonable estimates in a flow phantom. Additionally, both techniques can characterize the flow field of epithelial-driven ciliary flow in Xenopus, although DLS-based methods seem to have a high rate of artifact. DPIV can additionally take advantage of endogenous scattering material to quantify fluid flow in the ventricles of the brain, and can also be used for slower tracheal ciliary flow. Overall, we believe

Scattering Regime
Optimized for sparsely scattering media (<1 particle / focal volume), but possible to be used in speckled regime Optimized for densely scattering media ( 1 particle / focal volume) Axial and lateral spatial resolution dictated by ROI in processing (20x20 µm 2 -120x120 µm 2 ) Axial spatial resolution equivalent to imaging system (∼7 µm), lateral resolution limited by extent of lateral scanning (30-50 µm) that OCT-based approaches have great promise for elucidating important ciliary physiology.

SNR Optimization in DPIV
Prior to implementing our DPIV protocol, we investigated the signal-to-noise ratio (SNR) properties of our method to optimize acquisition. SNR in DPIV typically depends on many factors including shot noise of the original acquired images, diffusive motion, and the number of independent measurements averaged [15]. In the context of this work, we focused specifically on the effects of measurement angles θ 1 and θ 2 relative to the flow direction θ f (Fig. 2(b)). As described in Eqs. (3) and (4), two measurements v 1 and v 2 are required to reconstruct two orthogonal flow components v x and v y . The SNR of this reconstruction in turn depends on (1) the SNR of each of the two independent measurements v 1 and v 2 , and (2) the degree of linear independence of the two measurements. As outlined in Fig. 8(a), we will show how both of these quantities depend on the measurement angles, first by calculating the dependence of the error of v 1 and v 2 on measurement angle, and then by propagating that error through equations (3) and (4)  In DPIV, two adjacent frames are used to calculated the spatiotemporal cross-correlation function g(x 0 , y 0 ; δ x, δ z, τ = δt) at a variety of spatial lags δ x and δ z for a fixed temporal lag δt, equal to the interframe time ( Fig. 1(a)). The location of the peak of this function gives the displacement (∆x, ∆z) between the two frames, which is then used to calculate the velocity. It has been previously described that the SNR of velocity measurements in DPIV is proportional to the height of this correlation peak relative to the height of the second highest peak [15]. At low particle concentrations (<1 particle per focal volume on average), OCT-based correlation methods reduce to that of incoherent image formation [41]. Additionally, from image correlation spectroscopy, it is known that the spatiotemporal correlation function of an incoherent image takes the form [42]: where w xy and w z are the beam waist of the Gaussian focus in the xy-plane and z-axis, respectively; D is the diffusivity of the sample; N is the number of particles per focal volume; and v = (v x , v y , v z ) is the flow velocity of the sample. Equation (5) gives the height of the correlation peak. The physical interpretation of Eq. (5) is as follows. There is perfect correlation when an image is compared to itself with zero temporal and spatial lag, i.e. g(δ x = 0, δ z = 0, δt = 0) = 1 because we are multiplying the signal with a duplicate of itself. A maximum in correlation is found between two adjacent frames (δt = 0) when δ x = v x δt and δ z = v z δt, as these terms correspond to in-plane translation. In this case, however, the correlation may no longer be unity because of particle diffusion and out-of-plane (v y ) motion. In other words, diffusion and out-ofplane motion cause image patterns to be lost instead of translated purely in-plane, thus leading to a decrease in the magnitude of the correlation.
We know that the SNR is proportional to the height of this correlation (assuming that the noise floor is constant) but also depends on other factors including signal intensity. We introduce a pre-factor α to encapsulate other non-angle dependent factors. Additionally, we neglect the diffusive contribution 1/(1 + 4Dτ/w 2 xy ) and 1/(1 + 4Dτ/w 2 z ) in the exponentials, noting that based on the diffusivity of 5 µm beads and an interframe time of 80 ms, this term is 0.999, or approximately 1. Lastly, based on the angle of flow θ f , the angle of the plane of imaging θ i (i = 1, 2), and the total speed v = (v 2 x + v 2 y + v 2 z ) 1/2 , we define v y = v sin(θ i − θ f ). These definitions leave us with an expression for SNR: Equation (6) now relates how SNR in 2D2C DPIV depends on the relative flow angle. Plotted in Fig. 8(b) is Eq. (6) for two values of flow as a function of scan angle, with the flow angle arbitrarily set to be 90 • . The SNR has a peak both where θ 1 = θ f and θ 1 = θ f +180 • . These two cases occur when flow has only in-plane components. Equation (6) reaches a minimum where the imaging plane is 90 • (perpendicular) to flow. The degree to which SNR is decreases from out-of-plane flow depends on the quantity vδt/w xy , or the magnitude of out-of-plane motion.
To confirm these predictions, we performed DPIV on our parabolic flow phantom and varied the angle of the imaging plane with respect to the flow angle, defined to be 90 • . We calculated an empirical SNR as follows. We recorded 20 frames, leading to n = 19 independent measurements of v x at each location. SNR was calculated at a point by point basis in the channel as µ v x /σ v x , where µ v x and σ v x were the mean and standard deviation, respectively, of the 19 measurements at a single location. The SNR was averaged over all points in the channel to come up with a single SNR . Shown in Fig. 8(c) is the SNR(θ ) curve, which shows the same qualitative behavior as predicted by Eq. (6), reaching a peak near 0 • where flow is nearly coincident with the imaging plane, and dropping sharply, and increasing again at 180 • . Notably, the SNR points on Fig. 8(c) have been scaled to reach a maximum of unity in order to better compare with Fig.  8(b).
Once we know the SNR of each individual measurement, the next step is to predict how that error due to our measurements of v 1 and v 2 , which we denote as δ v 1 and δ v 2 , propagates to our estimation of v x and v y . We denote the errors of our global velocity measurements as δ v x and δ v y . Assuming relatively small error, this step can be predicted based on standard error propagation formulas: Applying Eq. (6) to Eqs. (7) and (8), we derive the formula: where θ 1 and θ 2 are the angles of the imaging planes. The dependence of δ v 1 and δ v 2 on flow angle and plane angle has been highlighted by writing them as a function of θ f and θ i . Importantly, the error in each of these terms diverges when θ 1 = θ 2 . This large increase in error is a statement of the fact that in order to measure two orthogonal velocity components, two linearly independent measurements must be made. Measurements that are linearly independent but non-orthogonal will give a unique answer, but the degree of error propagation depends on how orthogonal the two angles are. In order to optimize our acquisition, we optimize not δ v x or δ v y alone, but rather the sum of the two. Shown in Fig. 8(d) are contour plots of the relative SNR v xy , defined as 1/(δ v 2 x + δ v 2 y ) as a function θ 1 , θ 2 for low and high speed, with θ f now set to be 0 • . SNR is very low when θ 1 = θ 2 or θ 1 = θ 2 + 180 • due to the linear dependence of the two measurements. For low flow speeds, maximal SNR occurs near θ 1 = ±45 • , θ 2 = −θ 1 . (θ i and θ i + 180 • are equivalent, because the imaging plane is coincident with flow, and thus there are 8 total equivalent solutions that maximize SNR.) At higher flow speeds, because SNR of v 1 and v 2 drops significantly when θ 1 and θ 2 are not coincident with θ f , it becomes more advantageous to align each independent plane closer with the direction of flow, even if the two planes are no longer orthogonal. In the case of high flow, then, SNR is still maximized when θ 1 = −θ 2 , but now at some angle that is less than 45 • .
The maximum at each of these contour plots can be calculated numerically, and the location of the maxima can be plotted as a function of flow speed. Figure 8(e) shows these results. As the magnitude of speed approaches zero, the optimal imaging angle approaches ±45 • , while when magnitude of speed is high, the optimal imaging angle moves closer towards being coincident with the direction of flow.
The interpretation of Fig. 8(e) is as follows. At low flow speeds, the absolute error δ v 1 and δ v 2 are relatively low, so these two measurements should be orthogonal in order to minimize error propagation from v 1 and v 2 to v x and v y . As such, the planes are placed optimally at ±45 • oblique to the direction of flow. Even at these slow speeds, however, the SNR is less than unity because both imaging planes cannot be placed coincident with flow simultaneously. As the flow speed is increased the absolute value of δ v 1 and δ v 2 begins to increase. As that error increases, the gains from having two orthogonal measurements starts to be outweighed by the fact that the absolute error in v 1 and v 2 can be reduced by bringing both planes closer to coincident with the angle of flow. In the most extreme scenario, when flow is very fast, both planes must be nearly coincident to it to get any reasonable value of v 1 and v 2 . As v 1 and v 2 begin to become more collinear, however, small errors propagate more strongly. Thus, at higher flow speeds, even though we optimize SNR by bringing the two imaging planes coincident with flow, the total SNR of the process is still much worse than when we can make completely orthogonal measurements in the low flow speed regime.